FT²DP: large atomic model fine-tuned machine learning potential for accelerating atomistic simulation of iron-based Fischer-Tropsch synthesis

DFT calculations

All spin-polarized DFT calculations were carried out by using Atomic-orbital Based Ab-initio Computation at UStc (ABACUS) package ^{[44, 45]}. The SG15-optimized Norm-Conserving Vanderbilt (SG15-ONCV) multi-projector pseudopotentials ^{[46, 47]} were employed and the valence configurations were [H]1s$$ ^1 $$, [C]2s$$ ^2 $$2p$$ ^2 $$, [O]2s$$ ^2 $$2p$$ ^4 $$ and [Fe]3d$$ ^6 $$4s$$ ^2 $$. The generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) variant ^[48] was adopted for the exchange-correlation functional. The second generation of numerical atomic orbitals (NAOs) in the double-$$ \zeta $$ plus polarization function (DZP) form ^[49] was used as the basis set. The periodic boundary condition (PBC) and the $$ \Gamma $$-centered Monkhorst–Pack scheme^[50] for sampling the Brillouin zone were adopted in the DFT calculations, with an automated mesh determined by k-spacing = 0.14 Bohr$$ ^{-1} $$ and only one k-point for the direction without and with vacuum layers, respectively. The dipole correction perpendicular to the surface was applied for all DFT calculations of surfaces. The electron density criterion for electron self-consistency convergence was set at $$ 1 \times 10^{-7} $$, and the first-order Methfessel-Paxton (MP) smearing ^[51] was used for the occupation of orbitals. In geometry and transition state (TS) optimizations, the convergence criterion for the largest force among all atoms was set to 0.05 eV/Å.

DPA-2 and fine-tuning methodology

DPA-2 is a multi-task pre-trained LAM originating from the DP architecture and evolved from the DPA-1 model ^[52]. The DPA-1 descriptor has introduced an element-type embedding for encoding the elemental information covering the whole periodic table, and a gated self-attention mechanism ^[53] excelling in modeling the importance of neighboring atoms and re-weighting the interaction among them, which also makes the model generalizable and pre-trainable. Inheriting the DPA-1 backbone, the DPA-2 descriptor further enhances its resolution and generalizability of atomic representation through stacking multiple transformer ^[53] layers called representation transformer, incorporating operators such as convolution, symmetrization, localized self-attention, and gated self-attention, which can be interpreted as an E(3) equivariant graph neural network (GNN) and offers greater capacity compared to conventional GNNs ^[42], ensuring the robust capability of DPA-2 for serving as a LAM assembling comprehensive knowledge from massive pre-training data.

Besides having a more sophisticated model architecture, DPA-2 employs a multi-task training strategy for pre-training in multiple datasets labeled with different DFT settings to extract multidisciplinary knowledge. In particular, the multi-task DPA-2 model has multiple heads, and each head is an identical fitting network used to fit the DFT labels of each pre-training dataset from different downstream domains. During the pre-training process, the parameters within the DPA-2 descriptor are concurrently optimized through back-propagation using all pre-training datasets, while the parameters of the fitting network are updated exclusively with the specific pre-training dataset to which they are associated ^[42].

The pre-trained descriptor and fitting networks can be fine-tuned on specific downstream tasks, and the multidisciplinary knowledge learned from the multiple upstream datasets can help to reduce the consumption in model training and the amount of training data. In the fine-tuning process, the descriptor of the downstream model will be initialized with the pre-trained parameters, and the fitting network could also be initialized by choosing a fitting network in the pre-trained model. The energy bias in the fitting network will be aligned to the labels of the downstream dataset subsequently, and then the typical model training process will proceed with the initialized parameters incorporating upstream knowledge. Our FT$$ ^2 $$DP is constructed by this pre-training-fine-tuning methodology based on the publicly available pre-trained upstream DPA-2 model. All fine-tuning and model validation in this work were performed using DeePMD-kit software in version 3.0.0 beta 3 ^[30].

Double-to-single workflow for TS optimization

To investigate the catalytic reaction mechanism from a theoretical point of view, a crucial task is the optimization of TS structures of each elementary reaction for acquiring the activation energy. The climbing-image nudged elastic band (CI-NEB) ^{[54, 55]} and the dimer method ^{[56, 57]} are the two most popular TS optimization methods in heterogeneous catalysis simulation, representing two types of TS optimization methods respectively: double-ended methods that start from combining optimized initial state (IS) and final state (FS), and single-ended methods that are based on only one state, usually a guessed TS. It is well-known that the efficiency of single-ended methods highly relies on the quality of the input structure, but the optimization will converge quickly when the optimization reaches the quadratic region around TS. On the contrary, double-ended methods have no reliance on any guessed TS structure, but they tend to have convergence problems. In our work, we combine these two types of methods together as a workflow, named double-to-single (D2S), implemented in the atomic simulation environment (ASE) package ^[58] (as illustrated in Figure 1A) to make good use of them for accelerating the TS optimization process. The D2S workflow uses CI-NEB first to generate a rough reaction pathway with a relatively loose convergence criterion (usually 1.0 eV/Å for maximum of atomic forces), and then a single-ended method, such as the dimer, or a better choice, Sella algorithm based on iterative Hessian diagonalization and partitioned rational function optimization (P-RFO)^{[59, 60]}, is utilized by starting from the maximum point of the NEB pathway for strict TS optimization with target convergence criterion (usually 0.05 eV/Å for maximum of atomic forces). The free energy corrections, including zero-point energy (ZPE) and thermal corrections (translational, rotational and vibrational) for gas-phase molecules, as well as ZPE and vibrational contribution for adsorbates, are also computed using ASE. All these codes are open-sourced in the ATST-Tools suite ^[61], which supports using ABACUS and our FT$$ ^2 $$DP model as a property calculator. All structures in this part were visualized using ASE.

Figure 1. Scheme of model applications. (A) TS optimization, starting from a map of probable reaction network providing the reaction patterns of elementary reactions, and each reaction (taking CO dissociation as an example) is calculated by optimizing its IS, FS, and TS, where the TS is optimized by the D2S workflow; (B) Conventional and local schemes in optimizing the edge sites of the $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) slab using a GA. Orange: Fe atoms; Grey: C atoms; Red: O atoms. TS: Transition state; IS: initial state; FS: final state; D2S: double-to-single; GA: genetic algorithm.

Genetic algorithm for global optimization

We employed the genetic algorithm (GA) implemented in ASE ^[62], using FT$$ ^2 $$DP as the energy-evaluation calculator. In our global optimization process, the initial and subsequent populations each consist of 80 members, with each population exploring 80 new mutated candidates. If convergence is found difficult with this default setting, the population size and number of mutated candidates will be increased to 100. The mutation of candidates involves three operators implemented in the ASE-GA module: mirror (mirroring half atoms in a randomly oriented cutting plane), rattle (perturbing a part of atoms with a small random displacement), and permutation (switching the positions of a subset of atoms randomly), each with the same probability of 1/3. The stopping criterion in GA is regarded as met when the best candidate does not change for ten generations and the current generation number is at least 20. Following this, at least ten lowest-energy candidates are refined with PBE-DFT single-point calculations to identify the most stable structure. When using GA to obtain low-energy surface reconstructions, conventional implementation typically operates on an entire surface with constraints only on the z-axis^{[63, 64]}. This approach is unnecessary for predicting surface reconstructions that only involve local environments such as edge structures. In this work, we adopted a local setting implemented in ASE, specifying the involved atoms and where they are generated [Figure 1B]. All structures in this part were visualized using VESTA^[65].

Ab initio atomistic thermodynamics

To evaluate the stability of reconstructed surfaces with different compositions, the ab initio atomistic thermodynamics theory developed by Reuter and Scheffler^{[66, 67]} was used in this work. In this way, the surface energy ($$ \gamma $$) of a symmetric surface can be calculated from

where $$ G_\text{slab} $$ is the Gibbs free energy of a slab with two equivalent surfaces, $$ \mu_\text{Fe} $$ and $$ \mu_\text{C} $$ are the chemical potentials of Fe and C atoms, $$ N_\text{Fe} $$ and $$ N_\text{C} $$ are the numbers of Fe and C atoms, and $$ A $$ is the surface area. The Gibbs free energies were calculated by FT$$ ^2 $$DP under $$ T $$ = 523 K, which is a typical iron-based FTS temperature.

For a surface that is in equilibrium with a bulk with a fixed composition (e.g., FeC$$ _\text{x} $$), the chemical potentials of the contained elements are not all independent. In this work, $$ \mu_{\mathrm{Fe}} $$ and $$ \mu_{\mathrm{C}} $$ are related to the Gibbs free energy per formula unit of the FeC$$ _x $$ bulk ($$ \mu_{\text{FeC}_x} $$) as:

In our discussion, since not all structures contain the same amount of atoms, we defined a relative surface energy ($$ \Delta \gamma $$) as:

Here, the reference ($$ G_{\text{slab }} $$) is the clean slab without edge sites. The reason why the denominator "$$ 2A $$" in Equation (1) is replaced by "$$ A $$" in Equation (3) is that edge sites are only built on the top layer in this work, while the bottom layers remain fixed in structural relaxations. Additionally, $$ \Delta N_{\mathrm{Fe}} $$ or $$ \Delta N_{\mathrm{C}} $$ represents the differences in the number of Fe or C atoms between the reconstructed structure and the clean surface, respectively.

For convenience, we used the electronic energy of an isolated carbon atom ($$ E_\text{C} $$) as the reference for the carbon chemical potential, that is $$ \Delta \mu_\text{C} = \mu_\text{C} - E_\text{C} $$. Since the free energies and chemical potentials are relevant to temperature, pressure, and gas atmosphere. Here we used the results from Liu et al. to simulate a realistic iron-based FTS condition ($$ T $$ = 523 K, -6.60 eV $$ \leq $$$$ \Delta \mu _\text{C} $$$$ \leq $$ -7.45 eV)^[24].

FT²DP construction and validation

Our model, FT$$ ^2 $$DP, is constructed through fine-tuning on our downstream dataset from the upstream DPA-2.2.0 model ^[68], which is a pre-trained open LAM (OpenLAM) from the AIS Square website ^[69]. Thanks to the multi-task training protocol, this LAM was trained on more than 20 different datasets containing various physical and chemical systems including organic molecules, clusters, alloys, semiconductors, surfaces, and adsorbates through multi-task training mechanism, gathering multidisciplinary knowledge in one unified DPA-2 descriptor. Apart from the descriptor, the fitting model is a neural network containing three hidden layers with the typical numbers of neurons being (240, 240, 240) for all heads in the upstream DPA-2.2.0 model and our fine-tuned model.

There are 30, 656 frames in our FT$$ ^2 $$DP downstream dataset, including Fe-C-H-O element combinations and various types of structures, derived from the previous work by Liu et al. ^[25], and approximately 8, 000 structures were removed after the data cleaning procedures below to remove outliers and redundancies: (1) Removal of structures with identical DFT-calculated energy labels to eliminate redundant conformations; (2) Removal of the structures with fewer than 12 atoms per cell, which often represent isolated molecules or radicals in a cell. Such configurations are prone to DFT inaccuracies in PBCs or poor MLP generalizability; (3) Exclusion of structures having the absolute difference between model prediction and DFT results exceeds 80.0 meV/atom (energy) or 1.00 eV/Å (maximum atomic forces), where the model referenced here was fine-tuned from the upstream DPA-2 model on the original datasets following the same fine-tuning protocol detailed in the next paragraph. All DFT energies and forces were calculated by ABACUS following the computational settings mentioned above. A brief overview of this dataset is given in Table 1, showing the number of structures (with or without Fe) in different types, including three-dimensional (3D) bulks, two-dimensional (2D) surfaces, one-dimensional (1D) strings, and zero-dimensional (0D) clusters. Besides, a sketch-map visualization is shown in Figure 2 for illustrating the wide configuration distribution of our FT$$ ^2 $$DP dataset.

Figure 2. Sketch-map visualization of the FT$$ ^2 $$DP downstream dataset. Each point of this map represents an individual atomic configuration. The position of each point is determined by t-SNE of the learned descriptor of FT$$ ^2 $$DP model, and its color indicates the corresponding formation energy of the structure. FT$$ ^2 $$DP: Fine-tuned Fischer-Tropsch deep potential; t-SNE: t-distribution stochastic neighbor embedding

Our fine-tuning protocol was initialized by using the parameters of the global descriptor in pre-trained DPA-2.2.0 model and fitting network from the Domains_OC2M branch. The fine-tuning process on our dataset is done by following the default training process of the DPA-2 model with some setting modifications. In the default pre-training process of DPA-2.2.0 LAM, the learning rate starts from $$ 2 \times 10^{-4} $$ and gradually decreases to $$ 3.51 \times 10^{-8} $$ by an exponential decreasing scheme with each decrease performed at every 1/200 checkpoint of the total training step. The setting is usually effective for from-scratch training process, but the initial learning rate may be too large to be suitable for our fine-turning tasks owing to the loss of global knowledge in model descriptor after over-fitting in downstream dataset and possible gradient explosion in training practices. As a result, a relatively low initial learning rate $$ 2 \times 10^{-5} $$ was utilized together with one million training batches in our fine-tuning protocol.

To evaluate the generalizability of the FT$$ ^2 $$DP model, we performed a validation test where a model, named FT$$ ^2 $$DP-80p, was trained on a randomly selected 80% subset and tested on the remaining 20%. Validation results [Supplementary Figure 1, Table 2] present parity plots of formation energies and atomic forces, alongside R$$ ^2 $$ values, as well as the mean absolute error (MAE) and the root mean square error (RMSE) metrics comparing model predictions to DFT-calculated energies and forces. The FT$$ ^2 $$DP-80p model achieved comparable performance on training and validation subsets, confirming the robust generalizability. Furthermore, the accuracy of the final FT$$ ^2 $$DP model fine-tuned on the entire dataset is demonstrated in Figure 3A and B and Table 2 (final column), showing strong agreement between FT$$ ^2 $$DP prediction and DFT calculation for energies and forces. Moreover, the violin plots in Figure 3C and D illustrate the distribution of energy difference and atomic force difference between FT$$ ^2 $$DP prediction and DFT results, showing that although there are some outliers in a relatively wide range, our FT$$ ^2 $$DP model still achieves a good accuracy with energy deviation of less than 3 meV/Atom and a force deviation of less than 0.09 eV/Å for 75% of the FT$$ ^2 $$DP dataset. All these results indicate the model's promising reliability for direct usage in high-efficiency structural optimization and TS searching tasks. The accuracy of the final FT$$ ^2 $$DP model will be further demonstrated below by comparing FT$$ ^2 $$DP and DFT results for the lots of structures that emerged from atomistic modeling practices.

Figure 3. Validation results of the FT$$ ^2 $$DP model after fine-tuning on the entire dataset. Figures on the left are the parity plots comparing formation energies (A) and forces (B) from FT$$ ^2 $$DP against those from DFT on the dataset, with R$$ ^2 $$ equal to 1.000 and 0.9570, respectively. Figures on the right illustrate the violin plots of (C) distribution of DFT formation energies and atomic forces of the FT$$ ^2 $$DP dataset and (D) distribution of difference in energies and atomic forces between FT$$ ^2 $$DP predictions and DFT results on this dataset, and a box-plot without outlier is used to show the detailed distribution inside the peak region. FT$$ ^2 $$DP: Fine-tuned Fischer-Tropsch deep potential; DFT: density-functional theory.

Additionally, it is valuable to identify outlier structures with high prediction errors to investigate their characteristics and similarities. Structures in the entire FT$$ ^2 $$DP dataset with the top ten highest absolute prediction errors in energies ($$ \vert \Delta E \vert $$) and maximum absolute prediction errors in atomic forces ($$ \vert \Delta F \vert_{max} $$) are presented in Supplementary Tables 1 and 2 and Supplementary Figures 2 and 3, respectively. We observe that most of these structures exhibit disordered configurations or unphysical features. For example, 60% of these structures correspond to the same chemical formula C$$ _{6} $$H$$ _{12} $$O$$ _{6} $$, displaying disordered molecular configurations in a relatively small cubic cell with lattice vector equals to 10Å, rendering their physical state (gas/liquid/solid) indeterminate. Moreover, the third- and sixth-ranked structures in $$ \vert \Delta E \vert $$ contain unphysical H$$ _6 $$ clusters located in the vacuum layer of Fe$$ _{16} $$H$$ _{7} $$ surfaces. While these outliers reflect the equilibrium diversity of the dataset, which could improve model stability in the non-equilibrium region^[70], they also highlight the presence of unphysical configurations that necessitate data cleaning to remove the outliers for preventing training instability and downgraded model performance in practical atomistic modeling applications.

Reaction pathways

Many previous studies have demonstrated that certain surfaces exhibit much higher FTS activity, identifying them as the active surfaces^{[6, 71]}. A well-known example is the $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface, which has been shown to present relatively low CO dissociation and C-C coupling barriers^{[16, 24, 25]}. In this section, the FT$$ ^2 $$DP model and our D2S TS optimization workflow are utilized to investigate the key elementary reactions of FTS reaction pathways on the A-P5 site of $$ \chi $$-Fe$$ _5 $$C$$ _2 $$ (510) surface, known as an important active site for FTS process due to the participation of the lattice carbon and carbon vacancy on it in the surface reaction through a Mars-von Krevelen (MvK) mechanism revealed by previous studies ^{[6, 24, 25]}. In particular, we consider the dissociative adsorption of H$$ _2 $$ by Langmuir-Hinshelwood (L-H) mechanism, the dissociative adsorption of CO by the MvK mechanism, the competition between chain growth (C-C coupling) and carbon hydrogenation (C-H coupling) for carbon adsorbates, and the desorption of hydrocarbon compounds such as CH$$ _4 $$. The MvK reaction mechanism similar to previous works is also revealed in our investigation, especially for CO dissociation and chain growth process, illustrating the importance of the A-P5 site on iron carbide for FTS.

To further validate the accuracy of the FT$$ ^2 $$DP model in the TS optimization, we compared the results from three types of calculations: (1) purely FT$$ ^2 $$DP-based; (2) DFT single-point calculation after TS optimization by FT$$ ^2 $$DP, denoted as DFT@FT$$ ^2 $$DP; and (3) purely DFT-based, in Figure 4, illustrating that most reaction pathways calculated directly by FT$$ ^2 $$DP agree well with DFT benchmarks. However, a subset of cases exhibits notable energy discrepancies between FT$$ ^2 $$DP and DFT, particularly for TS. For example, the TS5 and TS6 states in Figure 4B show energy differences up to 0.30 eV. Such errors could propagate into significant uncertainties in identifying the rate-determining step (RDS) and microkinetics modeling of reaction networks. This highlights the inherent challenges in training MLPs to capture reactive dynamics and the limitations of relying solely on MLPs for reaction simulation, especially for out-of-distribution configurations. Further improvement can be attained in the DFT@FT$$ ^2 $$DP scheme. For instance, the energy differences of TS5 and TS6 states in Figure 4B are reduced to 0.01 and 0.03 eV, respectively, when comparing DFT@FT$$ ^2 $$DP results and pure DFT results. This hybrid approach aligns with established practices in MLP-based studies^[24], following the spirit of using low-fidelity (or rough) computational methods for geometry evaluation while performing high-fidelity (or precise) computational methods to refine the energy and electronic structure of crucial states^{[72, 73]}. While FT$$ ^2 $$DP can qualitatively map reaction networks without requiring posterior DFT energy corrections, quantitative accuracy cannot be fully guaranteed without further model refinement and dataset expansion to span the target reactive chemical space.

Figure 4. FTS reaction pathways on $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface, including (A) dissociation of H$$ _2 $$ and CO on A-P5 site; (B) *CH-*C coupling and following steps towards long-chain (C$$ _{2+} $$) hydrocarbon products; (C) *CH-*H coupling steps towards CH$$ _4 $$ product; and (D) *CH-*CH coupling steps towards CH$$ _2 $$CH$$ _2 $$ product. In these plots, red line denotes that the PES calculations are done merely by FT$$ ^2 $$DP, blue line denotes that the calculated data are acquired by single-point DFT calculation after FT$$ ^2 $$DP optimization, purple line denotes that the PES calculations are purely performed by DFT calculation, and C$$ _v $$ stands for carbon vacancy on A-P5 site. All the top views are for DFT-optimized structures. Orange: Fe atoms; Grey: C atoms; White: H atoms; Red: O atoms. FTS: Fischer-Tropsch synthesis; PES: potential energy surface; FT$$ ^2 $$DP: fine-tuned Fischer-Tropsch deep potential; DFT: density-functional theory.

Global optimization of reconstruction of edge sites

In addition to the $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface, some other surfaces have also been proposed as catalytically active in Fe-FTS process, including $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021), $$ \chi $$-Fe$$ _5 $$C$$ _2 $$($$ \overline{4} $$11), $$ \eta $$-Fe$$ _2 $$C(111), $$ \theta $$-Fe$$ _3 $$C(010), and $$ \theta $$-Fe$$ _3 $$C(031)^{[18, 25]}. In recent years, some studies have suggested that the FTS activity of low-coordinated atoms (typically at edges and corners) may differ significantly from that of regular surface atoms^{[12, 74]}, highlighting the importance of these sites in mechanistic research. In this section, four potential FeC$$ _x $$ active surfaces in iron-based Fischer-Tropsch, i.e., $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510), $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021), $$ \eta $$-Fe$$ _2 $$C(111), and $$ \theta $$-Fe$$ _3 $$C(010) [Supplementary Figure 4], were considered to investigate the surface morphology and stability of edge sites. Each surface was modeled using a slab with a thickness of more than 10 Å, with the bottom two or three layers of FeC$$ _x $$ fixed during structural relaxations. A vacuum layer of 20 Å was introduced to separate two neighboring slabs for all surfaces. To model the edge sites, a 'half-surface' scheme was employed: half atoms of the top layer were removed, maximizing the distance and simultaneously minimizing the interaction between two adjacent edge sites. This approach is particularly useful for GA calculations, as it allows for two independent local searches on two edges to identify the most stable reconstructed structure.

Before performing GA calculations, we first checked the accuracy of FT$$ ^2 $$DP on large clean surfaces and unreconstructed edge structures. We used the FT$$ ^2 $$DP model to optimize the structures and then calculated their single-point energies on the PBE level. In Supplementary Tables 3 and 4, we report the energy differences between PBE single-point energies for FT$$ ^2 $$DP and PBE-optimized structures. The energy differences for all tested structures are less than 1.2 meV/atom, indicating that FT$$ ^2 $$DP is sufficiently accurate for investigating surface reconstructions of iron carbides.

There are two possible terminations of $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) ^[75]: Fe- and C-terminated surfaces, and their relative stability is dependent on the chemical potential of carbon ($$ \Delta \mu _\text{C} $$). In modeling the edge sites, we chose to remove four consecutive rows of Fe atoms (24 Fe atoms in total), along with the corresponding C atoms. These unreconstructed structures with edge sites may contain many under-coordinated Fe and C atoms, making global optimization essential to obtain low-energy reconstructed structures. A representative low-energy reconstructed structure of $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface with edge sites is shown in Figure 5A. In this configuration, only the relocation of carbon atoms is observed, while the positions of iron atoms remain largely unchanged, which is consistent with previous findings ^[24]. In contrast, reconstructions involving the migration of both iron and carbon atoms on Fe-terminated surfaces were identified, together with some newly formed [Fe$$ _4 $$C] squares [Supplementary Figure 5]. Since the atom number of different structures is not identical, we used Equation (3) to characterize their stabilities and plotted the relative surface energies ($$ \Delta \gamma $$) against $$ \Delta \mu _\text{C} $$ in Figure 5B. Notably, reconstructed C-terminated surfaces, such as C-4-7 in Figure 5B, exhibit higher stability than other clean surfaces at the upper limit of $$ \Delta \mu _\text{C} $$, including the reconstructed $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface, which is built according to the previous work of Liu et al.^[24]. Previous studies ^[25] have demonstrated the important role of [Fe$$ _4 $$C] squares (or A-P5 sites) in facilitating both CO dissociation and C-C coupling reactions. We obtained a reconstructed structure containing a row of inclined [Fe$$ _4 $$C] squares as shown in Figure 5A, which implies that the reconstructed edge sites on $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surfaces could serve as active sites for iron-based Fischer-Tropsch reactions.

Figure 5. Structure of the lowest-energy reconstruction of each surface and relative surface energies of low-energy reconstructions with respect to $$ \Delta \mu _\text{C} $$. (A) and (B) $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface, (C) and (D) $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021) surface, (E) and (F) $$ \theta $$-Fe$$ _3 $$C(010) surface. When referring to $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surfaces, "Fe-" or "C-" denotes the Fe-terminated or C-terminated surface, respectively. The numbers in the names of surfaces [e.g., 6-15 in $$ \theta $$-Fe$$ _3 $$C(010) surface] represent the range of column numbers (counted from left to right) from which iron atoms have been removed. The cyan curve in (B) marked by "clean-rec" corresponds to the reconstructed $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(510) surface reported by the previous work of Liu et al.^[24].

$$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021) and $$ \theta $$-Fe$$ _3 $$C(010) were also considered as potential active surfaces in iron-based FTS ^{[17, 25]}. Among them, the $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021) surface exhibits greater structural complexity compared to the other three surfaces, consisting of both [Fe$$ _4 $$C] squares and [Fe$$ _5 $$C] pentagons [Supplementary Figure 4]. Despite this complexity, our workflow successfully performed global optimizations and identified the reconstructed surfaces with edge sites. The most stable reconstructed structure is illustrated in Figure 5C, where newly formed [Fe$$ _4 $$C] squares are also observed (also see Supplementary Figure 6 for additional structures). The relative surface energies ($$ \Delta \gamma $$) against $$ \Delta \mu _\text{C} $$ are plotted in Figure 5D. Although all of them are less stable than the clean $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021) surface, with relatively small energy differences, they may still play certain roles in the FTS at a high temperature.

Figure 6. Results of variable-composition calculations on the $$ \eta $$-Fe$$ _2 $$C(111) surface. (A) A schematic representation of the variable-composition search, where "-3C" indicates three fewer carbon atoms compared to the unreconstructed structure; (B) Side and top views of unreconstructed and reconstructed structures of the most stable configuration at the lower $$ \Delta \mu _\text{C} $$ limit; (C) Plot of the relative surface energies of the most stable reconstructions with different total carbon numbers against $$ \Delta \mu _\text{C} $$.

The clean $$ \theta $$-Fe$$ _3 $$C(010) surface exhibits a highly ordered structure. The most stable reconstructed structure of $$ \theta $$-Fe$$ _3 $$C(010) surface with edge sites is presented in Figure 5E. This behavior resembles the phenomenon of "island decay" in surface science^[76], where a rough surface evolves toward a smoother surface to achieve higher stability (also see Supplementary Figure 7 for two additional reconstructions). Similar to the $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021) surface, newly formed [Fe$$ _4 $$C] are observed in these reconstructed $$ \theta $$-Fe$$ _3 $$C(010) surfaces, while the original surface morphology is still preserved. The relative surface energies ($$ \Delta \gamma $$) as a function of $$ \Delta \mu _\text{C} $$ are also plotted in Figure 5F. Notably, none of the reconstructed $$ \theta $$-Fe$$ _3 $$C(010) surfaces are found to be more stable than the clean surface. Moreover, the energy differences are significantly larger than those on $$ \chi $$-Fe$$ _5 $$C$$ _2 $$(021) surfaces, indicating that the formation of edge sites on $$ \theta $$-Fe$$ _3 $$C(010) is much less favorable under a typical Fischer-Tropsch reaction condition.

Compared to other surfaces, the $$ \eta $$-Fe$$ _2 $$C(111) surface exhibits a relatively simple structure and all surface Fe-4fold square sites are fully occupied by carbon atoms [Supplementary Figure 4]. Here, we demonstrated variable-composition global optimization calculations involving different numbers of carbon atoms. When performing fixed-composition calculations, two edges of each surface were optimized using GA separately. While performing variable-composition calculations, the number of atoms at each edge can be changed and these searches will be repeated several times under different total numbers of atoms [Figure 6A]. One example of a reconstructed surface is shown in Figure 6B, where carbon atoms migrate from one edge to the other, forming C-C dimers to fully bond with neighboring atoms. The relative surface energies ($$ \Delta \gamma $$) of the reconstructed surfaces against $$ \Delta \mu _\text{C} $$ are plotted in Figure 6C. Structures with varying carbon numbers are detailed in Supplementary Figure 8. Although all of the reconstructed surfaces are less stable than the clean $$ \eta $$-Fe$$ _2 $$C(111) surface, the energy differences between them are relatively small, especially for carbon-rich structures under a high $$ \Delta \mu _\text{C} $$. For example, at the lower limit of $$ \Delta \mu _\text{C} $$ = -7.45 eV, the most stable reconstructed surface is identified as "-3C". This means the surface energy is minimized when three fewer carbon atoms are present compared to the unreconstructed structure. Conversely, at the upper limit of $$ \Delta \mu_\text{C} $$ = -6.60 eV, the relative surface energy decreases monotonically with an increasing carbon number. In the end, the stability of the surface with most carbon atoms ("+C") is very close to the clean surface. This suggests carbon deposition on the iron-carbon surface is presumed to be thermodynamically favorable under a high $$ \Delta \mu _\text{C} $$, typically near the reactant-catalyst equilibrium in a carbon-rich gas environment.

Acknowledgments

The authors acknowledge Dr. Duo Zhang from AI for Science Institute for helping with model training and dataset visualization and Dr. Yike Huang from AI for Science Institute for sharing template code on utilizing ASE-GA in surface systems. The authors also appreciate the support of the High-performance Computing Platform of Peking University and National Supercomputer Center in Tianjin for the computational resources.

Authors' contributions

Contributed equally to this work: Liu, Z.Q.; Deng, Z.

Dataset preparation, model training and evaluation, data analysis, writing: Liu, Z. Q.; Deng, Z.

Discussion of results, revision: Zhao, H.; Wang, H.; Chen, M.

Project conceptualization, methodology, supervision, revision, funding acquisition: Jiang, H.

Availability of data and materials

The dataset and model used in this study are both available on AIS Square (https://www.aissquare.com/). In detail, the dataset with the structural information, DFT label information, and DFT input setting files can be accessed at https://www.aissquare.com/datasets/detail?pageType=datasets&name=FT2DP-dataset-FeCHO-v1&id=306, and the model with its input scripts is available at https://www.aissquare.com/models/detail?pageType=models&id=307. Additionally, the codes for TS optimization are all available in the ATST-Tools repository (https://github.com/QuantumMisaka/ATST-Tools). All other data and codes supporting the findings presented in this work are available from the corresponding author upon reasonable request.

Financial support and sponsorship

This work is financially supported by the National Key Research and Development Program of China (Project no. 2022YFB4101401) and National Natural Science Foundation of China (Project no. 22273002).

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

® The Author(s) 2025.

Supplementary Materials