Three Dimensional Activity Volcano Plot under External Electric Field

Much like Earth’s climate and spin glass, heterogeneous catalysis is an archetypal complex system and plays a pivotal role in our society. Electric field (EF) can manipulate the chemical reactivity and selectivity by tuning the relative stability of polar or ionic reactants and transition states in chemical reactions [1, 2, 3, 4, 5, 6], opening up exciting opportunities to control catalytic reactions in an on-demand manner using EFs, hereinafter referred to as electrostatic catalysis. To take full advantage of electrostatic catalysis, a fundamental challenge is to predict the response of heterogeneous catalysis, a complex physicochemical system often involving a network of chemical reactions and physical processes, to collective reaction parameters, specifically field strength and temperature.

The presence of an EF under realistic reaction conditions can have both intrinsic and extrinsic origins [7, 2, 4, 6]. It is more convenient in practice to directly apply an external EF (EEF) to a reactor as the strength and direction of the field is readily tunable. To apply EEFs in large-scale reactions, several reactor systems have been developed [8, 9, 10], i.e., probe-bed-probe reactors that have a conductive catalyst bed placed in the gap between two external probes and continuous-circuit reactors in which the conductive catalyst bed is integrated as an electrical circuit. The catalyst bed often consists of catalytically active metal nanoparticles that can generate very high surface EFs in the presence of EEFs due to the curvature-induced field enhancement effect [11]. Additionally, surface EFs are always perpendicular to the metal surface, potentially allowing for a fine control of the relative orientation between the EEF and adsorbed molecules. Therefore, understanding the effects of EEFs on heterogeneous catalysis by metals is a major focus of electrostatic catalysis [12, 13, 11].

Methanol, as an important hydrogen-storage liquid fuel [14], can be used for on-board hydrogen production through steam reforming of methanol (SRM) [15], in which methanol reacts with steam and produces hydrogen, CO, and CO₂ at about 300 ^∘C. However, the extensively used catalysts in industry often suffer from catalyst deactivation due to coke deposition and sintering [16]. Since the elementary reactions of SRM involve adsorbates possessing large electric dipole moments, we expect that an EEF will bring appreciable impact on the kinetics and thermodynamics of the SRM process as the energetics of molecules of large dipole moments are susceptible to EEFs [17, 5]. Considering the critical role of SRM played in the green methanol economy [14, 18], it is imperative to design next-generation catalysts superior to current ones, particularly for low-temperature SRM. In this work, we propose to use EEF as a “smart agent” to improve both the catalytic activity and coke resistance.

First-principles-based computational catalyst design is emerging as a promising approach to obtain high-quality catalysts [19, 20, 21]. Pertinent to electrostatic catalysis, finite-field density functional theory (DFT) calculations have been used to quantify the activation energies ($\Delta E_{a}$) of elementary reactions on metal surfaces in the presence of EEFs of varying magnitudes, revealing important atomistic insights that helped the understanding of EEF effects [3, 13, 22]. However, most previous studies only focused on a few representative electric fields (e. g., negative versus positive) and their impact on the catalytic activity of a specific metal (e. g., Ni, Ag) [13, 22], partly due to the expensive computational cost of finite-field DFT calculations. For a complex process like SRM comprising a network of elementary reactions, the EEF dependence of the overall catalytic performance is likely nonlinear such that the search of an optimal EEF should be performed in order to maximize the efficiency but has never been done. Moreover, the design of high-quality electrostatic catalysts targeting specific operation conditions requires detailed understandings of the collective effects of EEFs and various reaction parameters such as temperature, gas composition, and partial pressure. All these difficulties essentially forbid a DFT-based high-throughput computational catalyst screening in the presence of EEFs.

Here, using SRM as an example, we develop a deep-learning-assisted first-principles-based multiscale method that leads to the first three-dimensional (3D) activity volcano plot under EEFs, quantifying the field-dependent SRM activity for thousands of metallic alloys. The multsicale framework consists of three steps, as shown in Fig. 1a. First, we derive a general analytical relationship between the electric field ($\mathbf{F}$) and $\Delta E_{a}$ of a surface reaction within the harmonic approximation. This enables a rapid and accurate predication of $\Delta E_{a}[\mathbf{F}]$ at any given field strength. The analytical relationship is then incorporated into a microkinetic mode of SRM on Ni(111) surface, resulting in a continuous stirred-tank reactor model to simulate SRM under EEFs. The microkinetic modeling reveals a highly nonlinear temperature-dependent EEF effect: a positive EEF can increase the conversion of CH₃OH at high temperatures ($>$ 350^∘C) but suppress the conversion at low temperatures ($<$ 250^∘C), highlighting the necessity of multiscale modeling for catalyst design under EEFs. Finally, by combining the linear scaling relationship [23, 24], a DFT-based deep neural network (DNN) that rapidly predicts C* and O* adsorption energies [19], and a simplified kinetic model derived from the microkinetic model, we are able to construct the 3D activity volcano plot to identify high-quality catalysts for low-temperature SRM under EEFs.

MethodDFT calculations are performed using Vienna ab initio Simulation Package (VASP) [25, 26]. The interaction between core ion and electrons is described by the projector augmented wave (PAW) method [27]. The Perdew-Wang-91 (PW91) functional is chosen as the exchange-correlation functional [28]. The plane-wave kinetic energy cutoff is 400 eV and a 5$\times$5$\times$1 Monkhorst-Pack grid is used to sample the Brillouin zone. A three-layer 3$\times$3 Ni slab with a vacuum region of more than 10 Å thick is constructed to simulate the Ni(111) surface. The method proposed by Neugebauer and Scheffler [29] is employed to apply an EEF normal to the metal surface. Climbing image nudged elastic band (CL-NEB) method [30] and dimer method [31] are used to identify transition states of surface reactions, and the electronic energy and atomic force are converged respectively to 10^-7 eV and 0.01 eV/Å. More computational details can be found in Supplementary Materials.

Activation and adsorption energy under EEFs within the harmonic approximationFor an elementary surface reaction step, the EEF-induced changes in adsorption and activation energy have been attributed to various factors such as modified charge transfer between adsorbates and metal surfaces [32], bond elongation/contraction of adsorbates as well as the shift of metal work function [33]. These important yet intimately coupled atomistic mechanisms seem to suggest a quantitative determination of the EEF-induced energy change can only be achieved via genuine ab initio modeling that captures all the factors. We now show that $\Delta E_{a}[\mathbf{F}]$ at a given field ($\mathbf{F}$) can be readily obtained using only zero-field parameters within the harmonic approximation.

Without loss of generality, a simple elementary reaction, R*$\rightarrow$TS*$\rightarrow$P* is considered, where R*, TS*, and P* are reactant, transition state, and product adsorbed on the metal surface, respectively, and the EEF is applied along the surface normal. For an isolated species $i$ under an EEF, its free energy in the form of Taylor series up to the second order is

where $E^{0}$ is the zero-field energy, $\mathbf{d}$ is the dipole moment, and $\alpha$ is the polarizability. Equation 1 applies to both isolated molecules as well as the slab model with adsorbates (denoted as $s$+$i$). It is straightforward to derive

where $\Delta E^{0}_{a}=E^{0}_{s+\rm TS^{*}}-E^{0}_{s+\rm R^{*}}$ is the zero-field activation energy. We note that $\mathbf{d}_{s+\rm TS^{*}}$ and $\alpha_{s+\rm TS^{*}}$ ($\mathbf{d}_{s+\rm R^{*}}$ and $\alpha_{s+\rm R^{*}}$) are the electric dipole moment and polarizability, respectively, of the whole bounded system comprised of the adsorbed TS* (R*) and the slab used to model the surface. It is evident that all reaction-specific parameters ($\Delta E_{a}^{0}$, $\mathbf{d}$, $\alpha$ ) in Eq. 2 can be computed from DFT and perturbation theory at zero field. A generalization of Eq. 2 gives the adsorption energy under EEFs,

where $E_{\rm ad}^{0}=E^{0}_{s+i}-E^{0}_{s}-E^{0}_{i}$. We note that our definition of adsorption energy captures the effect of EEF-driven gas diffusion and is subtlety different from the conventional definition (see details in Supplementary Materials, Sect. II).

Expressions similar to Eqs. 2–3 have been derived previously[3, 22] but, surprisingly, have never been used to calculate $\Delta E_{a}[\mathbf{F}]$ or $E_{\rm ad}[\mathbf{F}]$, likely because of the general assumption that EEF-induced structural and electronic changes of metal surface-adsorbate complex would be highly anharmonic thus beyond the applicability of harmonic approximation. We compare the EEF-induced energy change, $\Delta\Delta E_{a}[\mathbf{F}]=\Delta E_{a}[\mathbf{F}]-\Delta E_{a}^{0}$ and $\Delta E_{\rm ad}[\mathbf{F}]=E_{\rm ad}[\mathbf{F}]-E_{\rm ad}^{0}$, obtained with Eqs. 2–3 using only zero-field parameters (see values in Supplementary Materials, Sect. III) and those from direct DFT calculations in Fig.1b. With a mean absolute error of only 1 meV and 2 meV for $\Delta\Delta E_{a}[\mathbf{F}]$ and $\Delta E_{\rm ad}[\mathbf{F}]$, respectively, the simple analytical relationships described by Eqs. 2–3 have remarkable accuracy over a wide range of field strengths. We further compare the analytical and DFT values of $\Delta E_{\rm ad}[\mathbf{F}]$ for three representative molecules, H₂O, CH₃OH, and CH₃O in Fig.1c, and find a satisfying agreement for $-0.4<\mathbf{F}<0.7$ eV/Å. Consistent with previous DFT investigations [3, 34], a negative field will induce desorption ($\Delta E_{\rm ad}[\mathbf{F}]>0$) while a positive field may facilitate the adsorption ($\Delta E_{\rm ad}[\mathbf{F}]<0$) for a molecule with negative electric dipole moment on the metal surface. It is only when the molecules become desorbed due to a strong negative field (e. g., $-0.5$ V/Å for H₂O), the analytical value of $\Delta E_{\rm ad}[\mathbf{F}]$ starts deviating from the DFT result.

The demonstrated quantitative nature of a rather simple analytical theory linking field strength to EEF-induced energy change seems counterintuitive. From Eqs. 2–3, it is clear that only the whole bounded system (the slab with an adsorbed molecule) is relevant to EEF-dependent terms. We argue that for the the whole bounded system, the impact of an EEF on the total energy is a small perturbation that could be well described within the harmonic approximation, despite highly anharmonic local changes between the adsorbate and the metal surface (see numerical verification in Fig. S3).

Microkinetics of SRM under EEFsThe reaction network of SRM, illustrated in Fig.2a, includes CH₃OH and H₂O dehydrogenation and CO* oxidation. Based on Eqs. 2–3, we construct a microkinetic model that simulates the Ni-catalyzed SRM process in the presence of EEFs under realistic conditions that take into account the impacts of gas composition, temperature, and inlet velocity. The validity of the microkinetic model is confirmed by comparing theoretical predictions with available experimental results (see details in Supplementary Materials, Sect. IV). For a 3:1 steam to methanol ratio (S:M = 3:1), the temperature- and field-dependence of the conversion of CH₃OH is presented in Fig. 2b where we define an “optimal field line” (OFL) that tracks the electric field resulting in the highest conversion. The OFL reveals several interesting characteristics of SRM under EEFs. At high temperatures $\geq 310^{\circ}$C, a positive field can promote the overall reaction rate compared to the zero-field case, and the degree of enhancement scales with the field strength. Interestingly, at an intermediate temperature such as $300^{\circ}$C, a large positive field that is to the right of the OFL will suppress the conversion. When the OFL crosses the zero field and enters into the negative field region, it means a negative field now can facilitate the conversion; this happens at $\approx$250 ^∘C, but a too strong negative field also becomes unfavorable.

We perform sensitivity analysis by computing the partial derivatives of methanol conversion rate with respect to the rate constant of the elementary reaction [35]. As shown in Fig. 2c, there are two RDSs in SRM, the dehydrogenation of CH₃O*, CH₃O*+* $\rightarrow$ CH₂O*+H*, and the consumption of CO*, CO*+O*$\rightarrow$CO₂. The nontrivial evolution of OFL is due to a delicate balance between the two RDSs with the dehydrogenation of CH₃O* being more EEF sensitive while the consumption of CO* being more temperature sensitive. Specifically, as shown in Fig. 2d, at a high temperature of $350^{\circ}$C, the overall reaction rate of SRM is dictated by the dehydrogenation of CH₃O* as the consumption of CO* is fast enough at such high temperature. Therefore, a stronger positive EEF that reduces the activation energy of CH₃O* dehydration leads to increased methanol conversion. When the temperature is reduced to $300^{\circ}$C, the consumption of CO* remains fast such that applying a positive EEF below $\mathbf{F}_{\rm OFL}$ can still facilitate the conversion by increasing the rate of CH₃O* dehydrogenation. However, after passing $\mathbf{F}_{\rm OFL}$, the CO* generation will outspeed the consumption, causing CO poisoning of Ni. At a low temperature of $250^{\circ}$C, the CO* consumption becomes rather slow thus requiring a negative field to suppress CH₃O* dehydration to inhibit CO poisoning; a too strong negative field that severely slows down CH₃O* dehydration will unsurprisingly cause a low methanol conversion rate. Simply put, the OFL essentially marks the boundary between CH₃OH-dehydrogenation-controlled region and CO-poisoning-controlled region.

It is well known that coking is the most common deactivation mechanism of industrial SRM catalysts [36, 37], and increasing the steam pressure is the general approach to suppress coking. Though a high S:M ratio (e. g., 3:1) is beneficial for enhancing coke resistance, it will raise energy consumption for steam heating, reduce the methanol conversion rate, and decrease the H₂ content of final products. It would be helpful to have an efficient tool to identify the optimal conditions for balanced SRM rate and coke resistance. For S:M = 1:1, we construct a contour plot of methanol conversion rate as a function of temperature and EEF strength (Fig. 2e), showing both OFL and constant carbon concentration lines (3CLs). For a target level of coke resistance (a selected 3CL, denoted as 3CL_t), the lower part of the merged OFL and 3CL_t gives the best operating EEF for a wide range of temperatures.

3D activity volcano plot under EEFsA detailed study of the microkinetic model of SRM shows that the methanol conversion rate ($r_{\rm M}[\mathbf{F}]$) can be quantified with a simplified pathway [38],

in which CH₃O*+* $\rightarrow$ CH₂O*+H* is the RDS. The corresponding rate equation reads

Here, $P_{\rm CH_{3}OH}[\mathbf{F}]$, $P_{\rm H_{2}}[\mathbf{F}]$, $\theta[\mathbf{F}]$ are the field-dependent partial pressure of CH₃OH, H₂, and the coverage of active sites that depend on the adsorption energies of CO*, O*, and H* (see Section VI in Supplementary Materials); $A$ is a temperature dependent parameter deduced from DFT calculations. The simplified kinetic model based on Eq. 4 is capable of reproducing the results of the microkinetic model based on the complete reaction network but with a fraction of the computational cost.

To construct the activity volcano plot under EEFs for a large number of metallic catalysts, accurate and rapid determinations of $E_{\rm ad}[\mathbf{F}]$ and $\Delta E_{a}[\mathbf{F}]$ are prerequisite. We find that the EEF-induced energy change turns out to be rather insensitive to the metal type, which greatly simplifies the problem (see Supplementary Material Sect. VII). That is, for an adsorbate $i$ on a given metal surface ($m$), the field-dependent adsorption energy is reduced to $E_{{\rm ad}}[\mathbf{F},m]=E_{{\rm ad}}^{0}[m]-(\mathbf{d}_{s_{0}+i}-\mathbf{d}_{s_{0}})\cdot\mathbf{F}-\frac{1}{2}(\alpha_{s_{0}+i}-\alpha_{s_{0}})\mathbf{F}^{2}$, where $E_{{\rm ad}}^{0}[m]$ is the zero-field adsorption energy on metal $m$ that can be estimated using C* and O* adsorption energies on metal $m$ with the well-known linear scaling relationship [19] while the last two terms only depend on $\mathbf{F}$ and are computed using a reference metal surface $s_{0}$. Equipped with a DNN that can quickly predict the zero-field adsorption energies of C* and O* [19] for various metallic alloys and the simplified kinetic model (Eq. 4), we finally construct a 3D activity volcano plot that quantitatively predicts the menthol conversion in the presence of EEFs for 1711 metallic alloys (Fig. 3a) at a low temperature of 200 ^∘C. This 3D volcano plot contains rich information and is worthy of detailed investigations. Here we only highlight the catalytic performance of NiZn. Despite having a low SRM activity at zero field, NiZn can promote the methanol conversion in the presence of an EEF of 0.5 V/Å at 200 ^∘C. Combined with its proved coke resistance [19], EEF-assisted NiZn is recommended as an efficient and environment-friendly catalyst for SRM at 200 ^∘C and S:M = 1:1.

Conclusion.In summary, we demonstrate that a simple analytical theory within the harmonic approximation can accurately predict finite-field energetics using only zero-field parameters, eliminating the needs of expensive finite-field DFT calculations for field-dependent microkinetic modeling. Focusing on an important process in green methanol economy, steam reforming of methanol, we reveal nontrivial collective effects of EEF and temperature: a positive EEF enhances the methanol conversion at higher temperatures while suppressing the overall reaction at lower temperatures. The introduction of OFL and 3CL_t in the temperature–EEF parameter space allows for facile determinations of best operating EEF for a wide range of temperatures. Finally, using a chemistry-based simplified kinetic model and a first-principles-based DNN, we construct a 3D activity volcano plot under EEFs for 1711 metallic alloys and identify EEF-assisted NiZn for low-temperature steam reforming of methanol. We expect the multiscale approach developed in this work can be readily applied to other heterogeneous catalysis by metal under EEFs.