The computational design of catalytic materials is a high dimensional structure optimization problem that is limited by the bottleneck of expensive quantum computation tools. Current implementations of first principles computational models for catalyst design are datahungry, problemspecific and confirmatory in nature. However, they can be made less datadependent, more transferable and exploratory by developing both forward and inverse catalyst mapping tools that are either inexpensive correlations, like scaling relations, or regression models that are based on relevant descriptors analysis. This work reviews the current application and the possible landscape for future advancements of such tools for developing generalized schemes for catalyst design and optimization.
Computational design of catalytic materials is a high dimensional structure optimization problem that is limited by the bottleneck of expensive quantum computation tools. An illustration of interaction of different factors involved in the design and optimization of a catalyst.
Today heterogeneous catalysis facilitates highly energyefficient selective molecular transformations for over 90% of the chemical manufacturing processes and contributes towards 20% of all industrial products.^{1,2} Given the importance of heterogeneous catalysts in chemical processes, developing rational heterogeneous catalyst design rules is one of the most fundamental goals in reaction engineering and is the key towards developing future sustainable chemical technologies. However, catalyst design is a complex process, involving numerous interacting factors, such as active species (metal, oxide,
With the development of computeraided data storage and processing techniques, vast amounts of accumulated knowledge on catalysts can be employed in analysis simultaneously.^{3} Additionally, the multidimensionality of the catalyst design problem can be represented in a machine readable format for developing prediction models. For example, ML algorithms have been trained on data corresponding to catalyst performance metrics (
Developing a generalized design scheme for heterogeneous gasphase reactions that is computationally effective, accurate, and transferable (valid for different reaction systems and catalyst compositions), without relying on dataheavy techniques, is among the major research goals in catalysis. In this review, we discuss the implementation of popular theoretical tools, approaches and workflows used in materials investigation for catalyst screening and optimisation under the listed constraints. We compare their applicability for different problem systems while highlighting the recent advancements made in order to bypass the current computational bottleneck in catalyst design. We start with the tools used for computing reaction energetics over catalysts surfaces. These tools include firstprinciples methods, their semiempirical/empirical adaptations and molecular dynamics. We compare different methods reported in the literature for developing reaction networks, catalyst optimization and screening. Lastly, we summarize stateoftheart in theoretical catalyst design and discuss its potential future landscape based on the existing developments.
The kinetics on a given catalyst surface are driven by the energetics of elementary steps occurring on active sites. In this section, we discuss the most common approaches for computating reaction energetics (adsorption energies, activation barriers) over catalyst surfaces. Density functional theory (DFT) is the most popular choice for quantum mechanical computations. Other tools are essentially semiempirical/empirical implementations derived from or trained on DFT but are relatively less expensive. This also includes tools for molecular dynamics. We further compare them based on three key factors: transferability, accuracy, efficiency, applicable system size, and highlight their respective limitations, see
In the realm of reaction chemistry computation, quantum mechanical calculations have emerged as indispensable tools.^{14–17} These calculations have revolutionized our understanding of the behavior of atoms, molecules and materials by providing insights into their electronic structure and properties at the microscopic level. Hartree–Fock (HF) calculations are one of the oldest quantum mechanical methods based on wavefunctions that is used to study electronic structure.^{18–20} They are based on a selfconsistent field approach, where the electrons are treated as noninteracting particles within the mean field generated by the other electrons. While HF calculations provide a good starting point for understanding the electronic structure, they have limitations, particularly in accurately describing electron correlation effects.^{19} Several postHartree–Fock methods have been developed to account for electron correlation effects beyond the HF level. Examples include configuration interaction (CI),^{21} coupled cluster (CC),^{22} Møller–Plesset perturbation theory (MP2),^{23} and density functional theory (DFT).^{24,25} These methods offer improved accuracy by including electron correlation contributions, but some of them can be computationally demanding, especially for large systems. Among the postHartree–Fock methods, density functional theory (DFT) is the only method based on density function instead of wave functions. It has become a widely employed computational tool for studying the electronic structure and reaction energetics, especially for catalytic systems.^{26} It offers a versatile framework to investigate a wide range of systems and phenomena.
DFT is a widely used method for electronic structure calculations. At the core of DFT calculations is the concept of electron density. Instead of solving Schrödinger equation for the system's wave function, DFT focuses on determining the electron density distribution, which is a more tractable quantity.^{24,25} This approach has proven to be highly effective in studying a wide range of catalytic systems.^{26} The most widely used form of DFT is Kohn–Sham density functional theory (KSDFT). It was proposed by Walter Kohn and Pierre Hohenberg in 1964 and has become the cornerstone of modern DFT calculations.^{24,25}
The Kohn–Sham approach introduces a set of auxiliary noninteracting electrons, represented by Kohn–Sham orbitals, to approximate the real system's electron density.^{24,25} These orbitals are derived from a fictitious system in which the electron–electron interactions are neglected. However, the electron density obtained from the Kohn–Sham orbitals should match the electron density of the real system. In KSDFT, the electronic structure problem is mapped onto a set of noninteracting electrons moving in an effective potential. To solve the Kohn–Sham equations, a basis set is used to expand the Kohn–Sham orbitals. The choice of the basis set is critical in accurately representing the electronic structure of the system.^{27} The basis set represents a set of functions that span the space in which the electronic wave functions are defined. In KSDFT calculations, the choice of a basis set and the package used to solve the equations depends on the nature of the system being studied: cluster/isolated system, bulk (periodic) system, or condensed media.
In cluster calculations, where a small group of atoms or molecules is isolated from its surroundings, the basis set is typically chosen to describe the electronic structure of the cluster accurately.^{28} The commonly used basis sets in cluster calculations include Gaussiantype orbitals (GTOs).^{29,30} Gaussian,^{31} ORCA,^{32} and GAMESS^{33} are widely used software packages for calculations involving cluster/isolated molecules.
In bulk or periodic system calculations, the electronic structure is determined for extended periodic structures such as crystals and surfaces. In these calculations, the basis set is expanded to include periodic boundary conditions.^{34–36} Planewave basis sets are commonly used in periodic system calculations, along with pseudopotentials or projectoraugmented wave (PAW) potentials to efficiently treat the electron–ion interactions.^{36,37} Since the potential is periodic, the plane wave basis sets are wellsuited for the description of periodicity in the crystal lattice. VASP,^{34,38} Quantum Espresso^{39,40} and CASTEP^{41} are popular packages for periodic DFT calculations.
Condensed media calculations involve studying systems where the electronic structure is influenced by the surrounding environment, such as liquids, solutions, or solids with embedded solvents. CPMD, which stands for Car–Parrinello molecular dynamics, is particularly useful for simulating condensed phase systems and exploring reaction mechanisms with atomistic detail.^{42} CPMD is a computational method that combines molecular dynamics simulations with DFT. CPMD utilizes plane wave basis sets to describe the electronic structure and employs the Born–Oppenheimer approximation to treat the nuclear motion.^{42} This approach allows for the study of dynamic processes, such as chemical reactions and materials transformations, at the quantum mechanical level. In short, the choice of basis set in DFT calculations significantly impacts the accuracy and computational cost of the calculations. Larger and more sophisticated basis sets can provide a more accurate description of the electronic structure but at the expense of increased computational resources and time.
In addition to the selection of an appropriate basis set, one of the other major challenges in DFT is choosing an appropriate exchange–correlation functional, which accounts for the electron–electron interactions.^{43,44} However, the exact form of the exchange–correlation functional is not known. Many approximations exist for the exchange–correlation functional, and selecting an appropriate functional is crucial to achieve reliable results. Generalized gradient approximation (GGA)^{45} functionals, such as PBE,^{46} PW91,^{47} and revPBE,^{48} are commonly employed in DFT calculations. Benchmarking of functionals is necessary to ensure their reliability in DFTbased catalyst design procedures. For instance, many studies employ GGAbased PBE and PW91 functionals for studying CO_{2} conversion reactions on solid catalysts. However, inconsistencies arise between the calculated values using these functionals and experimental data.^{49–51} Notably, the binding energies of CO_{2} on metal surfaces, often deviate significantly.^{52} Additionally, GGA functionals are inadequate for describing weakly bound systems, where dispersion interactions play a crucial role.^{53–55} Incorporating van der Waals (vdW) interactions through specific correlation functionals can improve the accuracy of dispersionbonded systems.
The Hubbard U is yet another DFT parameter based on the Hubbard model (DFT+U method)^{56,57} that is crucial for accurately describing the electronic properties of transition metal oxides by employing DFT. Metal oxides, with their transition metal elements, exhibit strong electron–electron interactions and localized electron behaviour. Standard DFT functionals may not adequately capture these effects, leading to inaccurate predictions. The Hubbard U term represents the onsite Coulomb repulsion between electrons and helps describe the physical behaviour of transition metal oxides. Incorporating the Hubbard U parameter improves the agreement between DFT calculations and experimental observations, enabling accurate predictions of structural, electronic, and magnetic properties. Determining the appropriate Hubbard U value is challenging and often relies on empirical or theoretical approaches. By including Hubbard U, DFT calculations provide valuable insights into metal oxide properties, including bandgaps, energy levels, charge localization, and magnetic behaviour.^{58,59} The Hubbard U parameter is essential for understanding metal oxides' electronic structure and behaviour within the DFT framework. For further detailed overviews on DFT, its parameters and application we direct the readers to review articles in ref.
Although DFT calculations are accurate and reliable for studying different reactions on a small set of catalysts, they cannot be directly applied to problems that require high throughput computations, particularly in catalyst design problems involving the calculations of potential energy surfaces (PES) for reaction chemistry. Although the computational cost of determining minima in the PES is reasonable, determining the maxima (saddle points) can increase the computation time significantly depending on the number of images used. For complex systems like nanoclusters, single atom alloys (SAAs) and singleatom catalysts (SACs), this increase can be even higher.
Given that the method is based on first principles, it is transferable, but the accuracy is highly affected by the choice of functional and other computational flags. The speed/efficiency of the method is also way below the requirement for a comprehensive analysis. This can be observed from
Type of DFT calculation  No. of calculations  No. of cores  CPU time  Total core hours 

Geometry optimization  600  64  5  192 000 
Transition state searches (minimum energy path)  135  96  15  194 400 
Transition state optimization  45  64  2  5760 
Transition state confirmation 
45  64  10  28 800 
Adsorbate vibrational frequency  24  64  10  15 360 
Grand total of core hours  436 320  
Number of catalyst surface investigated (Ni(111), Ru(0001), NiB(111))  3  
Total CPU (core years)  145 
The density functionalbased tight binding method or DFTB was initially based on a secondorder expansion of the DFT total energy with respect to charge density fluctuations.^{62} However, the most common implementation, DFTB_{3},^{63} includes a thirdorder expansion of the DFT total energy.
DFTB methods can be twotothree orders of magnitude faster than
However, without proper benchmarking, the accuracy of the method can be heavily compromised,^{69} making the use of DFTB tradeoff between speed
UBIQEP is a generalization of the BOCMP (bond order conservationMorse potential) method that is used for modelling chemisorption energetics and reaction mechanisms on metal surfaces. It has a fast and easy computational implementation and the UBIQEP projections of reaction energetics are usually more accurate. It should be noted that the UBIQEP modelling in particular is not competitive but complementary to quantum mechanical modelling.^{71} As input parameters, the method employs thermodynamic observables such as gasphase bond energies and atomic chemisorption energies mostly obtained from DFT. Its output is the surface reaction energetics for all elementary steps in a mechanism thus improving the overall efficiency. Although these energies cannot be transferred to new adsorption systems without the initial DFT input.
The BOCMP and the UBIQEP methods have been applied successfully to analyse mechanisms of many reactions of practical importance such as methanol synthesis,^{72} Fischer–Tropsch synthesis,^{73} and methane reforming chemistry on metal surfaces.^{74} There are also UBIQEP analytical formalisms for bimetallic surfaces however, similar developments cannot be found for systems like single atom catalysts and nanoclusters.
The reactive forcefields (ReaxFF)^{75} are based on a bondorder formalism that implicitly describes chemical bonding without expensive quantum mechanical calculations. Its bondorder parameters are derived from computationally intensive DFT derived methods. Once these parameters are known, computing the energy matrix could be more than 100 times faster than DFT for gasphase heterogeneous reaction systems. ReaxFF potential parameters have already been reported for hydrocarbons chemistry,^{76} H_{2} dissociation on transition metals,^{77,78} carbon interactions with transition metals^{79} and lastly hydrocarbon reactions catalysed on transition metals such as nickel^{80} and vanadium oxides.^{81,82}
In short, ReaxFF can model reactions at the gas–solid interface and assess the stability of SACs,^{83} nanoclusters^{84} and SAAs,^{85} which is pertinent to our problem of heterogeneous catalysts design. Moreover, the available parameters for elements in the literature are also transferable, as long as the aqueous phase reactions are not involved.^{75} However, proper benchmarking needs to be performed to determine the accuracy prior to any calculations and ReaxFF cannot be transferred to new systems without reparameterizing the method.
The major drawback accompanying these forcefield methods is the low accuracy of predictions which can be limited by proper parameterization of interatomic potentials. Recently there have been some developments in this area with the use of machinelearned (ML) force fields. These MLcorrected force fields have been found to be much more accurate. Further details on ML force fields can be found in section 5.1.4.
Method  Transferability  Accuracy  Efficiency  Approx. system size (# atoms) 

DFT  Based on first principles hence more transferable than the semiempirical/empirical methods  Relative comparison of catalysts can be performed accurately  Computational cost increases cubically with the increasing number of atoms in the system  10^{2} 
Absolute catalyst predictions are also possible if the exchange correlation functionals are correctly identified using experimental data  
DFTB  Limited for heavy metal systems  Benchmarking is required to compute the accuracy  Computational cost increases cubically with the increasing number of atoms. However each single point computation is 2–3 three orders of magnitude faster than DFT  10^{3} 
UBIQEP  Needs recalibration in case of a new adsorption system ( 
Accuracy depends on the DFT data used for calibration  Computational cost increases linearly with the increasing number of atoms  10^{2} 
Implementations limited to mono and bimetallic catalyst systems  
ReaxFF  The ReaxFF potentials need to be evaluated and validated every time a new element and/or is added to the adsorption system  Benchmarking is required to compute the accuracy  Computational cost increases linearly with the increasing number of atoms  10^{5}–10^{6} 
Not suitable for aqueous phase reactions 
A reaction mechanism describes a network of elementary reactions (reaction network) which corresponds to the reaction coordinate from reactants to products and byproducts. The hypotheses about possible reaction intermediates and elementary reaction steps are based on experimental observations, literaturebased dissociation and association routes, and/or autoreaction generators algorithms^{89} that employ reaction energy data,^{92–94} templates^{88,95} or heuristics. The methods for developing a reaction network can be broadly classified as:
(1) Automated reaction mechanism generation that is based on reaction rates, and.
(2) Userdefined reaction rules for mechanism generation that are based on literature and experiments.
This section discusses the tools for developing reaction mechanisms
An automatic reactionmechanism generation (ARMG) is based on an algorithm that can perform the following tasks:^{89}
1. Recognize when two or more species in the mechanism are equivalent.
2. Predict all the possible elementary reactions for each species and pair of species.
3. Determine which of these possible reactions are important.
4. Estimate accurately all the necessary thermodynamic and kinetic parameters.
5. Ensure that the mechanism is thermodynamically consistent.
ARMG can be simultaneously used to update the reaction mechanism with catalyst surfaces during the optimisation of overall catalyst performance. This is particularly useful since reaction mechanism identification is based on the surface energetics over a catalyst. Therefore, surface energetics and by extension, the mechanism can vary from one catalyst to another.
A challenge with ARMG is determining the kinetically relevant reaction steps in the initial reaction network. A prior assumption about the list of intermediates and reactions can lead to biased results or limited solutions. If all possible intermediates and reactions are included, then the resulting reaction network quickly becomes unsolvable and, therefore, useless.
Gao
Reaction networks can also be constructed using userdefined rules around a mechanism that is based on experimental evidence and literature. For example, Rangarajan
This method can be categorized as reaction specific, partly intuitive and based on heuristics observed in the literature/experiments. The common framework is to look in the literature for elementary steps corresponding to reactant dissociation, product association and redox reaction steps in similar reactions with known mechanisms. For example, a study on methane reforming^{99} copies the mechanism for dry reforming of methane (DRM) directly from CH_{4} steam reforming with the addition of CO_{2} dissociation steps. It is found that for small molecules reaction systems like DRM, the literature and automated mechanism generators agree on the same reaction network.^{89,99} There are also approaches in literature that automate the rule based reaction network generation^{100} by including a language compiler that can convert an Englishlike reaction language into internal representations and instructions.
Designing a simple, yet representative reaction network is important for limiting the cost of evaluation and ensuring model solvability during mechanism identification and subsequent catalyst optimization. A way to achieve that would be
A reaction network and corresponding energetics observed over a catalyst surface are used for building kinetic models for evaluation of catalyst' performance and their subsequent optimization. There are topdown and bottom up approaches for building a kinetic model. The topdown kinetic model includes approaches based on power law expressions and Langmuir–Hinshelwood–Hougen–Watson (LHHW) mechanism. The assumptions in these approaches are more implicit. For example, assuming certain ratedetermining elementary reactions, quasiequilibrated elementary reactions, and the most abundant surface intermediates.^{102} Bottomup approaches are theoretical models that do not make such initial assumptions about the nature of ratedetermining elementary reactions but can be designed to enable different simplifications as per user requirements. This section discusses the theoretical bottomup approaches used for quantifying catalyst performance.
Microkinetic models (MKM) are a vital tool in catalyst research and design.^{102} Microkinetic models predict and compare the performances of different catalyst surfaces for a given reaction in terms of their reactivity, selectivity and stability against deactivation mechanisms like coking. MKMs are often designed and validated for a specific reaction system based on certain reaction conditions. As a result, complete MKMs are not adopted in subsequent studies, except for maybe certain model fragments, such as reaction networks and energy data. Different assumptions while creating an MKM also lead to different model predictions and the validity of these assumptions can vary. For example, using collision theory instead of transition state theory to model gasphase adsorption, kinetic Monte Carlo instead of meanfield approximation to include the effect of surface coverage, and assuming quasiequilibrium of certain reaction steps.
Regardless of the assumptions, formulating MKMs are computationally expensive. Subsequent catalyst design using MKM is therefore described as an expensive headon approach for performance optimization.
(1) An initial catalyst surface is fed into the microkinetic model.
(2) Reaction energetics are computed based on thermodynamics.
(3) The reaction kinetics are written followed by the reactor massbalance equations.
(4) The model is then solved to predict the corresponding catalyst performance.
The MKMaided catalyst optimization employs the model to develop a function mapping between the catalyst structure and its performance. It then maximises the performance by closely studying the inverse function and its derivative,^{103}
The literature also reports another study^{104} employing the MKMaided catalyst optimization framework but in a less expensive inverse approach. The reaction energies are optimized independently, that is without subsequent optimisation of catalyst structure (see
The reaction energetics of a given reaction system can be very complex with interlinking energy values. A descriptor identification is important to make catalyst screening feasible. Herein, the descriptor would correspond to the energetics of the most relevant elementary reaction step (reader is directed to section 4.3 for further details).
Contrary to the MKMaided optimization approach, the quantitative structure–activity/property relationships (QSAR/QSPR) are based on developing a direct mapping between the microscopic (
Despite the straightforward implementations of QSAR/QSPR, they are often limited to a specific problem system due to the lack of a diverse database, limited data points and reproducibility issues during experiments. A robust and reliable property prediction model is not possible if either of the following is true:
(1) The experimental measurements have high uncertainties.
(2) The chemical diversity beyond training sets.
(3) The range of measured property values is too small.
The descriptorbased catalyst search is the most common approach towards theoretical catalyst optimization. It reduces the cost of exact model evaluation to descriptor evaluation and catalysts with descriptor values corresponding to maximum performance can be screened using the volcano plots^{106} (reader is directed to section 4.3.3 for further elaboration on volcano plots). The overall optimization will depend on the chosen catalyst descriptors and is limited by the volcano plot relation, thus identifying relevant catalyst descriptors is a very important step for catalyst design. We discuss methods used to identify reaction descriptors,
In this approach, the sensitivity of the catalyst's performance w.r.t. individual reaction steps/energies are evaluated (see
Although DRC is an extraordinarily useful concept for reaction mechanisms analysis, the energies of intermediates and transition states are linked through scaling relations (see section 5.2) and so are the rate constants. Thus, they cannot be changed independently during catalyst design. Therefore, instead of employing reaction energies with the highest DRC as catalyst descriptors, another metric is applied (
Based on the DCC, the activation/binding energy of the most sensitive reaction step/species is then chosen as the reaction descriptor. This descriptor is specific to the reaction investigated and can change if the reaction conditions are drastically changed.
Reaction network and pathway analysis is a graphbased approach that is frequently used to identify highly interconnected intermediates and reaction steps that are solely connecting different reaction chemistries. Given a list of elementary reaction steps, reaction network construction of complex reactions outlines all available pathways for product formation. Further ratebased flux analysis on these pathways can be performed using algorithms, like Dijkstra's and variants,^{109} to identify dominant pathways, respective contributions and ratedetermining steps,
Volcano plots are based on Sabatier's principle,^{110} which states that a catalyst should bind a substrate neither too strongly, else it will destabilize the catalyst, nor too weakly, as that will reduce the activity. Thus, the maximum performance can be identified somewhere in between. The volcano plots provide an estimate of catalytic performance in terms of turnover frequency, product yield or over potential, w.r.t. a descriptor variable,
Despite its popularity, volcano plots limit the maximum attainable catalytic performance corresponding to bulk catalytic structures and cannot be generalized to structures like single atom alloys (SAAs) that are known for breaking the scaling relations based on bulk materials.^{111,112} They are also prone to error depending upon the linear scaling relation (see section 5.2) used in their construction. Section 5.2 highlights a specific scenarios where a deviation from volcano plots and scaling relations in observed. The reader is referred to ref.
The ability to solve either the forward or the reverse optimization problems essentially boils down to the catalyst screening method. An ideal catalyst screening tool should be capable of performing highthroughput computation for a range of catalyst structures and composition variations.
In this section, we discuss the most common methods for inexpensive highthroughput screening of catalysts based on identified reaction descriptors,
Machine learning models can be employed to reduce the cost of reaction thermochemistry computation, thus enabling efficient catalyst screening. They can be trained on energies evaluated from DFT or other semiempirical/empirical tools depending on the required accuracy and available data. ML models can be implemented at different stages in the overall computation. Implementations have been identified in the literature for predicting DFT wave functions or electronic density,^{113} predicting DFT energies^{114} and accelerating transition state searches.^{115,116} Further implementations include feature construction from relevant descriptors for initial catalyst screening^{117} and reactive forcefields parameters optimization for a more rigorous screening (see section 5.1.4).^{118}
Given the inherent complexity of catalytic reactions, these machine learningbased implementations are relevant tools to efficiently navigate through this highdimensional catalyst optimization problem. Although their performances rely on the training algorithm and feature sets, their predictions are not transferable beyond the training set. In the following sections, we list common catalyst features, available catalyst databases and learning algorithms used for developing these models.
Heterogeneous catalysts for gas–solid reactions are exclusively inorganic materials derived from transition metals. Thus, relevant features for such materials can be identified based on the band theory of chemisorption. However, the following points should be considered when deciding on features for catalyst surfaces:^{119}
(1) The features must be unique in representing the electronic and geometric structures of an active site.
(2) They must be easily computed or readily available from databases to enable rapid screening.
(3) They should be physically intuitive to ensure model robustness and direct inference of chemical insights.
Class  Name  Abbreviations  References 

Stoichiometry 

‖ 

Atomic properties of components  Atomic number 


Atomic weight 



Group 



Period 



Mendeleev number  MN 


Atomic radius 



Covalent radius 



van der Waals radius 

—  
Pauling electronegativity  PE 


Ionization potential  IP 


Electron affinity  EA 


# s valence electrons 



# p valence electrons 



# d valence electrons 



# f valence electrons 



Total # valance electrons 



Magnetic moment/atom at 0 K ground state 



Bulk properties of components  fcc nearest neighbour distance  bulk_{nnd} 

Partial radial distribution function 



Radius of dorbitals 



Radius of porbitals 



Coupling matrix element squared 



Specific volume 



Band gap energy of 0 K ground state 



Space group number of 0 K ground state  — 


Melting point 



Boiling point 



Enthalpy of fusion  Δ 


Surface  Work function 


Cohesive energy 



Surface energy 



Site  Number of atoms in ensemble  site_{no} 

Coordination number  CN, GCN 


Orbital wise coordination number  CN^{α} ( 


Bondenergyintegrated orbital wise coordination number 



Nearest neighbour distance  site_{nnd} 


dband centre 



pband centre 



dbandwidth 



dband skewness 



dband kurtosis 



dband filling 



spband filling 



Antibonding states (e_{g}) filling 



Density of dstates at Fermi level  DOS_{d} 


Upper dband edge 



Density of spstates at Fermi level  DOS_{sp} 


Thermodynamic stability of active sites  BE_{M} 


Crude estimate of the property being predicted  — 

GCN corresponds to the generalized coordination number.
These primary features can be implemented directly or can be used tohandcraft secondary features using dimensionality reduction methods like LASSO, SISSO,
From
Even so, these expensive
Noh
Most of the studies employing MLbased reaction energy prediction models use inhouse developed catalyst databases specific to the problem. This is because universal catalyst databases are difficult to realize given the enormous compositional and structural variations that are feasible. Even so, there are a few openaccess databases that have been commonly employed across different studies and have consistently been improved to include more compositional and structural variations, see
The CatApp database^{140} includes reaction energies for approximately 3000 surface reactions on closepacked as well as stepped fcc and hcp on the following metal surfaces; Ag, Au, Co, Cu, Fe, Ir, Mo, Ni, Pd, Pt, Re, Rh, Ru, Sc, V. All values have been calculated with the same code (DACAPO), the same exchange–correlation energy functional (GGARPBE^{46}). Therefore, one adsorption energy or reaction barrier can be compared to another with some confidence. The limitation of CatApp is that it does not store the atomic structures that are an output of the electronic structure calculations. Because of this data reproducibility is limited.
The CatalystHub database^{141} contains more than 100 000 electronic structure geometries and results from more than 50 publications, including the ones present in the CatApp database. The datasets stem includes a direct link to their respective publications. The database contains a large variety of alloy surfaces and oxides. A large part of these reaction energies stems from the highthroughput study of chemical adsorption and hydrogenation on more than 2000 bimetallic alloy and pure metal surfaces.^{142} The calculations are performed using different tools and functionals and, hence, cannot be directly compared. However, the geometries could easily be used to recompute the energies as per requirement.
The Open Catalyst database (OC20)^{143} is the most recent and comprehensive database that is inclusive of previous databases and has further calculations based on the Materials Project database^{144} that are consistent in their computation parameters like the DFT functional. OC20 consists of 1 281 121 density functional theory (DFT) relaxation calculations for 5243 different material compositions and 82 different adsorbates (small adsorbates, C1/C2 compounds, and N/Ocontaining intermediates). Although the adsorption energies most likely do not correspond to the lowest energy configuration.
It should be noted that despite the knowledge accumulation over the years for creating heterogeneous databases, they are far from being representative of the vast number of possible catalytic systems. Even the most recent OC20 database seems to include only 18.9% of the total systems and only 0.07% of the possible calculations.^{143} Databases are also biased in the sense that they mostly include systems that have been considered hot by the computational heterogeneous catalysis community.
Materials Project^{145} is another notable database that was constructed to accelerate material property investigation in general but also guides catalyst structure investigation. It has structural and
In addition to these databases, there are also several studies^{114,127,137,146} that have published data on intermediate binding energies, mostly CO* and H* adsorption energies. Although these data are mostly specific to their corresponding catalytic systems under consideration.
The learning algorithms for ML prediction models are selected based on the training set distribution and size. Training set size is often the limiting factor given the expensive data point evaluations. Most of the studies employ problemcentric databases where compositional variations are very limited, and catalysts are structurally similar with defined binding sites. Catalyst features, as shown in
There are structural variations within the same catalyst composition and stable binding sites could also be different on the same surface depending on the adsorbate. To incorporate these variations, graphbased representations of catalyst surfaces have been developed.^{8,150} These representations do not require extensively handcrafted features (
Convolution neural networks (CCNs) and variants^{151} are the most common approaches for preserving this information along the different training layers of neural networks and are found to outperform other implementations.^{150} However, these models have a huge data requirement, more than 10^{4} data points for a single adsorbate. These limitations bring us to the different training algorithms that have been specifically developed for these scenarios.
Active learning is a userinteractive learning technique that is applied in cases where obtaining data labels is expensive. The algorithm identifies new training data dynamically based on a query strategy. The technique is helpful in the effective utilization of resources for data point evaluation and model training. It is found to increase the model prediction RMSE of binding energies from 0.18 eV to 0.05 eV and from 0.65 eV to 0.4 eV in bimetallic systems^{127} and perovskites respectively.^{7} Active learning algorithms have also been employed on datasets consisting of inhouse generated data coupled with the abovementioned available databases.^{127,146} Although there is a drawback to this supervised learning technique,
Transfer learning is a popular technique where a model trained on one task is repurposed on a second task,
MLbased force fields are one of the many applications of ML in heterogeneous catalysis for evaluating free energy values and generating reaction mechanisms. They are developed to bridge the gap between the accuracy of
An ML potential can be generated for a given material using a database of reference structures (see section 5.1.2), a mathematical representation of the atomic structure and a regression model.^{155} The representations are based on the local environment of the active site, centered on an atom and encoding information about its neighbors, ranging from simple two or threebody terms all the way to complex “manybody” formalisms. The most commonly used representation is the “smooth overlap of atomic positions” (SOAP) as a manybody descriptor.^{156} The regression approaches, on the other hand, can be either artificial neural networks (NNs), kernelbased methods or linear regression.^{154,155}
Recently MLbased interatomic potentials have been actively generated for elements in different bonding states^{157–159} including long range interactions.^{160} Studies have also been focusing on developing a training workflow for rapidly generating machinelearned force fields (MLFFs) for investigating reaction mechanisms over catalyst surfaces. These workflows are automated and iterative where the training set eventually expands using different adaptive sampling and/or query strategies.^{154,161} The accuracy of reaction energetics predicted
Scaling relations are theoretical constructs derived from the dband model^{162} that relates the binding energies of a wide variety of catalytic intermediates to atomic adsorption energies for a range of catalyst surfaces.^{163} There are two types of scaling relations in heterogeneous catalysis: thermodynamic scaling relations^{162,164} which describe correlations between adsorption properties of chemically related species (
The development of scaling relationships has provided a rigorous theoretical basis for enabling automate reaction mechanism generation (see section 3.1), identifying independent catalyst descriptors (see section 4.3.1) and further developing volcano plots (see section 4.3.3) for catalyst screening based on DFT calculations. These relationships are extensively used for heterogeneous reactions on monometallics,^{106} bimetallics,^{166,167} oxides,^{168} and metal/support interfaces^{169} and have even been extended for application in nanocatalyst developments.^{135} Scaling relations only require a handful of data point evaluation (approx. 20–30 data points) and can thus be developed on the fly for screening catalysts or even assessing the formation of nanocatalysts and SAAs.
It should however be noted that scaling relations with prediction errors beyond 0.2 eV are not considered appropriate for facilitating predictions, especially when using scaling relations consecutively,
Volcano plots arise from correlations in reaction energetics, known as scaling relationships. They were first reported for the NH_{3} synthesis reaction^{170} and have since guided many material searches. However, there could be deviations from the known scaling relations and volcano plots under the scenarios listed in ref.
Adsorbing intermediates prefer to coordinate to a metal surface. Sometimes the need for higher coordination of adsorbate could lead to spillover to host metal,
There are also other factors reported in literature^{112} that could cause a deviation from volcano plots and scaling relation predictions. For example, intermediates adsorbed over high entropy alloys experience sitespecific Pauli repulsion interactions that can lead to deviations from scaling relations.^{171} We refer readers to ref.
The reviewed ML models and scaling relations are developed using data points generated from either DFT or DFTbased semiempirical/empirical methods. They have been found to demonstrate reasonable accuracy and can be used for highthroughput screening of catalyst materials based on reaction energetics. They can also be developed for mapping material compositions that demonstrate the required activity, stability and selectivity for a given reaction system. However the corresponding property descriptor needs to be identified, for example using convex hulls to assess thermodynamic stability of catalyst materials.^{172} Their implementation can be both exploratory and confirmatory in nature, depending on the descriptor's applicability.
This review discusses different stateoftheart tools that allow to bypass the current challenges in catalyst design, screening and optimisation. Starting with the computational bottleneck,
Once the reaction network is known, there are different approaches for theoretical catalyst optimization
Identified descriptors can be used for quick screening of relevant catalysts materials. However, highthroughput computational screening of catalysts can be elusive because descriptor evaluation is expensive, especially when the descriptor corresponds to activation barriers rather than binding energies and the material search space is large. However, inexpensive screening tools based on machine learning, secondary descriptors and onthefly scaling relations can be developed to solve this problem. ML models in particular have a range of implementations in heterogeneous catalysis,
Although, it should be noted that these databases only represent 18.9% of the total possible compositions and the perfect universal catalyst database does not exist. This is where ML techniques like active learning and transfer learning become relevant. These techniques are implemented for minimizing the number of expensive data evaluations, whereas the later focuses on bridging the accuracy of
Another challenge associated with catalyst design is the functional mapping between the required rate/selectivity and the material composition in 3D. Scaling relations and generative ML have been found to be relevant tools for screening materials and asses their synthesizability based on stability factors like convex hulls/machinelearned accelerated molecular dynamics. A potential future catalyst design approach, therefore, combines multiple ML models and provides a list of candidate 3D structures with their stability/rates/selectivity/synthesizability scores.
We are still some way away from this scenario as individual segments of such a workflow are in place,
There are no conflicts to declare.
Shambhawi acknowledges the research scholarship funding from Science and Engineering Research Board, India and Cambridge Trust. This work was in part supported by National Research Foundation (NRF), Prime Minister's Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) program as a part of the Cambridge Centre for Advanced Research and Education in Singapore Ltd (CARES), and Engineering and Physical Sciences Research Council