Van der Waals density functional study of formic acid adsorption and decomposition on Cu(111)

The needs to discover renewable energy sources and to reduce fossil fuel consumption have driven the research and development of fuel cell technology. One of the leading developments of fuel cells is proton-exchange membrane fuel cells (PEMFCs).^1–3 The PEMFCs convert chemical energy into electrical energy by the electrochemical reaction from hydrogen and oxygen to water. While hydrogen gas can be used directly in PEMFCs, the storage of compressed hydrogen suffers from a loss of hydrogen, safety issues, and low volumetric capacity. Therefore, alternative material for hydrogen storage is required to overcome those problems.

Recently, formic acid (HCOOH) has been considered as a potential material for hydrogen storage. The density of HCOOH is 1.22 g/cm³ with volumetric capacity of 53 gH₂/l at standard temperature and pressure (STP), which surpasses that of other storage materials.^4,5 Moreover, storage of HCOOH is easier and safer compared with other hydrogen storage material because HCOOH is a nonflammable liquid at STP. Hydrogen stored in HCOOH can be released by the decomposition process. HCOOH can be decomposed on metal surfaces through the dehydrogenation process into CO₂ and H₂, or the dehydration process into CO and H₂O.^1–3,6 As a hydrogen source, a catalyst with high selectivity to the dehydrogenation process is preferable.^7–11 Therefore, elucidation of HCOOH decomposition on metal surfaces has been attracting substantial attention.

Decomposition of HCOOH has been studied experimentally on several transition metal catalysts, e.g., Cu,^12–14 Ni,¹⁵ Au,¹⁶ Co,¹¹ Pd,¹⁷ and Pt.¹⁸ Cu catalysts have been reported to selectively decompose HCOOH through the dehydrogenation process. Bowker et al.,¹⁴ Hayden et al.,¹³ and Sexton¹² observed decomposition of HCOOH on the Cu(110) and Cu(100) surfaces when exposed at 300 K. On the other hand, the decomposition of HCOOH does not occur on Cu(111) when exposed to the gas phase HCOOH at room temperature.^19–23 On the Cu(111) surface, HCOOH was reported to form polymeric structures when exposed to the gas phase HCOOH at low temperatures, and then the decomposition occurs by heating those polymeric HCOOH.^21–23 

A number of theoretical studies have been done to investigate HCOOH adsorption and decomposition on several transition metal surfaces.^{11,16–18,24–27} In these studies, generalized gradient approximations (GGAs) to the exchange-correlation functional were used to describe the HCOOH adsorption and decomposition processes on the surface. Those functionals do not describe the van der Waals (vdW) interaction or dispersion forces properly, and as a consequence, they fail to describe weakly adsorbed systems. HCOOH is weakly adsorbed on Cu(111) and the vdW correction is indispensable. In the present study, we have investigated the adsorption and decomposition mechanisms of monomeric HCOOH on Cu(111) by using density functionals (DFs) that account for the vdW forces and compared them with room temperature experimental results.

Our calculations were performed using the STATE (Simulation Tool for Atom TEchnology) code, which has been used for formate adsorption,²⁸ formate hydrogenation,²⁹ CO₂ hydrogenation,³⁰ CO₂ dissociation,³¹ CO₂ adsorption,³² and formate decomposition.³² Here, we compared results obtained using the Perdew-Burke-Ernzerhof (PBE)³³ functional with those using the van der Waals density functionals (vdW-DFs), i.e., optB86b-vdW³⁴ and rev-vdW-DF2³⁵ functionals. We also included the dispersion correction proposed by Grimme with PBE (PBE-D2).³⁶ The implementation of the self-consistent vdW-DF^37–39 in the STATE code is described in Ref. 40. The electron-ion interaction is described using the Vanderbilt’s ultrasoft pseudopotentials.⁴¹ Wave functions and charge densities are expanded in terms of a plane wave basis set with cutoff energies of 36 and 400 Ry, respectively. The Cu(111) surface was modeled using a repeated slab model with three atomic layer-thickness with the bottom most atomic layer fixed to its bulk configuration. By increasing the layer thickness from three to six, we confirm that the error in adsorption and activation energies due to the finite slab thickness is within 0.09 eV, as shown in Tables S1, S2, and S5 in the supplementary material. We used a (3 × 3) surface unit cell of Cu(111) and the surface Brillouin zone was sampled using a 4 × 4 × 1 uniform k-points mesh. We set the convergence criterion for forces acting on atoms in the structural optimization to 5.14 × 10⁻² eV/Å (10⁻³ hartree/bohr).

The HCOOH adsorption energy E_ads on the surface is defined by

where E_tot, E_surf, and E_HCOOH represent the total energies of the adsorbed system, the isolated clean surface, and the isolated molecule in gas phase, respectively. We employed the climbing image nudged elastic band (CI-NEB) method^42,43 to evaluate the activation energies for diffusion, rotation, and decomposition of HCOOH on Cu(111). A transition state was identified by finding a replica with the highest energy in between the initial and final states with the maximum force below 5.14 × 10⁻² eV/Å and was further confirmed by obtaining a single imaginary frequency in the vibrational frequency calculations. In the adsorption and activation energies calculations, we included the zero-point energy (ZPE) correction, which was calculated by assuming a quantum harmonic oscillator as the sum of contributions from all vibrational modes of the adsorbed system.

We performed the reaction rate analysis for desorption and decomposition processes. The desorption rate is calculated based on the assumption that the desorption rate in the vacuum condition is the same as that in the equilibrium condition between gas phase molecules with temperature (T) and pressure (P) and adsorbed molecules on the surface.^44–46 In the equilibrium condition with gas phase molecules, the adsorption and desorption rates per area are equal. We can express the rate of desorption as the flux (F) of impinging molecules on the surface per area and time, multiplied by the sticking probability S(θ, T)

where P_eq, m, k_B, and T represent equilibrium pressure, molecule molar mass, Boltzmann’s constant, and temperature, respectively. The equilibrium pressure can be derived from the condition that the chemical potentials of gas phase and adsorbed phase are equal. The chemical potential of the gas phase molecule (μ_g) can be expressed by the ideal gas model. On the other hand, the chemical potential for the adsorbed phase can be expressed by that of the 2D lattice gas model (⁠$μadslattice$⁠) or that of the 2D van der Waals gas model (⁠$μads2DvdW$⁠). The lattice gas model assumes that the adsorbates are chemically bonded to the adsorption sites of a surface.⁴⁴ On the other hand, the 2D van der Waals gas model assumes that the adsorbates can move freely on the surface and interact with each other via the pairwise hard-core repulsive potential. The nature of the adsorbate phase is determined by the relationship between the activation energy for diffusion (E_diff) and the thermal energy (k_BT). Doll and Steele showed that using a simple model, the transition between the two models takes place at T_c ≃ E_diff/5k_B, namely, for temperatures lower than critical temperature (T_c), the lattice gas model is more suitable to describe the adsorbate phase, while the 2D van der Waals is more suitable when temperatures are higher than T_c.⁴⁷ Those chemical potentials are defined by the following equations:

where h is the Planck constant, m is the molar mass of HCOOH, A_s is the surface area per substrate surface atom, and θ represents the coverage. q_v–g, q_v–ads, q_rg, and q_r–ads represent the gas phase vibration, adsorption phase vibration, gas phase rotation, and adsorption phase rotation partition functions, respectively, which are defined by

where ω_j represents the vibrational frequency, I_A, I_B, and I_C are the principal moments of inertia of HCOOH, and σ is the symmetry number. Here, σ is equal to 1 for gas phase HCOOH molecule.

By equating the chemical potential of the gas phase and adsorbed phase, we obtain the equilibrium pressure as a function of coverage, temperature, partition functions, and adsorption energy. Then, we substitute the equilibrium pressure into Eq. (2). Finally, the kinetic constants for the desorption, which correspond to the desorption rates per surface adsorbate for two models, namely, the 2D lattice gas model and 2D van der Waals gas model, are expressed as

where n_s is the number of surface atoms per unit area [for Cu(111), n_s = 1.74 × 10¹⁹ m⁻²]. E_des represents the activation energy for the desorption. In this work, we assume that there is no additional barrier except the energy increase from the adsorbed molecule to the gas phase, i.e., E_des = −E_ads. In both cases, we assume that the sticking probability (S(θ, T)) is equal to unity.

The decomposition rate is defined by using the transition state theory based on the Eyring formulation with an assumption that the activated complexes are in quasi-equilibrium with reactants.^48–50 The rate can be written proportional to the concentration of the complexes multiplied by the kinetic constant. The concentration is expressed by molar partition function, which contains the translation, rotation, and vibration partition functions. Thus, the kinetic constant for the decomposition, which corresponds to the decomposition rate per surface adsorbate can be expressed as

where E_dec is the energy difference between the transition state and the adsorbed HCOOH on the surface (E_dec = E_TS − E_ads). The Q_TS and Q_ads represent the product of the translational, rotational, and vibrational partition functions of the transition state and adsorbed phase, respectively.

We considered one parallel and two perpendicular configurations of monomeric HCOOH on the Cu(111) surface with several adsorption sites. The most stable adsorption site of the parallel configuration is shown in Fig. 1(a), in which HCOOH adsorbs with the carbon (C) atom located at the fcc hollow site and both oxygen (O) atoms are on top sites. This parallel configuration is the same as that reported in a previous theoretical study by Chutia et al.⁵¹ In the case of perpendicular configuration, we considered two different geometries, namely, CH-perpendicular and OH-perpendicular, as shown in Figs. 1(b) and 1(c), respectively. One of O atoms of HCOOH in CH-perpendicular configuration is located on top of a Cu atom and its CH bond pointing toward the surface as mention in the previous studies.^27,52 On the other hand, the HCOOH in the OH-perpendicular configuration [Fig. 1(c)] is bonded to a surface Cu atom through an O atom with its CH bond pointing away from the surface, and the OH bond pointing toward the surface. The OH-perpendicular configuration has been reported by previous theoretical studies on Cu,^26,51,53 Ni,¹⁵ Au,¹⁶ Pd,^17,27 and Pt¹⁸ surfaces.

As seen in Table I, E_ads of the OH-perpendicular configuration using PBE functional is −0.24 eV, which agrees with previous theoretical results,^24,26,53,54 but much weaker than experimentally reported values of −0.55 eV at low coverage.²³ On the other hand, E_ads calculated with vdW inclusive functionals are about −0.50 to −0.60 eV, in good agreement with the experimental result, indicating that the vdW correction is indispensable to describe the present system. We also found that the OH-perpendicular configuration is more stable than the other two configurations by 0.20–0.25 eV.

In this section, we investigate the activation energies of diffusion and rotation processes for all the three adsorption configurations of HCOOH on Cu(111) by using the rev-vdW-DF2 energy functional. The diffusion and rotation processes are calculated to clarify the behavior of the adsorbed monomeric HCOOH on the surface at finite temperature. The diffusion process starts from the most stable adsorption site as the initial state (IS) to the neighboring most stable adsorption site as the final state (FS). On the other hand, the rotation process is simulated by rotating the adsorbed HCOOH at the equilibrium adsorption site. Then, the diffusion and rotation activation energies were obtained by taking the energy difference between the initial state (IS) and the transition state (TS), and the results are summarized in Table II.

Figure 2 shows the diffusion and rotation processes for HCOOH on Cu(111) in the parallel configuration. We found that the activation energy of the diffusion process for parallel configuration is 0.03 eV. In the case of rotation, we considered the same IS by gradually rotating the adsorbed molecule from the IS up to 120° within the molecular plane, and the activation energy is calculated to be 0.02 eV. These small activation energies for diffusion and rotation indicate that HCOOH in the parallel configuration is easy to diffuse and rotate on the surface.

The diffusion and rotation processes for HCOOH in the CH-perpendicular configuration are shown in Fig. 3. We obtained that the activation energies for diffusion and rotation are 0.06 eV and 0.01 eV, respectively. Therefore, HCOOH in the CH-perpendicular configuration is also easy to diffuse and rotate on the surface.

Figure 4 shows the diffusion and rotation processes for HCOOH on Cu(111) in the OH-perpendicular configuration. The diffusion process starts from the most stable adsorption site and ends at the most stable nearest adsorption site, as shown in Fig. 4(a). We found that the activation energy of the diffusion process for the OH-perpendicular configuration is 0.19 eV, which is higher than those of the other configurations. On the other hand, in the case of rotation, HCOOH is rotated by 60° around the O–Cu bond [Fig. 4(b)] or rotated by 180° around the O–H bond [Fig. 4(c)]. We found that the activation energies for the rotation are 0.05 eV and 0.19 eV for (b) and (c) cases, respectively. In the case of (b), rotation around the O–Cu bond gives lower activation energy compared with that for the (c) case, where the O–Cu bond is broken. Therefore, the O–Cu bond between HCOOH and the Cu(111) surface restricts the movement of adsorbed HCOOH in the OH-perpendicular configuration.

In Sec. II, we introduced two different chemical potentials for the adsorbed phase, namely, 2D lattice gas model (⁠$μadslattice$⁠) and 2D van der Waals gas model (⁠$μads2DvdW$⁠). According to Doll and Steele,⁴⁷ the transition between the two models takes place at T_c ≃ E_diff/5k_B. We obtain that the T_c for parallel, CH-perpendicular, and OH-perpendicular configurations are 70 K, 139 K, and 441 K, respectively. For the OH-perpendicular configuration, T_c is above room temperature, while for the parallel and CH-perpendicular configurations, T_c is far below room temperature. Therefore, the lattice gas model is more suitable to calculate the chemical potential of adsorbed monomeric HCOOH in the OH-perpendicular configuration, while the 2D van der Waals gas model is suitable for the parallel and CH-perpendicular configurations.

We investigated the coverage (θ) of parallel, CH-perpendicular, and OH-perpendicular configurations as a function of temperature by using the rev-vdW-DF2 energy functional. The coverages of HCOOH in the three configurations are calculated by equating the chemical potentials between an ideal gas and the adsorbed phase described by Eqs. (3)–(5). Here, we assumed that the Cu(111) surface is exposed to the gas phase HCOOH at the equilibrium pressure of P_eq = 1 × 10⁻⁶ Torr, as reported by Nakano et al.²⁰ 

Figure 5 shows θ of the adsorbed HCOOH in all configurations. Below ∼280 K, the order of coverages of the three adsorption configurations is θ^∥HCOOH < θ^CH-⊥HCOOH < θ^OH-⊥HCOOH. This order comes from the adsorption energies $Eads∥HCOOH>EadsCH−⊥HCOOH>EadsOH−⊥HCOOH$ (see Table I), namely, the OH-perpendicular configuration, which is the most stable, dominates below ∼280 K. Above 300 K, the order of the coverages is completely reversed, namely, θ^∥HCOOH > θ^CH-⊥HCOOH > θ^OH-⊥HCOOH. This re-order is due to the entropic effect in the adsorbed molecule. Since the diffusion of OH-perpendicular is restricted, the entropic contribution becomes lower compared with the other two configurations, leading to the suppression of the OH-perpendicular configuration near room temperature. On the other hand, the increase of θ for the parallel configuration compared with the CH-perpendicular configuration is mainly due to the vibrational entropy, especially the metal-molecule stretching mode (see Table S4 in the supplementary material).

Note that the coverages of the three molecularly adsorbed states are rather small at room temperature, θ ≃ 10⁻⁸. At P_eq = 1 × 10⁻⁶ Torr, in average, nearly one gas phase molecule impinges to the surface per second per surface molecule. This indicates that at room temperature, a HCOOH molecule adsorbed on the Cu(111) surface stays on the surface temporally, and it is desorbed or decomposes within ∼10⁻⁶ s. We will discuss the desorption and decomposition rates in more detail later.

The reaction path and activation energies for the decomposition (E_dec) of the OH-perpendicular configuration on the Cu(111) surface are shown in Fig. 6 and Table III, respectively. It had been suggested that the formate (HCOO) mediated pathway is the most preferable decomposition pathway for the HCOOH on Cu(111) compared to carboxyl (COOH) and formyl (HCO) mediated pathways.^{24–26,55,56} The decomposition takes place by elongating the O–H bond in HCOOH. Here, we assumed that a bidentate formate (bi-HCOO) and hydrogen atom are formed as final product of the decomposition process [shown in Fig. 6(a)]. This final state configuration was also reported in previous theoretical studies.^26,51,53 The calculated activation energies of decomposition (E_dec’s) with the ZPE correction are 0.37 eV, 0.31 eV, 0.32 eV, and 0.30 eV for PBE, PBE-D2, rev-vdW-DF2, and optB86b-vdW functionals, respectively, indicating that the E_dec are almost independent among the considered functionals in this work. The E_dec for the PBE functional (0.37 eV) is slightly lower than the previously calculated value by 0.1 eV.²⁶ Moreover, the calculated E_dec of OH-perpendicular configuration using PBE is 0.24 eV larger than the activation energy for its desorption (E_des). On the other hand, PBE-D2 and vdW-DFs produce lower E_dec than E_des by (0.16–0.27 eV).

The decomposition mechanism from the parallel configuration is shown in Fig. 7(a). In this case, we considered the decomposition through HCOO-mediated pathway, the same as the OH-perpendicular configuration. In this path, transformation from the parallel configuration to the OH-perpendicular perpendicular configuration is inevitable. We notice that there is almost no activation energy for this transformation and the intermediate state (IM) becomes more stable than the parallel configuration by 0.27 eV. Thus, the effective activation energy of decomposition from the parallel configuration is calculated by taking the energy difference between the IM and transition state (TS), resulting in the same value as that from the OH-perpendicular configuration, namely, 0.31 eV (Table IV).

In the case of the CH-perpendicular configuration, Wang et al.²⁷ proposed that the decomposition through COOH-mediated pathway is more preferable rather than the other decomposition processes on Pd(111), and here, we consider a similar process as shown in Fig. 7(b). We obtain that the E_dec with ZPE correction for the CH-perpendicular configuration is 0.69 eV, which is higher than those of other configurations.

Therefore, we conclude that the decomposition starting from the OH-perpendicular configuration to the bidentate HCOO is the most dominant path among the three decomposition paths considered and accordingly, in the decomposition rate analysis which will be shown in Subsection III E, we focused on the decomposition path from the OH-perpendicular configuration to bidentate HCOO on the Cu(111) surface.

By using the vdW inclusive density functionals, the calculated desorption energy of the OH-perpendicular configuration is in very good agreement with the experiment, while the activation energy for the decomposition becomes lower than the desorption energy. Our results indicate that molecular decomposition takes place faster than desorption, seemingly in contradiction to the experimental report, in which no decomposition is observed when the Cu(111) surface is exposed to the gas phase HCOOH at room temperature.^19–23 To predict which process is faster than the other, we should not discuss based only on the activation energies, but we should take the entropy contribution to the reaction rate into account, as frequently discussed in the literature.^45,57 Therefore, in order to resolve the seeming contradiction between the room temperature experiment and calculated results, we performed the reaction rate analysis by assuming the first-order reaction process for HCOOH.

The desorption rate from the OH-perpendicular configuration is calculated by Eq. (8), while the decomposition rate is calculated by Eq. (10). We assumed that the translation and rotation modes are hindered for the transition state. For the OH-perpendicular configuration, the diffusion is also hindered, while the rotation is not. Therefore, we only consider rotation and vibrational partition functions for the adsorbed molecule. Figure 8 shows the calculated desorption and decomposition rates per surface adsorbate. In this analysis, we used the calculated results by the rev-vdW-DF2 energy functional. At low temperatures, the decomposition process becomes faster than the desorption process, while at higher temperatures, the desorption process becomes dominant due to large pre-exponential factor for the desorption process, and the transition temperature between the two processes is ∼187 K. These results are in very good agreement with the experimental situation that there is no HCOOH decomposition on Cu(111) by exposing the surface to the gas phase HCOOH at room temperature.

We found the calculated pre-exponential factors are 3 × 10¹⁶ s⁻¹ and 2 × 10¹² s⁻¹ for the desorption and decomposition, respectively. The calculated pre-exponential factor for the decomposition is four orders of magnitude smaller than that of the desorption. We obtain that the q_TS/q_ads term in decomposition process is ∼1, indicating that the molecule at the transition state is constrained at a specific site, leading to smaller pre-exponential factor compared to the desorption process. Our calculated pre-exponential factors are in reasonable agreement with the previously reported values by Sims et al.¹¹ Moreover, in the case of CO desorption from the Ni(111) surface, the pre-exponential factors are 10¹⁶ s⁻¹ to the 10¹⁷ s⁻¹, as reported by Ibach^44,58 and Elliot and Ward.⁵⁹ 

Note that our calculated desorption and decomposition rates per surface adsorbate for low temperature regime may contain error because we do not consider polymer formation, while in reality polymeric HCOOH is formed at low temperature. On the other hand, at room temperature, monomeric HCOOH is adsorbed and temporarily stays on the surface and it may decompose or be desorbed without forming polymeric HCOOH. Therefore, we do not intend to clarify the behavior of HCOOH at low temperature. However, we emphasize that our reaction rate analysis from first principles can clarify the behavior of HCOOH at room temperature very well.

We have investigated the adsorption and decomposition of HCOOH on Cu(111) by means of density functional theory (DFT) calculations. We show that by using the vdW inclusive density functionals, the adsorption energy is significantly improved over GGA. However, using vdW functionals, the activation energy for HCOOH decomposition becomes smaller than the adsorption energy, seemingly in contradiction to the experiment, in which no decomposition is observed when the Cu(111) surface is exposed to the gas phase HCOOH at room temperature. This contradiction is resolved by the reaction rate analysis with pre-exponential factors calculated from first-principles. The pre-exponential factor for the desorption process becomes much larger than that for the decomposition process, leading to faster desorption process at temperatures higher than ∼187 K, in very good agreement with the experiment. Our study demonstrates the importance of the kinetic theory along with accurate energetics for a more accurate and precise description of catalytic reactions.

See supplementary material for detailed data of the convergence test, diffusion and rotation processes, vibrational frequencies, and the effect of energy functionals on the desorption and decomposition rates.

We would like to thank to Professor Jun Yoshinobu (the University of Tokyo) for valuable discussions. S.E.M.P. is supported by Indonesia Endowment Fund for Education (LPDP) Scholarship from the Ministry of Finance of Indonesia. The present work was supported by the Advanced Catalytic Transformation Program for Carbon utilization (ACT-C) (Grant No. JP226MJCR12YU) of the Japan Science and Technology Agency (JST) and partly supported by Grants in Aid for Scientific Research on Innovative Areas “3D ActiveSite Science” (Grant Nos. 26105010 and 26105011) and “Hydrogenomics” (Grant No. JP18H05519) from the Japan Society for the Promotion of Science (JSPS), the Elements Strategy Initiative for Catalysts and Batteries (ESICB) supported by the Ministry of Education, Culture, Sports, Science, and Technology, Japan (MEXT). Numerical calculations were partially performed using Super Computers in the Institute for Solid State Physics (ISSP), the University of Tokyo.