An approach to calculate the free energy changes of surface reactions using free energy decomposition on ab initio brute-force molecular dynamics trajectories

Abstract

To understand the mechanisms and kinetics of catalytic reactions in heterogeneous catalysis, ab initio molecular dynamics is one of the powerful methods used to explore the free energy surface (FES) of surface elementary steps. The most significant aspect of performing such calculations is to choose the specific collective variable (CV) of the reaction. Here, we take CO oxidation on Pt(111) at 300 K as an example to demonstrate the protocol of selecting CVs guided by the free energy decomposition which quantifies individual bond free energy contributions. The basic concept is to conduct the brute-force molecular dynamics initiated from the transition state on the FES, which is refined from the one on the potential energy surface, to generate the reaction path at a finite temperature. The validity of this reaction path is further demonstrated by a 2-D free energy landscape spanned by the path-CV. By choosing CVs including other bond distances, we find that CO oxidation cannot be well understood by umbrella sampling or constrained molecular dynamics (CMD) solely along the OC–O bond distance. The free energy decomposition analysis suggests that not only the OC–O bond but also two O–Pt bonds are responsible for the free energy change. The further CMD simulations along selected CVs based on the insights from our protocol capture different reaction stages and give solid estimations of free energy barriers.

---

1. Introduction

Understanding reaction mechanisms at the atomic level has always been one of the main tasks of computational heterogeneous catalysis. In the last few decades, the entire field has devoted itself to exploring a large number of reactions by searching for local minima and saddle points on the potential energy surface (PES). The transition-state theory bridges the gap between equilibrium properties and kinetics with the applicability of free energies (Gibbs and Helmholtz). For surface reaction calculations,^1–3 a specific free energy is usually estimated as the potential energy plus the statistical mechanics correction.⁴ However, the most inviting heterogeneous catalytic reactions are carried out under situations far from the ideal conditions that we adopt in calculations. Some reactions such as oxidation and reduction can only happen at a relatively high temperature which leads to surface disorder from the 0 K lattice structure.^5,6 Other reactions can take place at ambient temperature but have to occur at solute–solid interfaces especially electrocatalysis^7,8 and photocatalysis.^9,10 Moreover, the structures of nanoclusters change rapidly during reactions, which cannot be treated as frozen solids.^11,12 In the studies related to the above reactions, ab initio molecular dynamics (AIMD) based on the density functional theory (DFT) is usually adopted, which offers novel understanding with acceptable accuracy. Not only new insights into reaction mechanisms but also free energy values with much better accuracy can be obtained from AIMD.^13,14

Since the difference between the Gibbs free energy and the Helmholtz free energy is very small for surface reactions, the latter is usually evaluated using AIMD. To explore the free energy surface (FES) using molecular dynamics simulations, plenty of strategies have been developed, namely constrained molecular dynamics (CMD),¹⁵ umbrella sampling (US),^16,17 metadynamics^18,19 and others.^20–22 In the framework of the above methods, the free energy change is associated with a particular coordinate known as the collective variable (CV). The free energy at a given value of the CV ξ₀ can be determined by

where P(ξ₀) is the probability density of finding the system with the CV ξ equal to ξ₀. The free energy change along the progress of the CV is the free energy profile which is often referred to as the potential of mean force (PMF). The CV selection for some reactions is quite obvious. Normally, they are structural parameters (bond distances, angles or torsions) changing drastically throughout the reaction. The linear combination of a few forming and breaking bonds can also be used for reactions involving many atoms.²³ In addition, the path-CV^24,25 devised by Branduardi et al. is suitable for studying the free energy change along a preassigned path by offering a set of reference structures.

Nowadays, the CV selection deeply relies on how we understand the reaction of interest based on chemical insights. Various methods have been proposed to efficiently extract suitable CVs for a given reaction. Most of them employ mathematical analysis of metastable state fluctuations²⁶ or MD trajectories.²⁷ When it comes to surface reactions, many bonds can be affected due to the interaction between the adsorbates and the surface. As a consequence, choosing bonds responsible for the free energy change is a non-trivial task. Haas et al. tried to understand biochemical reactions by calculating individual atom free energy contributions along the reaction path.²⁸ As bond distances frequently serve as CVs, the bond free energy contribution is a better descriptor than the atomic one. Therefore, we suggest a protocol which analyses individual bond free energy contributions by utilizing the information of free energy gradients (FEGs) along the reaction path and further guides the CV selection in the current work. The elementary step of CO oxidation on the Pt(111) surface is chosen and tested for the following reasons. This reaction is perhaps the most typical model surface reaction. As one of the noxious pollutants resulting from power generation and automobiles, it is necessary to remove CO. Catalytic CO oxidation as an a posteriori process answers the demands for environmentally friendly technology in today's society.^29–31 Among the steps in the overall process, the elementary step from the co-adsorbed CO and O to the adsorbed CO₂ is arguably the most important one, the chemical equation for which is

CO* + O* → CO₂* (2)

where * denotes the adsorption state. Pioneering DFT exploration of its potential energy profile²⁹ and later theoretical work^32,33 unravelled CO oxidation on the Pt(111) surface at the atomic scale. Recently, a few AIMD studies focusing on the post-transition-state dynamics investigated the process of the CO₂ desorption.^34,35 In the current work, we study the above elementary reaction to demonstrate the above protocol, aiming to gain more insights into this process from the perspective of the FES.

2. Computational methods

2.1 Electronic structure calculation and molecular dynamics simulations

All electronic structure calculations were performed using the Vienna Ab initio Simulation Package (VASP)³⁶ in the framework of DFT at the level of generalized gradient approximation (GGA) applying the Perdew–Burke–Ernzerhof (PBE)³⁷ functional. The projector-augmented-wave (PAW)³⁸ pseudopotential with a kinetic energy cut-off of 400 eV was adopted to describe core-electrons, which is the largest among the three elements in the system (400 eV for O). The first order Methfessel–Paxton³⁹ scheme was used with a smearing width of 0.2 eV. The Brillouin zone was sampled using a 12 × 12 × 12 Gamma-centred k-points scheme⁴⁰ for Pt bulk and 4 × 4 × 1 for the Pt(111) surface. The optimized lattice length of the FCC Pt bulk is 3.97 Å which is close to the experimental result of 3.91 Å.⁴¹ The Pt(111) surface model was constructed as a four-layer p(2 × 2) unit cell with a ∼12 Å vacuum slab along the Z direction. The bottom two layers were frozen during all ionic position updates. For geometric optimizations, all the forces on flexible atoms were minimized to less than 0.05 eV Å⁻¹.

The Born–Oppenheimer MD (BOMD) was used for AIMD simulations. Adsorbates and the top two layers were updated using the Verlet algorithm with a timestep of 1.0 fs, coupled with the Nosé–Hoover⁴² thermostat to maintain the target temperature of 300 K. Thus, the free energy mentioned below is the Helmholtz free energy corresponding to the NVT ensemble. The standard brute-force MD and enhanced sampling methods including the CMD and the US were realized using the code Bucko implemented in the VASP.⁴³ For the CMD, each structure was first equilibrated for 0.5 ps and later 3 ps to average the FEG. The block average algorithm⁴⁴ was used to assess the statistical uncertainty of the FEG, which is implemented by Spencer.⁴⁵ The integration of FEGs was performed simply using the trapezoidal method. In the US, the Gaussian function was taken as the biasing potential with a height of −5.0 eV and a width of 0.5 Å. The samples were collected during a 3 ps simulation after a 0.5 ps equilibration. The code of the WHAM⁴⁶ analysis used to reconstruct the PMF was developed by Grossfield.⁴⁷ For free energy calculations along the path-CV, our homemade code was utilized to carry out the US along the path-CV according to velocities and forces given by the VASP. Two components of the bond referred path-CV s and z took the form of⁴⁸

where Ξ is path vector whose entries are primitive CVs, P is the number of reference structures and λ is the coefficient which smooths the path. More specific parameters for the US and the CMD and mathematical derivations of the path-CV can be found in the S2 and S5 of the ESI† separately.

2.2 Protocol of the free energy decomposition

In order to understand the reaction from the free energy perspective, we put forward a protocol which decomposes the total free energy change into the free energy contribution of individual bonds. This method can be readily transferred to other surface reactions because it utilizes the advantage of the CMD and begins with the transition state (TS) on the PES. More details will be provided by applying it to the reaction of the CO oxidation on the Pt(111) surface below. The procedure is organized as follows:

(i) Find the TS on the PES with frequently used search methods such as the constrained optimization^49,50 or the nudged elastic band (NEB) method.⁵¹ Refine the above TS to the one at the target temperature according to the FEG given by the CMD.⁵² The FEG on this CV should be minimized to a threshold near zero.

(ii) Generate the reaction path by using the brute-force MD starting from the TS forwards and backwards to the final state (FS) and the initial state (IS), respectively. The outcome path intrinsically comprises the reaction environment change due to the finite temperature including the surface relaxation.

(iii) Carry out the CMD on selected equidistance structures along the reaction path to obtain individual bond free energy changes, which guides the selection of CVs for various stages of the reaction. Later, enhanced sampling methods can be performed to acquire better interpretable PMFs along selected CVs.

3. Results and discussion

3.1 Potential energy profile of the CO oxidation on the Pt(111) surface

We started our investigation by calculating the potential energy profile of the CO oxidation on Pt(111). To better understand the reaction, here we recap the structure characteristics of the Pt(111) surface. Fig. 1 denotes a few adsorption states and key bond distances used in our work. There are four adsorption sites on the surface, namely the hexagonal close packed (HCP) site, the face centred cubic (FCC) site, the bridge site and the top site. The IS with the lowest potential energy is that CO is on the top site while O co-adsorbs on the FCC site in the p(2 × 2) unit cell.

Fig. 1 Notations used through the article. (a) Four Pt atoms in the first layer are marked as Pt_A to Pt_D. The FCC, HCP, bridge and top sites are where CO and O can be adsorbed on. (b) The structure of the initial state with the lowest potential energy is illustrated. Bond distances R_OC–O, R_O–PtA, R_O–PtB, R_O–PtC and R_C–PtD are important structural features during the reaction.

This reaction is driven by the interaction between the co-adsorbed CO and O and Fig. 2 shows four key structures during the reaction. Beginning from the co-adsorbed IS, O moves to the bridge site which provides the opportunity for CO to interact with O. The TS is located at the HCP site where CO weakly bonds to O. Finally, the formed CO₂ will move from HCP site to the bridge site, which then desorbs into the gas phase. Intuitively, the distance of the forming OC–O bond (R_OC–O) can serve as the CV. However, microscopically, several bonds are generated and cleaved consecutively due to the interaction between the adsorbates and the surface. It may be inappropriate to attribute all of the free energy change to the sole R_OC–O change and the interpretation of this reaction may not be completely convincing without taking other bonds into consideration. Obviously, R_OC–O changes explicitly but the other four bonds also vary during the whole reaction. R_O–PtA and R_O–PtB have to elongate until these two bonds break when O moves from the FCC site to the top site at Pt_C. If the CO₂ finally desorbs, R_O–PtC will also change. In addition, CO is restricted by the C–Pt_D bond around the top site, which varies to allow C to interact with O during the whole reaction.

Fig. 2 Key structures of the CO oxidation on the Pt(111): (a) the initial state; (b) the transition state; (c) the final state; and (d) the physisorbed CO₂. In this study, the elementary reaction from (a) to (c) is extensively investigated by the density functional theory calculation and ab initio molecular dynamics. The p(3 × 3) unit cell is used only for the illustrative purpose. All calculations were carried out using the p(2 × 2) unit cell.

Then we calculated the potential energy changes along the path, starting with the co-adsorbed CO and O and ending with the adsorbed CO₂. Fig. 3 illustrates the potential energy profile along the minimum energy path (MEP) via nudged elastic band (NEB)⁵¹ calculations. Direct thermo-corrections were applied to estimate the corresponding free energies at 300 K (the mathematical formula of which can be found in the ESI† S1). Since AIMD here does not retrieve the nuclear quantum effect, the zero point energy is excluded from the free energy correction for the further comparison. The potential energy barriers are similar to those reported in recent articles,^34,53 which assures the accuracy of the electronic structure calculation setting in this study. The reaction process has a forward potential energy barrier of 0.92 eV and 0.97 eV for the backward step while free energy barriers are 0.95 eV and 1.00 eV, respectively. The energy profile is almost symmetrical except for the turning point located in the last half reaction, where the process of the CO₂ rotation from the HCP site to the bridge site takes place.

Fig. 3 Potential energy profile of the CO oxidation on the Pt(111) surface along the minimum energy path. Seven structures were linearly interpolated between the initial state and the final state. Forces perpendicular to the path were minimized smaller than 0.05 eV Å⁻¹. The forward potential energy barrier is 0.92 eV while the backward step is 0.97 eV. Thermodynamic corrections are also applied to obtain the free energies for the IS, the TS and the FS. The resulting forward and the backward free energy barrier is 0.95 eV and 1.00 eV respectively.

3.2 Potential of mean force along the OC–O bond distance

Towards the free energy calculation of this reaction, we accomplished the US and the CMD by taking R_OC–O as the CV first. These two methods are most frequently used to calculate the PMF in the AIMD. Both methods sampled R_OC–O from 3.3 Å to 1.2 Å with an interval of 0.1 Å. Detailed FEGs with their standard errors from the CMD and free energies with errors of the US are listed in ESI† S3.

The US gives a forward free energy barrier of 0.83 eV, which is about 0.1 eV lower than the direct thermodynamically corrected results 0.95 eV. In the thermodynamic corrections the harmonic approximation is used, while the US actually samples all the dynamic motions at the finite temperature. The deviation between these two values can be attributed to the existence of anharmonicities in the system.⁵⁴ Although the free energy curve seems smooth, the sampling along R_OC–O from the TS to the FS appears to be peculiar. For the PMF from the US, the FS has a much lower free energy as shown in Fig. 4. Since the US estimates the free energy from the relative probability of sampled CV values, we cannot determine the stable states on the FES when unaware of the FEGs. The R_OC–O in the optimized adsorbed CO₂ on the PES is 1.290 Å, and the R_OC–O in the FS on the FES is expected to be about 1.3 Å. However, the histogram of samples from the US illustrates that the number of samples at 1.3 Å is significantly lower than those in other regions. As its barrier is about 0.2 eV,⁵⁵ the desorption of chemisorbed CO₂ over the bridge site can be easily observed within few hundred MD steps at 300 K depending on the initial structure and velocities. Meanwhile, the restraining potential was only applied to the reacting OC–O bond while the C–Pt_D and O–Pt_C bonds that anchor the molecule on the surface are rather flexible in our simulations, and CO₂ would form readily and fleet into the gas phase during the simulations initiated from the structure with the R_OC–O between 1.2 Å and 1.4 Å. As a consequence, the most samples with R_OC–O of about 1.3 Å are contributed by CO₂ in the gas phase, which unlikely describe the behavior of the adsorbed CO₂. The FS with the free energy of −0.42 eV is actually the gas CO₂ ensemble. Moreover, the interpretation of the reaction based on R_OC–O is ambiguous, which will be discussed below.

Fig. 4 Potential of mean force and numbers of samples collected in various bins for the CO oxidation on the Pt(111) surface along R_OC–O by the umbrella sampling. A Gaussian function with 5 eV depth and 0.5 Å width worked as the restraining potential. The blue curve is the PMF while colourful bars are numbers of samples collected during simulations. The forward free energy barrier is 0.83 eV. This method fails to determine the free energy of the final state because samples representing the adsorbed CO₂ are difficult to collect during simulations initiating from structures with R_OC–O between 1.2 Å and 1.4 Å.

On the other hand, the CMD presents a similar PMF to the one obtained from the US by integrating FEGs at sample points with a little larger forward free energy barrier of 0.86 eV and the backward free energy barrier of 0.92 eV. The FEG decreases linearly with the decrease of R_OC–O from 3.3 Å to 2.2 Å but suddenly drops to zero with only 0.2 Å reduction in R_OC–O. Moreover, intermediate states on the FES should have zero FEG which is analogous to the case of the PES.⁵⁶ The requirement of the zero FEG for intermediate states indicates that key states of the reaction can be located along R_OC–O, which is 3.3 Å, 2.0 Å and 1.28 Å for the IS, the TS and the FS, respectively. Fig. 5 illustrates the PMF with the key states marked. The abnormally large FEG at 1.2 Å is the result of the large tension of the C–O bond in the gas CO₂. From the structures sampled in the CMD simulations, we find that CO₂ could feasibly rotate from the HCP site to the bridge site when R_OC–O is smaller than 1.5 Å. Thus, the collected samples after the TS can be classified into two categories. When R_OC–O is between 1.5 Å and 2.0 Å, most samples are the OC–O complex over the HCP site. For R_OC–O being smaller than 1.5 Å, adsorbed CO₂ on the bridge site constitutes the majority of samples. Although the CMD gives a better performance in describing the FS than the US, the samples collected in simulations from the two categories imply the deficiency of R_OC–O in understanding this reaction.

Fig. 5 Potential of mean force for the CO oxidation on the Pt(111) surface along R_C–O by the constrained molecular dynamics. The SHAKE algorithm was used to constrain R_OC–O at given values. The generated forward free energy barrier is 0.86 eV. The key states, the initial state, the transition state and the final state are marked by stars, located at 3.3 Å, 2.0 Å and 1.28 Å respectively. The extremely large free energy gradient at 1.2 Å is majorly caused by the tension of the OC–O bond.

3.3 Collective variable selection by the free energy decomposition

The incapability of describing the adsorbed CO₂ state implies that R_OC–O alone may not be the most suitable CV. But it is extremely time-consuming to test whether a structural parameter does contribute to the reaction free energy change rather than fluctuations caused by its periodic thermal vibration. To demonstrate the protocol of the CV selection by the free energy decomposition along the reaction path, we first tested the first two steps in Section 2.2. The MEP of the CO oxidation displays that there are two bonds (R_O–PtA, R_O–PtB) cleaved and one bond (R_OC–O) formed during the reaction. Two additional bonds (R_C–PtD, R_O–PtC) are constantly being adjusted during the entire reaction. When O is close to CO, the O–Pt_D bond distance also becomes smaller. While the smallest distance during the reaction is always larger than 2.7 Å, we suggest that the change in this bond will not affect the free energy change. Besides, one may assume that the change in the triple bond of CO could also influence the reaction. Hence, we performed an on-the-fly FEG analysis on this bond which shows a synchronous change in the FEG and the bond distance (see the ESI† S4 for a discussion). Our result shows that the firstly mentioned five bonds could contribute to the free energy change in the current system.

Reaction path generation and evaluation. To obtain the reaction path along the FES, we used the following three-step approach. Step 1: To find the TS on the PES, which is key to reduce the computing cost, we utilized the constrained minimization, which has been successfully applied to the TS search in various reaction systems.^49,50,57 Step 2: The refinement of the TS would be relayed by the CMD during which the chosen CV updates according to the FEG.⁵² For the CO oxidation, though the reaction is driven by different bonds at different stages, the OC–O bond plays the most significant role when it comes to the TS either on the PES or on the FES. Since the structure on the PES has already captured the major characteristic of the reaction, the optimization of the TS on the FES can be accelerated by adopting the TS on the PES as the initial structure and also avoids the region on the FES which is not merely determined by the OC–O bond. The TS at 300 K was found with very few updates of the R_OC–O by minimizing the FEG to zero as shown in Fig. 5, the length of which is 2.007 Å, a little longer than 1.988 Å of the one on the PES. Since the simulated temperature is ca. 300 K, the bond distance difference between the 0 K-PES and the 300 K-FES is very small. The increase of the OC–O bond length at the TS on the FES relative to that on the PES can be expected due to thermal vibrations. Step 3: Beginning with the structure of the TS on the FES, brute-force MD simulations were carried out with initial velocities being set to make R_C–O shorten or elongated to push it to the IS or the FS, respectively, at 300 K. These two trajectories were connected to constitute a complete reaction path from the IS to the FS. The potential energy changes and the temperature fluctuations of selected trajectories are shown in the ESI† S6. Unlike the few powerful but expensive minimum free energy path (MFEP) optimization methods, e.g. the string method^58,59 and the CMD-based discrete path optimization,^60,61 the above procedure only requires a brute-force trajectory starting from the TS. In addition, the difference between the trajectory and those of intrinsic reaction coordinate trajectories is that the former can incorporate the surface relaxation at the finite temperature.

Before carrying out the free energy decomposition, we utilized the path-CV to evaluate the effectiveness of the reaction path generated above. With the difficulty in free energy calculations along multi-CVs by enhanced sampling methods due to the unaffordable computing time, it is easy for one to utilize the path-CV to comprise plentiful CVs. From the obtained brute-force MD trajectory, we selected 30 equidistant structures with a smoothing parameter λ of 100 Å⁻² to form the path-CV including the above five bond distances as its components. To sample structures around the reaction path, harmonic potentials with the spring constant of 500 eV and 80 eV are used to restrain the path-CV s at certain values between 0 and 1 and z at 0.05, respectively. Fig. 6 shows the 2-D free energy landscape from the co-adsorbed CO and O to the adsorbed CO₂ by the US and also the projected PMF along s. The PMF from the IS to the TS is really similar to ones from the US and the CMD along the sole R_OC–O, which yields a forward free energy barrier of 0.87 eV and the backward energy of 0.92 eV. The obvious difference lies in the region s between 0.5 to 0.7. Unlike the monotonic lowering of the free energy in R_OC–O from the US simulations, there is a plateau indicating the distinct behavior between the PES and the 300 K-FES. The reason for this is the rotation from the loosely bonded CO and O with a R_OC–O of about 1.3 Å passing the TS to the stabilized bent CO₂ as the FS. The outcome 2-D free energy landscape not only demonstrates the effectiveness of the brute-force MD trajectory but also captures the rotation process.

Fig. 6 Potential of mean force for the CO oxidation obtained by the umbrella sampling along the path-CV. The 3-D free energy landscape illustrates the reaction path from a selected brute-force molecular dynamics trajectory. The free energy barriers are taken as values along the projected curve at z = 0.05, which gives 0.87 eV for the forward and 0.92 eV for the backward. The plateau between s ∈ [0.5, 0.7] reflects the slow rotation of CO₂ from the HCP site to the bridge site.

Free energy decomposition. It is clear that the path-CV made up of five bonds can generate the proper PMF. To understand the importance of these bonds, we carried out the free energy decomposition analysis, in which the CMD of the 30 reference structures were obtained. Fig. 7 shows the constrained bond distances along the reaction path and the FEGs. As can be seen, R_OC–O changes drastically along the reaction path and then keeps fluctuating around 1.3 Å, which is expected to be observed for the formation of CO₂. R_C–PtD stays around 2.0 Å which is stretched compared to the adsorption bond length. Two of three O–Pt bonds are elongated at different stages of the reaction, which reflects the bond cleavage for O to react with CO. However, the understanding obtained from the bond distance changes in the reaction is limited. Thus, we utilized the FEG from the CMD to gain insights into the free energy contributions of each bond.

Fig. 7 Bond distance changes of five bonds along the selected brute-force molecular dynamics trajectory. Three bond distances, R_OC–O, R_O–PtA and R_O–PtB, change explicitly while R_C–PtD and R_O–PtC are only slightly adjusted. During the reaction, R_OC–O varies from 3.3 Å to 1.28 Å. The distances of two O–Pt bonds, R_O–PtA and R_O–PtB, ranges from 2.0 Å to more than 3.5 Å. The distances of rest two bonds, R_C–PtD and R_O–PtC, only fluctuates around 2.0 Å.

The free energy contribution of each bond is integrated from the FEG along the reaction path. Fig. 8 illustrates the free energy changes and the FEG, respectively. The free energy changes combine the information from both the bond distance changes and the FEG. These curves show that the entire reaction can be divided into four stages: (i) the O–Pt_A bond activation with the bond increase; (ii) the repulsive interaction between CO and O with the decrease of the OC–O bond; (iii) the attractive interaction between CO and O with the further decrease of the OC–O bond; and (iv) the second O–Pt_B activation with the just formed CO₂ rotation. The individual bond contribution to the total free energy changes is apparent from the free energy changes along the bond by integrating the FEG. It can be seen from both the bond distance changes and the free energy contribution that the CO oxidation reaction is ignited by the first O–Pt_A activation. When the image number N ∈ [1, 8], O would crawl from the FCC site to the bridge site by overcoming a free energy barrier of 0.60 eV. After the O–Pt bond increases to 2.7 Å, O is close enough to interact with CO, which is reflected by the increasing free energy of 0.24 eV from ∂A/∂R_OC–O for N ∈ [8, 12]. The second stage ends with the CO getting through the repulsive region to form a roughly bonded OC–O. These results demonstrate that the major difficulty in the CO oxidation is the O self-activation (0.60 eV). The total forward free energy barrier of 0.84 eV calculated by the path-CV US can be mostly attributed to the R_O–PtA (0.60 eV) and the R_OC–O (0.24 eV). After the TS, CO and O is getting closer which is in the region N ∈ [12, 15], during which the OC–O bond becomes attractive, thus lowering the free energy by 0.45 eV. At this stage, the R_OC–O of about 1.3 Å indicates that CO₂ can be considered to have been formed. Then the second activation of the O–Pt_B bond as in N∈ [15, 30] pushes CO₂ to move from the HCP site to the bridge site. This activation is a relatively long process and lasts for almost 15 images, which lowers the free energy of 0.36 eV. During the long rotation process, the OC–O bond does not contribute considerably to the free energy changes because of its periodic stretching. As discussed above,  the OC–O bond attraction lowers the free energy by 0.45 eV larger than that of the O–Pt_B activation (0.36 eV), which will dominate the reverse process of the CO oxidation i.e. the CO₂ activation. In addition, the C–Pt_D bond and the O–Pt_C bond have no significant effect on the CO oxidation process: their free energy contributions are at a magnitude of 0.05 eV.

Fig. 8 Free energy contributions of five bonds. The reaction can be divided into four stages dominated by various bonds: (i) the first O–Pt bond activation by the O–Pt_A bond; (ii) the repulsive interaction of the OC–O bond; (iii) the attractive interaction of the OC–O bond; (iv) the second O–Pt bond by the O–Pt_B bond. Red curves are free energy changes while blue curves are respective free energy gradients.

3.4 Further verification on the collective variable selection

The free energy decomposition analysis guided us to include R_OC–O and R_O–PtA to describe the process before the TS while R_OC–O and R_O–PtB after the TS. The whole free energy decomposition pattern shows that only a few bonds contribute to the free energy changes at different stages of the reaction. The FEG of active bonds are evaluated under the conditions that other inactive bonds are constrained at the same time, which may underestimate the influence of the changes in those bonds. Although the four stages can be recognized from the free energy decomposition, it is difficult to accurately sample individual bonds for each stage. Recently, Hovan et al. suggested that this problem could be simplified by using the path-CV made up of weighted CVs.⁶² The method requires the weight coefficients of CV along the reaction progress to incorporate a multidimensional CV into one weighted path-CV. Despite its potential in capturing CV significance in driving the dynamics, it could be computationally intensive for obtaining these weights, and thus expensive for the AIMD. For simplicity, we only divided this elementary reaction into two parts which are separated by the TS, in which the CMD was performed along the CV by combining the indispensable bonds discussed above.

For the process from the IS to the TS, we took the linear combination of two bond distances S(R_OC–O, R_O–PtA), which incorporates the influence of these two bonds. R_O–PtA is used to describe the process of the O–Pt_A bond activation while R_OC–O is used for the interaction between CO and O near the TS. Fig. 9 shows the PMF along the selected CV. In this PMF, the near zero FEGs indicate that the IS is at S = 1.3 Å and the TS is at S = −0.9 Å. It gives a forward free energy barrier of 0.87 eV, which is similar to the ones obtained for other CVs. The output FEG decreases gradually with the increase of S(R_OC–O, R_O–PtA) and reduces mildly to zero, which is much smoother than those generated along the sole R_OC–O. Although the interaction between the adsorbates and the surface is reflected by samples implicitly no matter what CV is used, the CV should still be chosen carefully to better describe the surface relaxation, which is the major difference from gas reactions. In the CMD along the linear combination CV S(R_OC–O, R_O–PtA), there are no constraints applied on its components. Unlike the OC–O bond with a fixed length which will restrict the relaxation of the surface in the CMD along R_OC–O, samples from the CMD along S(R_OC–O, R_O–PtA) on the much more flexible surface will be collected in the simulations.

Fig. 9 Potential of mean force along S(R_OC–O, R_O–PtA) = R_OC–O − R_O–PtA. The sampling is from −1.0 Å to 1.2 Å with an interval of 0.1 Å. The smooth free energy gradient change yields a forward free energy barrier of 0.87 eV.

The process after the TS becomes a little complicated due to the rotation process of CO₂ from the HCP site to the bridge site. We utilized two CVs, R_OC–O and S(R_O–PtB, R_O–PtC). R_OC–O was used to describe the attractive interaction between C and O after the TS and the changes in S(R_O–PtB, R_O–PtC) can capture the rotation process. The resulting PMF and the FEGs are shown in Fig. 10. Checked by the CMD in advance, there are two metastable states where FEGs equal zero on this 2-D (S,R) FES, namely the TS at (0.0, 2.0) and the FS at (1.20, 1.28). The FEG indicates that the rotation, i.e. the increase in S(R_O–PtB, R_O–PtC), may encounter a small energy barrier when R_OC–O is larger than 1.4 Å (see ESI† S6). In addition, a few CMD simulations show that CO₂ on the HCP site can desorb readily when R_C–O is 1.3 Å. It is clear from our results that the OC–O complex is still bonded with the surface until R_OC–O is 1.4 Å. Then it goes through the rotation which is represented by the changes in S(R_O–PtB, R_O–PtC) from 0.0 Å to 1.2 Å while maintaining R_OC–O, and finally forms the adsorbed CO₂ with R_OC–O of 1.28 Å. The free energy changes of the above process were computed from the integral as

Fig. 10 Projected 2-D free energy surface (S,R) spanned by R_OC–O and S(R_O–PtB, R_O–PtC). The potential of the mean force along R_OC–O is separated into two parts by the CO₂ rotation indicated by the change in S(R_O–PtB, R_O–PtC) ranging from 0.0 Å to 1.2 Å. The rotation happens when R_OC–O is 1.4 Å due to the gradually decreasing free energy barrier along S(R_O–PtB, R_O–PtC). The free energy change in R_OC–O is −0.54 eV and −0.18 eV before and after the rotation. The free energy change in S(R_O–PtB, R_O–PtC) is −0.30 eV. Thus, the total backward free energy barrier is 1.02 eV.

The free energy change due to S(R_O–PtB, R_O–PtC) is −0.30 eV, which is similar to the C–Pt_B bond free energy contribution of −0.36 eV from the free energy decomposition analysis in Section 3.3. The total free energy contribution in R_OC–O is −0.72 eV, which can be separated into −0.54 eV and −0.18 eV due to the change in R_OC–O before and past the rotation. Thus, the backward free energy barrier is summed as 1.02 eV (0.30 eV + 0.72 eV). Although the value of the free energy is similar to above PMFs using other CVs, the 2-D FES is proved to be better at describing the rotation process. The standard errors of FEGs and information about other sampled points can be found in the ESI† S7.

Having thoroughly tested various choices of CVs and enhanced sampling methods, we are in the position to make general discussions. Since most free energy calculation methods are path-based and require a predefined collective variable, the protocol we devised above is able to guide the CV selection by quantifying the free energy contributions of various bonds especially for surface reactions. For the CO oxidation on the Pt(111) surface at 300 K, the free energy changes along various CVs are very similar. However, the interpretability of the reaction largely depends on the CV used in the free energy calculations. To determine free energy barriers of this reaction swiftly by the AIMD, we suggest to sample R_OC–O with the CMD. This choice gives a fair estimation of free energy changes if the collected structures are taken care of when the simulation comes close to the FS. In addition, the CMD assures the existence of local minima and saddle points on the FES by checking the FEG. On the other hand, the US requires much care in determining the CV sampling range with a prior knowledge of the IS and the FS. To effectively obtain the reaction path, our results suggest the use of S(R_OC–O, R_O–PtA) before the TS and R_OC–O in tandem with S(R_O–PtB, R_O–PtC) after the TS. These CVs will be much helpful for investigating the different reaction stages when the catalyst surface varies, which will provide fruitful insights into the tuning of the catalytic process. For the reaction with CVs entangled, our work shows that a path-CV is needed, which may be computationally intensive.

4. Conclusion

In this work, we have revisited the model catalytic reaction, CO oxidation on the Pt(111) surface, for heterogeneous catalysis at 300 K with enhanced sampling methods including the US and the CMD. Various choices of CVs for this reaction were tested, namely bond distances, their linear combinations and the path-CV. We found that the PMF along R_OC–O is not enough for understanding the whole reaction although it generates similar free energy values to that of others. To further select essential CVs for the free energy calculations, we proposed a method to quantify bond contributions on the reaction path extracted from brute-force MD trajectories initiated from the TS on the FES at the finite temperature. Before carrying out the free energy decomposition, the path-CV was utilized to construct the free energy landscape of this reaction by taking five bonds (R_OC–O, R_C–PtD, R_O–PtA, R_O–PtB, R_O–PtC) into consideration, which assesses the reaction path from a selected brute-force MD trajectory. As a significant step in the CV selection, the free energy decomposition analysis showed that R_OC–O and R_O–PtA are responsible for the forward free energy barrier while R_OC–O and R_O–PtB for the backward energy barrier. The further assessment on S(R_OC–O, R_O–PtA) for the first half reaction while R_OC–O and S(R_O–PtB, R_O–PtC) for the second yields fair estimations of free energy barriers and captures the various reaction stages. Regarding the free energy contributions, the O–Pt_A bond is the determined step for the CO oxidation while the OC–O bond is for the reverse process. These new insights into the CO oxidation could also advance the free energy calculation of carbon oxides reactions on noble metal catalysts. Above all, this reaction path based on free energy decomposition for the CV selection is remarkably feasible for the field of heterogeneous catalysis because it all starts from the TS on the PES.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We are grateful for the computational support from the UK national high performance computing service, ARCHER, for which access was obtained via the UKCP consortium and funded by EPSRC grant ref. EP/P022561/1. We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/P020194/1). We thank Chenxi Guo, Peter Price and Chenggong Hui for fruitful discussions. Jiayan Xu and Hao Huang thanks the Queen's University Belfast and China Scholarship Council for financial support.

Notes and references

1. X.-M. Cao, R. Burch, C. Hardacre and P. Hu, An Understanding of Chemoselective Hydrogenation on Crotonaldehyde over Pt(111) in the Free Energy Landscape: The Microkinetics Study based on First-principles Calculations, Catal. Today, 2011, 165(1), 71–79.
2. Z. Wang, X. M. Cao, J. Zhu and P. Hu, Activity and Coke Formation of Nickel and Nickel Carbide in Dry Reforming: A Deactivation Scheme from Density Functional Theory, J. Catal., 2014, 311, 469–480.
3. X. Sun, P. Wang, Z. Shao, X. Cao and P. Hu, A First-principles Microkinetic Study on the Hydrogenation of Carbon Dioxide over Cu(211) in the Presence of Water, Sci. China: Chem., 2019, 62(12), 1686–1697.
4. C. J. Cramer, Essentials of Computational Chemistry: Theories and Models, John Wiley & Sons Ltd, 2nd edn, 2004.
5. L. Foppa, M. Iannuzzi, C. Copéret and A. Comas-Vives, CO Methanation on Ruthenium Flat and Stepped Surfaces: Key Role of H-transfers and Entropy Revealed by Ab Initio Molecular Dynamics, J. Catal., 2019, 371, 270–275.
6. X. Zeng, Z. Qiu, P. Li, Z. Li and J. Yang, Steric Hindrance Effect in High-Temperature Reactions, CCS Chem., 2020, 460–467.
7. T. Cheng, H. Xiao and W. A. Goddard III, Free-Energy Barriers and Reaction Mechanisms for the Electrochemical Reduction of CO on the Cu(100) Surface, Including Multiple Layers of Explicit Solvent at pH 0, J. Phys. Chem. Lett., 2015, 6(23), 4767–4773.
8. T. Cheng, H. Xiao and W. A. Goddard III, Full Atomistic Reaction Mechanism with Kinetics for CO Reduction on Cu(100) from Ab Initio Molecular Dynamics Free-energy Calculations at 298 K, Proc. Natl. Acad. Sci. U. S. A., 2017, 114(8), 1795–1800.
9. K. Klyukin and V. Alexandrov, CO₂ Adsorption and Reactivity on Rutile TiO₂(110) in Water: An Ab Initio Molecular Dynamics Study, J. Phys. Chem. C, 2017, 121(19), 10476–10483.
10. D. Wang, T. Sheng, J. Chen, H.-F. Wang and P. Hu, Identifying the Key Obstacle in Photocatalytic Oxygen Evolution on Rutile TiO₂, Nat. Catal., 2018, 1(4), 291–299.
11. J. J. Sun and J. Cheng, Solid-to-liquid Phase Transitions of Sub-nanometer Clusters Enhance Chemical Transformation, Nat. Commun., 2019, 10(1), 5400.
12. H. Guo, P. Sautet and A. N. Alexandrova, Reagent-Triggered Isomerization of Fluxional Cluster Catalyst via Dynamic Coupling, J. Phys. Chem. Lett., 2020, 11(8), 3089–3094.
13. C. Guo, Z. Wang, D. Wang, H.-F. Wang and P. Hu, First-Principles Determination of CO Adsorption and Desorption on Pt(111) in the Free Energy Landscape, J. Phys. Chem. C, 2018, 122(37), 21478–21483.
14. R. Réocreux, C. Michel, P. Fleurat-Lessard, P. Sautet and S. N. Steinmann, Evaluating Thermal Corrections for Adsorption Processes at the Metal/Gas Interface, J. Phys. Chem. C, 2019, 123(47), 28828–28835.
15. E. A. Carter, G. Ciccotti, J. T. Hynes and R. Kapral, Constrained Reaction Coordinate Dynamics for the Simulation of Rare Events, Chem. Phys. Lett., 1989, 156(5), 472–477.
16. G. M. Torrie and J. P. Valleau, Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling, J. Comput. Phys., 1977, 23(2), 187–199.
17. J. Kästner, Umbrella Sampling, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1(6), 932–942.
18. A. Laio and M. Parrinello, Escaping Free-energy Minima, Proc. Natl. Acad. Sci. U. S. A., 2002, 99(20), 12562–12566.
19. A. Barducci, M. Bonomi and M. Parrinello, Metadynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1(5), 826–843.
20. E. Darve and A. Pohorille, Calculating Free Energies using Average Force, J. Chem. Phys., 2001, 115(20), 9169–9183.
21. C. Chipot and A. Pohorille, Free Energy Calculations: Theory and Applications in Chemistry and Biology, Springer-Verlag, Berlin Heidelberg, 2007.
22. E. Darve, D. Rodriguez-Gomez and A. Pohorille, Adaptive Biasing Force Method for Scalar and Vector Free Energy Calculations, J. Chem. Phys., 2008, 128(14), 144120.
23. P. P. Kumar, A. G. Kalinichev and R. J. Kirkpatrick, Dissociation of Carbonic Acid: Gas Phase Energetics and Mechanism from Ab Initio Metadynamics Simulations, J. Chem. Phys., 2007, 126(20), 204315.
24. D. Branduardi, F. L. Gervasio and M. Parrinello, From A to B in Free Energy Space, J. Chem. Phys., 2007, 126(5), 054103.
25. D. Branduardi, M. De Vivo, N. Rega, V. Barone and A. Cavalli, Methyl Phosphate Dianion Hydrolysis in Solution Characterized by Path Collective Variables Coupled with DFT-Based Enhanced Sampling Simulations, J. Chem. Theory Comput., 2011, 7(3), 539–543.
26. D. Mendels, G. Piccini and M. Parrinello, Collective Variables from Local Fluctuations, J. Phys. Chem. Lett., 2018, 9(11), 2776–2781.
27. S. R. Hare, L. A. Bratholm, D. R. Glowacki and B. K. Carpenter, Low Dimensional Representations along Intrinsic Reaction Coordinates and Molecular Dynamics Trajectories using Interatomic Distance Matrices, Chem. Sci., 2019, 10(43), 9954–9968.
28. K. Haas and J. W. Chu, Decomposition of Energy and Free Energy Changes by Following the Flow of Work along Reaction Path, J. Chem. Phys., 2009, 131(14), 144105.
29. A. Alavi, P. Hu, T. Deutsch, P. L. Silvestrelli and J. Hutter, CO Oxidation on Pt(111): An Ab Initio Density Functional Theory Study, Phys. Rev. Lett., 1998, 80(16), 3650–3653.
30. I. Chorkendorff and J. W. Niemantsverdriet, Concepts of Modern Catalysis and Kinetics, Wiley-VCH Verlag GmbH & Co. KGaA, 2003.
31. Z. Wang and P. Hu, Identifying the General Trend of Activity of Non-stoichiometric Metal Oxide Phases for CO Oxidation on Pd(111), Sci. China Chem., 2019, 62(6), 784–789.
32. K. Bleakley and P. Hu, A Density Functional Theory Study of the Interaction between CO and O on a Pt Surface-: CO/Pt(111), O/Pt(111), and CO/O/Pt(111), J. Am. Chem. Soc., 1999, 121(33), 7644–7652.
33. C. J. Zhang and P. Hu, Why Must Oxygen Atoms Be Activated from Hollow Sites to Bridge Sites in Catalytic CO Oxidation?, J. Am. Chem. Soc., 2000, 122(9), 2134–2135.
34. L. Zhou, A. Kandratsenka, C. T. Campbell, A. M. Wodtke and H. Guo, Origin of Thermal and Hyperthermal CO₂ from CO Oxidation on Pt Surfaces: The Role of Post-Transition-State Dynamics, Active Sites, and Chemisorbed CO₂, Angew. Chem., Int. Ed., 2019, 58(21), 6916–6920.
35. Q. Wu, L. Zhou and H. Guo, Steric Effects in CO Oxidation on Pt(111) by Impinging Oxygen Atoms Lead to an Exclusive Hot Atom Mechanism, J. Phys. Chem. C, 2019, 123(16), 10509–10516.
36. G. Kresse and J. Furthmüller, Efficiency of Ab-initio Total Energy Calculations for Metals and Semiconductors using a Plane-wave Basis Set, Comput. Mater. Sci., 1996, 6(1), 5–50.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77(18), 3865–3868.
38. G. Kresse and D. Joubert, From Ultrasoft Pseudopotentials to the Projector Augmented-wave Method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59(3), 1758–1775.
39. M. Methfessel and A. T. Paxton, High-precision Sampling for Brillouin-zone Integration in Metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40(6), 3616–3621.
40. H. J. Monkhorst and J. D. Pack, Special Points for Brillouin-zone Integrations, Phys. Rev. B: Solid State, 1976, 13(12), 5188–5192.
41. W. P. Davey, Precision Measurements of the Lattice Constants of Twelve Common Metals, Phys. Rev., 1925, 25(6), 753–761.
42. H. A. Posch, W. G. Hoover and F. J. Vesely, Canonical Dynamics of the Nosé Oscillator: Stability, Order, and Chaos, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 33(6), 4253–4265.
43. T. Bucko, Ab Initio Calculations of Free-energy Reaction Barriers, J. Phys.: Condens. Matter, 2008, 20(6), 064211.
44. H. Flyvbjerg and H. G. Petersen, Error Estimates on Averages of Correlated Data, J. Chem. Phys., 1989, 91(1), 461–466.
45. Spencer, J. Pyblock, https://github.com/jsspencer/pyblock/blob/master/docs/index.rst, accessed Feburary 28.
46. S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman and J. M. Rosenberg, The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method, J. Comput. Chem., 1992, 13(8), 1011–1021.
47. Grossfield, A. WHAM: the Weighted Histogram Analysis Method Version 2.0.9.1, http://membrane.urmc.rochester.edu/wordpress/?page_id=126, accessed Feburary 20.
48. K. Zinovjev, S. Martí and I. Tuñón, A Collective Coordinate to Obtain Free Energy Profiles for Complex Reactions in Condensed Phases, J. Chem. Theory Comput., 2012, 8(5), 1795–1801.
49. A. Michaelides and P. Hu, Catalytic Water Formation on Platinum: A First-Principles Study, J. Am. Chem. Soc., 2001, 123(18), 4235–4242.
50. Z.-P. Liu and P. Hu, General Rules for Predicting Where a Catalytic Reaction Should Occur on Metal Surfaces: A Density Functional Theory Study of C-H and C-O Bond Breaking/Making on Flat, Stepped, and Kinked Metal Surfaces, J. Am. Chem. Soc., 2003, 125(7), 1958–1967.
51. G. Henkelman and H. Jónsson, Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points, J. Chem. Phys., 2000, 113(22), 9978.
52. S.-Y. Yang, I. Hristov, P. Fleurat-Lessard and T. Ziegler, Optimizing the Structures of Minimum and Transition State on the Free Energy Surface, J. Phys. Chem. A, 2005, 109(1), 197–204.
53. X.-M. Cao, Z.-J. Shao and P. Hu, A Fast Species Redistribution Approach to Accelerate the Kinetic Monte Carlo Simulation for Heterogeneous Catalysis, Phys. Chem. Chem. Phys., 2020, 22(14), 7348–7364.
54. G. Collinge, S. F. Yuk, M.-T. Nguyen, M.-S. Lee, V.-A. Glezakou and R. Rousseau, Effect of Collective Dynamics and Anharmonicity on Entropy in Heterogenous Catalysis: Building the Case for Advanced Molecular Simulations, ACS Catal., 2020, 10, 9236–9260.
55. M. del Cueto, X. Zhou, L. Zhou, Y. Zhang, B. Jiang and H. Guo, New Perspectives on CO₂–Pt(111) Interaction with a High-Dimensional Neural Network Potential Energy Surface, J. Phys. Chem. C, 2020, 124(9), 5174–5181.
56. P. Fleurat-Lessard and T. Ziegler, Tracing the Minimum-energy Path on the Free-energy Surface, J. Chem. Phys., 2005, 123(8), 084101.
57. J. Xu, X.-M. Cao and P. Hu, Improved Prediction for the Methane Activation Mechanism on Rutile Metal Oxides by a Machine Learning Model with Geometrical Descriptors, J. Phys. Chem. C, 2019, 123(47), 28802–28810.
58. L. Maragliano, A. Fischer, E. Vanden-Eijnden and G. Ciccotti, String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces, J. Chem. Phys., 2006, 125(2), 024106.
59. A. C. Pan, D. Sezer and B. Roux, Finding Transition Pathways Using the String Method with Swarms of Trajectories, J. Phys. Chem. B, 2008, 112(11), 3432–3440.
60. C. Chen, Y. Huang and Y. Xiao, Free-energy Calculations along a High-dimensional Fragmented Path with Constrained Dynamics, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86(3), 031901.
61. C. Chen, Y. Huang, X. Ji and Y. Xiao, Efficiently Finding the Minimum Free Energy Path from Steepest Descent Path, J. Chem. Phys., 2013, 138(16), 164122.
62. L. Hovan, F. Comitani and F. L. Gervasio, Defining an Optimal Metric for the Path Collective Variables, J. Chem. Theory Comput., 2019, 15(1), 25–32.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp03852k

---