## Abstract

Rational catalyst design is one of the most fundamental goals in heterogeneous catalysis. Herein, we briefly review our previous design work, and then introduce a general optimization framework, which converts catalyst design into an optimization problem. Furthermore, an example is given using the gradient ascent method to show how this framework can be used for rational catalyst design. This framework may be applied to other design schemes.

## 1. Introduction

The computational design of functional materials is of great importance in material science. In particular, rational catalyst design using first-principles calculations is both fundamental and challenging, which requires not only a deep understanding of reaction mechanisms in heterogeneous catalysis, but also rational design schemes. In the last few decades, using first-principles calculations, some general relations and theories have been proposed towards an understanding of activity in heterogeneous catalysis. Hammer & Nørskov [1] found that the d-band centre of catalysts is a good indicator of the bond strength, which is useful to understand the bonding of adsorbates on transition metals. Furthermore, the Brønsted–Evans–Polanyi (BEP) relation, a linear relationship between the activation energy and the enthalpy change of an elementary step, has been found by several groups from density functional theory (DFT) calculations [2–5]. Combining the BEP relation and micro-kinetic modelling [6,7], the volcano curve/surface [4,8–11] was obtained between the catalytic activity and the bond strength, suggesting that, for a good catalyst, the bonding between the catalyst and the adsorbate should be neither too strong, nor too weak. Moreover, trends between the selectivity and the bond strength were also proposed [12,13], which offers a quantitative approach to selectivity evaluation.

As for the rational design scheme, only a few approaches have been proposed. Cheng & Hu [14] introduced the concept of chemical potential into heterogeneous catalysis, using which they found that the chemical potential of intermediates should stay in the window between the chemical potentials of the reactants and products for a good catalyst [15,16]. Examples of this design scheme will be reviewed below. Nørskov and co-workers [17,18] suggested a descriptor-based screening method. By screening a large material database, several catalysts were found for hydrogen evolution [19] and oxygen reduction reactions [20]. These approaches mentioned above are all based on screening a large database using thermodynamic properties as descriptors, and many successful examples have shown that these approaches are powerful and robust in screening catalysts. However, some disadvantages of these approaches may limit their utilization. Firstly, the performance of these approaches is closely related to the size and diversity of the screening database. Therefore, materials that are not yet discovered and stored in the database cannot be considered in the screening process. Secondly, some approximations, such as the BEP relation [21] and linear scaling [22], were involved, in which the parameters may vary for different surface structures, e.g. flat, stepped and kink surfaces. Thus, there may be some limitations in these approaches to be applied to any types of catalysts. Last but not least, these approaches are carried out using high-throughput screening, which to some extent lacks rationality. Therefore, a new framework is highly desirable for rational catalyst design in heterogeneous catalysis. In this contribution, we will briefly review our previous design work. After that, an optimization framework will be suggested along with a design example using the gradient ascent method.

## 2. Previous design work

In this section, two examples of our previous design work will be reviewed briefly, including the design of active counter-electrode (CE) materials for the triiodide reduction reaction and bi-phase catalysts for hydrogenation of tertiary amides. More details of these design examples were recently reviewed and can be found elsewhere [23].

Owing to their low cost, ease of fabrication and relatively high energy conversion efficiency, dye-sensitized solar cells (DSCs) [24–26] have attracted great attention as a photovoltaic device for solar energy conversion. However, the large-scale application of DSCs is greatly restricted by the current CE material for the triiodide reduction reaction, platinum, because of its high cost and low abundance. In order to search for potential non-Pt CE materials with high electrocatalytic activities, a general and efficient screening framework was suggested in our previous work [27]. There are three elementary steps in the triiodide reduction reaction (IRR), namely:
2.1
2.2
2.3where * stands for the free site on the electrode surface and sol indicates acetonitrile solution. Based on the chemical potential [14] and the optimal adsorption energy window [15,16] mentioned in the Introduction, the Gibbs free energies of the states in these elementary steps should decrease step by step to make the overall IRR proceed in the forward direction, from which the chemisorption energy of an I atom on good catalysts should be approximately in the range from 0.33 to 1.20 eV (derivation details can be found in [27]). We calculated the adsorption energies of an I atom on many materials and screened some promising CE materials according to the adsorption energy window mentioned above. Several known good catalysts for IRR were found within the adsorption energy range, such as Pt, WO_{3}, WC, NiS, CoS, MoC and MoN. Furthermore, two unreported surfaces were also found in the estimated range of suitable CE material for DSCs, namely the Fe_{2} O_{3}(012) and Fe_{2}O_{3}(104) surfaces of *α*-Fe_{2}O_{3}. To verify the theoretical predictions, pure-phase Fe_{2}O_{3} was synthesized via a simple modified hydrothermal method [27]. From photovoltaic data, the overall energy conversion efficiency of the *α*-Fe_{2}O_{3} was found to be 6.96% (more details can be found in [27]), which is remarkably comparable to that of DSCs with Pt CE (7.32%). These results suggest that *α*-Fe_{2}O_{3} is sufficient to catalyse the reduction of triiodide to iodide.

Amines are important organic compounds, which are widely used in many chemical industries, and the reduction of amides is a potentially simple method to produce amines. However, the low catalytic activity of amide reduction is the main obstacle for this process. In our previous work, a bi-phase catalyst for hydrogenation of tertiary amides was designed [28]. In the original work, the free energies of intermediates and transition states in catalytic hydrogenation of *N*-methylpyrrolidin-2-one to *N*-methylpyrrolidine and water were calculated using DFT and thermodynamic analysis [29,30] on four stepped transition metal surfaces, namely Re, Ru, Rh and Pt. We found that these surfaces can be divided into two groups: on Re surfaces, the hydrogenation of oxygen to water is the rate-determining step, while the barriers of the C=O bond activation are very high on Ru, Rh and Pt. Therefore, the mono-phase is limited by the balance between having sufficient oxophilicity to break the C=O bond and enabling facile oxygen hydrogenation [28]. Adding more phases may be possible to break this constraint and offer bi-functional sites to enhance both the bond cleavage of C–O and the hydrogenation of oxygen. Therefore, we systematically investigated the hydrogenation of tertiary amides on several bi-phase catalysts, namely Re/Pd, Re/Rh, Re/Ru and Re/Pt. The highest activity of amide hydrogenation was found on Re/Pt bi-phase surfaces. Furthermore, we prepared and tested a range of bimetallic Pt/Re catalysts, and over 90% conversion and more than 99% selectivity to the corresponding amine were found experimentally [28], suggesting that multi-phase catalysts are promising in the reduction of amides [28] and some other reactions [31].

## 3. Optimization framework

Designing the most active heterogeneous catalyst for a given reaction involves searching for one catalyst structure that possesses the highest overall turnover frequency (TOF) of this reaction than any other catalysts. Therefore, this process can be treated as an optimization problem as follows:

*Given*: a function *f*:*S*→TOF from a set of catalyst structures, *S*, to the corresponding set of TOF values, TOF

*Sought*: a catalyst structure **r**_{op} in *S*, such that *f*(**r**_{0})≥*f*(**r**) for all **r** in *S*.

With this definition, the catalyst design is transformed into a concrete optimization problem.

In this optimization process, function *f* calculates the value of TOF from the catalyst structure. However, the exact expression of this function is still unknown. With the help of DFT calculations and micro-kinetic modelling, for a given reaction, the TOF on a specific catalyst can be obtained routinely as shown in the real calculation part in figure 1. For example, in the dry reforming of methane on nickel-based catalysts [32], the catalyst structure is Ni(111) with four different reaction sites: top, bridge, fcc (face-centred cubic) and hcp (hexagonal close-packed) hollow sites (figure 1*a*). Using the DFT approach, the energies of all intermediates and transition states involved in the dry reforming were calculated, and, with some thermal corrections, the free energy profile can be plotted as shown in figure 1*c*. With the free energy changes and barriers of all the reactions, the TOF of dry reforming can be calculated to be 23.49 s^{−1} on Ni(111) using micro-kinetic modelling [32] (figure 1*e*). Therefore, function *f* can be divided into two functions as shown in the mathematical model part in figure 1:

—

*Function I*. A function*g*:*S*→*E*from a set of catalyst structures,*S*, to the corresponding set of energies of all the intermediates and transition states,*E*, i.e.**e**=*g*(**r**) in which**r**∈*S*and**e**∈*E*.—

*Function II*. A function*h*:*E*→ TOF from a set of energies of all the intermediates and transition states,*E*, to the corresponding set of TOF values, TOF, i.e. rate =*h*(**e**) in which rate ∈ TOF.

Thus, function *f* can be written as
3.1Using function *g*, the energies of all the intermediates and transition states involved in the reactions can be calculated, which can be obtained using DFT calculations, while function *h* combines all the thermodynamic and kinetic parameters together to calculate the TOF values using micro-kinetic modelling.

Concerning the two variables involved (**r** and **e**), they can be represented as matrices or vectors: as shown in figure 1, variable **r** for Ni(111) is written in the form of an *N*×4 matrix (*N* is the number of atoms in the catalyst), where the symbols and the positions of all the atoms are listed. For the atom positions, they can be written either in Cartesian coordinates or in polar coordinates. On the other hand, the energy variable, **e**, can be written as a matrix or vector containing the energies of all intermediates and transition states. With these notations, the calculation from a set of structures to the corresponding set of overall activities is converted to a mathematical model from matrix **r** to TOF values as shown in figure 1. Therefore, the process of designing the most active catalyst is converted to an optimization of searching for a matrix **r**_{op} so that *h*(*g*(**r**_{0})) is the highest than any other possible matrix **r**.

## 4. Optimization example

After introducing the optimization scheme for catalyst design, here we give an optimization example to show how this optimization framework can be used with one of the most simple and frequently used optimization methods, namely gradient ascent. A model reaction was considered as shown in figure 2, which contains two intermediates (IS1 and IS2) and three transition states (TS1, TS2 and TS3). Thus, the energies of seven states, i.e. the above five states plus the states of reactant and product (*E*(R) and *E*(P)), are involved in this system, among which *E*(R) and *E*(P) are constant and independent of catalyst surfaces. Therefore, there are five elements in variable **e**, namely *E*(TS1), *E*(IS1), *E*(TS2), *E*(IS2) and *E*(TS3). To optimize *f*(**r**) using the gradient ascent method, an initial guess of **r** is necessary, which can be in practice set to a known catalyst. We assume that the structure of the known catalyst is **r**_{0} containing *N* atoms in the unit cell, and the energy profile of this catalyst is shown in figure 2. According to the gradient ascent method, the TOF function *f*(**r**) increases fastest if variable **r** goes in the direction of the positive gradient of function *f* at point **r**_{n},∇*f*(*r*_{n}). Therefore, the surface structure can be updated in the sequence of **r**_{0},**r**_{1},**r**_{2},…, i.e.
4.1where *γ* is the step size allowing the change of surface structure **r**_{n} in the optimization, and the gradient of function *f* at point **r**_{n} is ∇*f*(**r**_{n}).

In equation (4.1), calculating the gradient of function *f* at point **r**_{n},∇*f*(**r**_{n}), is not easy due to the complexity of function *g* from surface structure **r** to the energies of all the states, **e**. We first consider the gradient of function *h*(**e**_{n}), in which **e**_{n}=*g*(**r**_{n}). In our model reaction shown in figure 2, **e** contains five elements, *E*(TS1), *E*(IS1), *E*(TS2), *E*(IS2) and *E*(TS3), as mentioned before. Therefore, the gradient of function *h*(**e**),∇*h*(**e**), can be written as follows:
4.2In order to increase the value of TOF, the five elements of energy vector **e** should move in the positive direction of ∇*h*(**e**). More specifically, for the catalyst with the energy profile in figure 2, the second step is the rate-determining step with the highest effective barriers along the reaction coordinate (more detailed definitions of effective barriers can be found in [10,11,32]). In this step, the energy of TS2 is too high, leading to the reaction being kinetically unfavourable, while the energy of IS1 is too low, which may poison the catalyst, indicating that the values of ∂rate/∂*E*_{TS2} and ∂rate/∂*E*_{IS1} should be negative and positive, respectively. Thus, in order to optimize the catalytic activity, the energies of TS2 and IS2 should move in the direction of ∂rate/∂*E*_{TS2} and ∂rate/∂*E*_{IS1}, respectively, which means that TS2 should be more stable, while IS2 should become less stable. Therefore, the optimization of TOF is transformed to optimizing the energy vector, **e**. It is interesting that the gradient expression of function *h*,∇*h*(**e**), is similar to the one in the general degree of rate control proposed by Campbell and co-workers [33].

To optimize the overall reaction rate, we can let the energy vector, **e**, move in the direction of ∇*h*(**e**). However, in some complex reactions, such as Fischer–Tropsch synthesis and methane dry reforming, several tens or even hundreds of energies are involved. Considering the current computation capacity, some approximations may be needed to simplify this optimization. Regarding the elements in the energy vector, **e**, they are in fact not independent of each other. Firstly, the energies of transition states are correlated with the enthalpy changes of the corresponding steps, known as the BEP relation [21], which can be introduced into the optimization to approximately halve the size of the energy vector. Secondly, Nørskov and co-workers [22] proposed a linear scaling, which states that the adsorption energies of many species scale approximately with the adsorption energies of the central atoms. Furthermore, Cheng & Hu also found that the adsorption energies of all the other intermediates can be represented by the adsorption energies of the key intermediate [10–12]. Therefore, the energies of the intermediates can be further simplified into one energy, the energy of the key intermediate. On the other hand, from the kinetic perspective, in most cases, the overall activity is determined by only one or two steps on a catalyst. Therefore, another approximation may be introduced, in which only the energies related to the rate-determining steps, i.e. the prominent elements in ∇*h*(**e**), are considered. With these approximations, we can write the gradient of function *h*,∇*h*(**e**), with only one energy term. For example, in the reaction shown in figure 2, the gradient of function *h*,∇*h*(**e**), can be written in terms of the partial derivative of the energy of IS1, which is both the key intermediate in this energy profile and the initial state of the rate-determining step.

Using the chain rule, the gradient of function *f*,∇*f*(**r**), with respect to catalyst structure **r** for the reaction in figure 2 can be finally expressed as follows:
4.3With the approximations mentioned above, we only need to consider the term of IS1, (∂rate/∂*E*_{IS1})(∂*E*_{IS1}/∂**r**). In this term, ∂rate/∂*E*_{IS1} is positive as shown above, suggesting that increasing the energy of IS1 can speed up the overall reaction. Therefore, ∂*E*_{IS1}/∂**r** can offer the direction, which points out that the structure should be changed to increase the energy of IS1 in order to enhance the overall reaction rate. However, calculating this derivative is very difficult for the following reasons. Firstly, the structure matrix includes *N*×4 (*N* is the number of atoms in the catalyst) elements, which are very time-consuming to calculate. Secondly, the structure matrix contains the symbol of the atoms, which is not continuous and differentiable. Thirdly, except for DFT calculations, there is hardly any general method to rapidly predict the chemisorption energies of adsorbates with reasonable accuracy.

Herein, we propose the following solution for the problems mentioned above. Concerning the matrix size issue, we can write the structure matrix using polar coordinates, in which the adsorption site of IS1 is set as the pole. Therefore, we can distinguish the symmetric atoms easily using the distance term in the polar coordinates, and calculate all the symmetric atoms. Moreover, the longer the distance between the catalyst atom and the adsorbate, the lower the contribution to the adsorption energy. Thus, we can set a cut-off for the distance and only take the atoms in the nearby atom shells into account, depending on the computation capacity. For the non-continuous problem of the atom symbol, we need to find some properties to replace the atom symbol term. The adsorption energy of the adsorbate on the close-packed slab of the atoms may be suitable to describe the activity of the atoms. For instance, in the reaction shown in figure 2, the atom symbol term of nickel can be represented by the adsorption energy of IS1 on Ni(111). Considering the adsorption energy is the overall contribution of all the nearby atoms, and, in most situations, changing one of them only slightly modifies the adsorption energy. Therefore, substituting one atom of the catalyst surface may be treated as approximately continuous. With the gradient of function *f*,∇*f*(**r**), the surface structure can be optimized in the direction of ∇*f*(**r**) to raise the energy of IS1 and thus increase the overall activity according to equation (4.1).

## 5. Design example

In this section, we propose a catalyst design example for CO oxidation to show how our optimization framework can be used to design better catalysts in real catalytic systems. Catalytic CO oxidation is chosen for the following reasons. On the one hand, it is of great importance in many technologies such as car-exhaust emission control and CO_{2} lasers [34], and thus designing catalysts for CO oxidation may greatly enhance these technologies. On the other hand, CO oxidation is a simple and typical reaction in heterogeneous catalysis, including adsorption, dissociation, association and desorption. Therefore, the design strategy shown in this section may be generally applied in some other reactions in heterogeneous catalysis. As suggested above, in order to optimize the activity of one catalytic system, an initial guess is needed. Nørskov and co-workers [35] reported that under high-temperature conditions (*T*=600 K, *P*_{O2}=0.33 bar and *P*_{CO}=0.67 bar) Pd is one of the most active catalysts for CO oxidation [35]. Therefore, the palladium surface used in their work, namely *p*(2×2)-Pd(111), was chosen as the initial guess in the optimization as shown in figure 3*a*. With the initial guess chosen, we then need a strategy of manipulating the surface structure to tune the energies and thus optimize the TOF. There are several possible manipulating strategies considering the structure matrix in figure 1. (i) Only the symbols in the matrix are allowed to be changed, which is widely used in heterogeneous catalysis known as alloying or doping. (ii) Only the positions in the matrix are allowed to be modified, which may be used to construct some different facets or cluster structures. One needs to note that in this strategy the surface structure is also constrained by stability. (iii) The elements in the matrix are allowed to be removed, or new elements are allowed to be added. This type of manipulation can be used to tune the adsorption energies by introducing a defect or changing the coordination number of the adsorption site. (iv) These three methods can also be combined together to build complex catalyst structures. In this work, the first method, namely the doping method, was chosen because it is very easy to implement experimentally, and several outstanding alloy catalysts were found to be very active in many reactions [13,19,36,37]. On the initial and optimized surface, we calculated the energies of intermediates and transition states involved in CO oxidation, and the activities of these surface were evaluated using micro-kinetic modelling. All the DFT calculations were carried out using a periodic slab model with the Vienna *ab initio* simulation program (VASP) [38–41], while the thermodynamic correction and micro-kinetic modelling were carried out using the CatMAP [42–44] code (more computational details can be found in the electronic supplementary material).

We first calculated the adsorption free energy of CO, and the free energy barriers of O_{2} dissociation and CO oxidation on the initial Pd surface (function *g*), and based on these energies the TOF was calculated as 2.49 s^{−1} using micro-kinetic modelling (function *f*). All these results are listed in table 1. As suggested above, we need to calculate the gradient of function *f* to determine the direction of manipulation. We first calculated the gradient of function *h*(**e**),∇*h*(**e**), as
5.1Because the energy is exponentially related to the TOF, in this real example we rewrote the partial derivative of species *i* (PD_{i}) in terms of free energy as
5.2where *R* and *T* stand for the gas constant and reaction temperature, respectively, *G*_{i} is the free energy of one state, while *G*_{j} are the energies of all the other states. From our micro-kinetic results, ∇*h*(**e**) was calculated to be
5.3Therefore, the surface should be modified to tune the energies in the direction ∇*h*(**e**). To simplify this problem, we only took the most significant partial derivative into account, namely PD_{CO}. After that, considering the overall function, *f*, the gradient of function *f*, ∇*f*(**r**), is represented as
5.4From the above equation, the last unknown term in the optimization is the partial derivative of CO energy, ∂*E*_{CO}/∂**r**. In the four-layer *p*(2×2)-Pd(111) surface, there are four different types of atoms, including the surface atom, subsurface atom, third-layer atom and fourth-layer atom, and any one of them may be substituted in the manipulation. In this work, we only considered the substitution of the surface atom, and the other types of substitution may be carried out in a similar approach. For a surface substitution, there are two directions, namely substitution with a more active metal or a less active metal. From the periodic table, a more active metal is Rh, which is on the left of Pd, while a more inert metal is Ag. However, from the segregation energies calculated by Nørskov and co-workers [45], only Rh is favourable to form an alloy with host Pd, while Ag is very likely to segregate with Pd, which may not be suitable in this optimization. Copper is another choice of the inert metal for gradient calculation, which is favourable to form an alloy in host Pd with a segregation energy of −0.20 eV/atom [45]. Therefore, we chose Rh and Cu as the solute metals to compute the partial derivative of CO energy. As suggested in §4, to solve the non-continuous atom symbol, the adsorption energy of CO on the close-packed surface of one metal is used to represent the symbol of this metal. Therefore, the Rh and Cu metals are quantified as −0.05 and 1.06 eV, respectively, which are the adsorption energies of CO on Rh(111) and Cu(111), with respect to the adsorption energy of CO on Pd(111) (negative value suggests stronger adsorption than that on Pd(111)). After the substitution, the modified surface structures are shown in figure 3*b* for Cu and figure 3*c* for Rh. On these surfaces, the adsorption energies of CO were calculated to be 0.12 eV on the Cu-doped surface and −0.18 eV on the Rh-doped surface. Therefore, based on the results above, the partial derivative of CO energy, ∂*E*_{CO}/∂**r**, is calculated as follows:
5.5Thus, the gradient of function *f*,∇*f*(**r**), is equal to 0.53 using equation (5.4), suggesting that the surface atom should be substituted with a more inert solute metal. According to equation (4.1), surface structure **r** should be updated using ∇*f*(**r**), depending on the step size, *γ* . To qualitatively check the feasibility of our framework, in this work we used the Cu-doped surface as the optimized surface, i.e. the step size is set to 2. A more systematic and quantitative optimization example will be presented in our future work.

To compare the activities of the initial and optimized surfaces, we calculated the adsorption free energy of CO, the activation free energies of O_{2} dissociation and CO oxidation on the optimized Cu-doped surface (table 1), and the corresponding geometries are shown in figure 4. From these results, it can be seen that the adsorption free energy of CO decreases from −0.64 eV on the initial Pd(111) surface to −0.46 eV on the optimized surface. Furthermore, both O_{2} dissociation and CO oxidation became more favourable kinetically. Using micro-kinetic modelling, the activity of CO oxidation on the optimized surface (figure 3*b*) was calculated to be 21 348.94 s^{−1}, which is several orders of magnitude higher than that of the initial surface. This simple example shows that our optimization framework can be used robustly and rationally in catalyst design. Even though the increase in activity is purely from DFT calculation and micro-kinetic modelling, it may offer some suggestions or catalyst candidate for experiment.

It may be worth noting that the optimization framework proposed in this work for rational catalyst design may cover some important features of most current schemes for catalyst screening or design. For example, the screening method can be treated as an optimization process, which uses the exhaustive search method as the optimization strategy, while the optimization example using the gradient ascent method suggested in the current work contains some rational features. Furthermore, the catalyst design is given as a typical optimization problem, and therefore other optimization methods can also be considered, including some of the state-of-the-art global optimization algorithms, such as evolutionary algorithms and swarm-based optimization algorithms. Moreover, even though some approximations are included in this work, the framework should also work without any approximations, which can offer a design scheme with the accuracy of DFT. Other approximations can be introduced into this framework, depending on the design system. We can also expect that some empirical or semi-empirical methods to rapidly predict the energy of adsorbates from a surface structure are highly valuable.

## Authors' contributions

P.H. designed the project. All the calculations were carried out by Z.W. The paper was drafted by Z.W. and revised by P.H.

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Acknowledgements

The authors gratefully acknowledge the UK’s national high-performance computing service ARCHER (for which access was obtained via the UKCP consortium) for computing time. P.H. thanks the National Natural Science Foundation of China (21333003) and Z.W. thanks Queen’s University Belfast for a PhD studentship.

## Footnotes

One contribution of 15 to a discussion meeting issue ‘Catalysis making the world a better place’.

- Accepted November 10, 2015.

- © 2016 The Author(s)