Trade-off between O₂ activation and active-site regeneration on biaxially strained Co-doped MoS₂ monolayers: a density functional theory study

Abstract

Strain engineering can tune O₂ adsorption, activation, and dissociation on two-dimensional transition-metal dichalcogenide catalysts; however, its synergistic impact on O₂ activation/dissociation and active-site regeneration, both of which are required for sustained turnover, remains unclear. Herein, spin-polarized density functional theory is used to examine O₂ activation and regeneration at a substitutional Co site in monolayer MoS₂ (Co@V_S). Phonon calculation results obtained for V_S and Co@V_S monolayers show the absence of imaginary modes at a 10% biaxial tensile strain (the highest strain examined), confirming dynamical stability at the upper bound of the studied strain window. With increasing strain, the O₂-adsorbed Co@V_S site undergoes a crossover between the S-preserved and S-reconstructed configurations near 4.5%, promoting the activation process. Consequently, O₂ dissociation preferentially follows a Mo-assisted pathway, yielding a deeply stabilized Co–O–Mo termination. The regeneration process assessed using an atomistic oxygen-migration proxy is increasingly hindered by strain: between 4.5% and 5.0%, oxygen penetrates the lattice deeply and disrupts the site, whereas at a strain of 5.5%, oxygen removal is rate-limited by a large lateral diffusion barrier of 1.63 eV, consistent with the strengthened Co–O interactions. Overall, the applied strain facilitates O₂ dissociation, but it can also deepen oxide-like product wells and kinetically impede oxygen removal along the lattice-site hopping pathway. More broadly, the obtained results highlight an activation–regeneration trade-off that may be overlooked when strain engineering is primarily evaluated by activation descriptors, suggesting that optimal strain windows for strain-tuned single-atom catalysts should balance dissociation facilitation against the site recovery feasibility.

---

1. Introduction

Two-dimensional transition-metal dichalcogenides have emerged as versatile and tunable platforms for oxygen electrocatalysis.^1–3 They are widely used in air cathodes for alkaline Zn–air batteries,^4,5 oxygen reduction reaction (ORR) electrodes operating under both acidic and alkaline conditions,^6–8 and in the on-site synthesis of hydrogen peroxide via the two-electron ORR pathway (2e⁻-ORR).⁹ Recent advances have further demonstrated the potential of engineered two-dimensional catalysts for ORR through heterostructure design and defect engineering.^10–12 Among these materials, MoS₂ has been extensively studied, and ORR-related oxygen electrocatalytic activity was previously observed for MoSe₂, WS₂, and WSe₂.^4–9

In MoS₂, the electronically inert basal plane contrasts with chemically active point defects, where substitutional modifications can create metal–oxygen binding sites capable of activating molecular oxygen (O₂).^13–15 These O₂-active sites in MoS₂ are generally categorized into edge sites and interior sulfur vacancies (V_S). In large-area films, sulfur vacancies are formed during chemical vapor deposition or liquid-phase exfoliation with areal densities ranging from 10¹⁰ to 10¹³ cm⁻², as confirmed by aberration-corrected scanning transmission electron microscopy combined with electron energy loss spectroscopy.^16,17 In small-area single crystals, edge sites typically predominate, whereas in wafer-scale films, V_S often serve as the primary active sites.^16,17 Theoretically, O₂ dissociation on the pristine MoS₂ basal plane requires a high energy barrier (∼1.6 eV); however, the presence of a V_S reduces this barrier by half to approximately 0.8 eV, making the dissociation process more feasible.¹⁸ Meanwhile, sustained O₂ turnover necessitates the continuous regeneration of clean vacancies, because a vacancy that has reacted with oxygen becomes difficult to reactivate. Density functional theory (DFT) and spectroscopic studies have demonstrated that an oxygen atom occupying a sulfur vacancy in monolayer MoS₂ is strongly stabilized, with adsorption energies around −4 eV and repair energies of approximately 2–4 eV relative to those of bare V_S. The resulting oxygen-filled defects eliminate vacancy-derived mid-gap states, restore a nearly pristine band structure, and yield long-lived photoluminescence recovery, all of which are consistent with a thermodynamically deep, electronically passivated state that is unlikely to participate in further O₂ activation cycles.^19–21

Transition-metal substitution—particularly Co substitution at a sulfur vacancy (Co@V_S)—has been proposed as an effective strategy for overcoming this limitation.^{14,15,22–24} In these Co-decorated vacancies, the Co atom forms a strongly bound single-atom site and activates the adsorbed O₂ species through the metal–oxygen hybridization and partially occupied Co d states near the Fermi level. This interaction promotes π* back-donation and weakens the O–O bond.^14,15,22,23 Among the 3d dopants, Co is especially attractive owing to its favorable ionic radius and electronic configuration, which ensure structural stability within the MoS₂ lattice. Additionally, its partially filled d orbitals provide a balanced oxygen-binding strength and experimentally verified catalytic activity in oxygen-involving reactions.^25,26 Besides substitutional doping, the further tuning of the local electronic environment near V_S offers a promising approach to controlling the O₂ reactivity. Previous studies on Co, Ni, Fe, and Cu dopants have demonstrated strong anchoring of transition metal atoms at V_S and significant reorganization of the local electronic structure.^{13–15,22,23} Notably, Co@V_S stands out because it (i) forms a thermodynamically stable single-atom site at the vacancy, (ii) hosts d states near the Fermi level that efficiently couple to the O₂ π* manifold, and (iii) weakens the O–O bond without causing excessive oxygen over-binding,^{13–15,22,23} making it an ideal platform for strain-tunable oxygen activation.

Biaxial tensile strain, which can reach approximately 5–6% in suspended monolayer MoS₂ membranes and 1% in polymer-supported devices, further modulates band edges and orbital hybridization.^27–30 Strain engineering can induce direct-to-indirect bandgap transitions and alter the Mo-d/S-p hybridization patterns near both the valence band maximum and conduction band minimum.^31,32 This pronounced strain sensitivity of the MoS₂ lattice suggests that in our Co@V_S–O₂ system, biaxial strain can modulate the three-way orbital coupling among the MoS₂ lattice, Co dopant, and O₂ π* manifold. However, the synergistic effect of the transition metal substitution (Co@V_S) and biaxial strain on O₂ activation and dissociation in MoS₂ has not been systematically explored yet. Most previous strain engineering studies focus on the activation-related descriptors, whereas the kinetic feasibility of active-site regeneration following O₂ dissociation has received significantly less attention from researchers. Therefore, in this work, we investigated how biaxial tensile strain influenced the electronic structure of Co@V_S and its interaction with adsorbed O₂ species. By analyzing structural, vibrational, and electronic responses under strain, we clarified how lattice deformation affected the O–O activation process and associated dissociation pathway and further examined the energetics of active-site regeneration to provide a comprehensive view of the full catalytic cycle. The present analysis focuses on the intrinsic, potential-independent structural energetics of the activation–regeneration balance; the effects of electrode potential, solvation, and proton–electron transfer, which are important under realistic electrochemical conditions, are beyond the scope of this work.

2. Computational method

We performed spin-polarized DFT calculations to evaluate the structural stability of Co-doped MoS₂ monolayers (1H phase) under tensile strain and explore the bonding characteristics of adsorbed O₂ species, including the energetics of the O–O bond dissociation. All calculations utilized the plane-wave projector augmented-wave (PAW) method^33,34 implemented in the Vienna Ab initio Simulation Package (VASP)^35,36 with exchange–correlation effects described by the GGA–PBE functional.³⁷ A plane-wave kinetic energy cutoff of 500 eV was applied. Geometry optimizations proceeded until residual atomic forces were below 0.01 eV Å⁻¹, and the electronic energy convergence was set to 10⁻⁶ eV. Stress–strain calculations for the pristine MoS₂ unit cell were performed using an 8 × 6 × 1 Monkhorst–Pack³⁸ k-point mesh. In contrast, adsorption and nudged elastic band (NEB) calculations for the 5 × 3 supercell employed a 2 × 2 × 1 mesh. For partial density of states (PDOS) calculations, a denser 6 × 6 × 1 mesh was utilized, while a 3 × 3 × 1 mesh was applied for Bader charge and crystal orbital Hamilton population (COHP)³⁹ analyses. Geometry optimizations and NEB barrier calculations were performed within an orthorhombic supercell under equal biaxial in-plane strain. The orthorhombic 5 × 3 cell is a standard rectangular representation of the hexagonal MoS₂ lattice (a_orth = a_hex, b_orth = × a_hex). Because the two lattice vectors are related by a fixed geometric factor, uniform scaling by (1 + ε) preserves the hexagonal symmetry of the underlying lattice; thus, the strain condition ε_x = ε_y in the orthorhombic frame is equivalent to isotropic biaxial strain in the hexagonal frame. To facilitate a comparison with previous studies and avoid artifacts from Brillouin zone folding, phonon dispersions were computed for the hexagonal primitive cell under the same biaxial strain and plotted along the Γ–M–K–Γ path. As shown in Fig. 1(a), a rectangular unit cell of the MoS₂ monolayer was adopted and expanded into a 5 × 3 × 1 supercell containing 90 framework atoms. Co doping was performed by substituting a Co atom at a sulfur vacancy site (Co@V_S), corresponding to approximately 1.7% doping on the chalcogen sublattice. Subsequently, an O₂ molecule was adsorbed at this Co site, resulting in a 92-atom Co@V_S–O₂ system (containing 30 Mo, 59 S, 1 Co, and 2 O atoms; see Fig. 1(a)). Periodic boundary conditions were applied, and a vacuum spacing of approximately 15 Å was introduced along the surface normal (z direction) to prevent spurious interactions between periodically repeated slabs. The in-plane stress of the MoS₂ monolayer was derived from the cell-averaged three-dimensional (3D) virial stress and reported as a 3D in-plane stress, using a fixed effective thickness d = 0.615 nm (half the interlayer spacing of bulk MoS₂) for consistency with a prior work.²⁷ Note that the stress values expressed in GPa are conventional 3D-equivalent quantities that depend on the assumed effective thickness and therefore do not represent an intrinsic 2D mechanical property. All physically meaningful analyses in this work, including the S-preserved (SP)–S-reconstructed (SR) transition (hereafter referred to as SP–SR transition), dissociation barriers, and regeneration energetics, are parameterized by strain (a dimensionless, thickness-independent quantity) rather than by stress. For reference, 2D stress can be obtained as σ_2D (N m⁻¹) = σ_3D (GPa) × d, where d = 0.615 nm is the effective thickness adopted herein. Isotropic biaxial tensile strain (ε_x = ε_y) was imposed by uniformly scaling the in-plane lattice parameters by (1 + ε), followed by the full relaxation of all internal atomic coordinates at each strain.

Fig. 1 (a) (5 × 3) supercell structure of Co-doped monolayer MoS₂ with an O₂ molecule adsorbed at the substitutional Co site; the dotted box indicates the pristine unit cell. (b) In-plane stress–strain curve of the pristine unit cell recorded under incremental biaxial tensile loading. The stress values are conventional 3D-equivalent quantities obtained using an effective thickness of 0.615 nm; see the Computational Method section for details. (c) Phonon dispersions obtained at strains ε = 0.18 and 0.19. In (a), the purple, yellow, blue, and red spheres represent Mo, S, Co, and O atoms, respectively.

Vibrational properties were calculated using a finite displacement method and subsequently processed with Phonopy software^40,41 to obtain phonon dispersion curves. Bader charge analysis was conducted according to the standard partitioning scheme.^42–44 COHP curves were generated with LOBSTER software⁴⁵ utilizing appropriate PAW-based projection sets, and the projected COHP (pCOHP) was analyzed. Following the standard convention, the negative values −pCOHP were plotted so that the positive and negative peaks corresponded to the bonding and antibonding contributions, respectively. Minimum-energy reaction pathways, including the O–O dissociation barrier, were determined using the NEB method with the climbing-image variant (CI–NEB).^46–48 For each elementary step, five intermediate images were inserted between the optimized initial and final states. The images along the reaction path were optimized until the maximum force became less than or equal to 0.02 eV Å⁻¹. To assess the sensitivity of the main conclusions to the choice of exchange-correlation functional, PBE + U calculations were performed for the Co@V_S–O₂ system using an effective Hubbard parameter U_eff = 3.3 eV for the Co d electrons, consistent with the value employed by Oh et al.⁴⁹ for Co 3d electrons in TMD-based heterostructures, originally adopted from Grimaud et al.;⁵⁰ the Hubbard correction was applied solely to the Co dopant, as the localized d states at the substitutional site are the primary concern for self-interaction errors (Table S3 in the SI).

3. Results and discussion

Fig. 1(a) shows the 5 × 3 supercell of Co-doped monolayer MoS₂ with an O₂ molecule adsorbed at the substitutional Co site; the dotted box highlights the pristine MoS₂ unit cell. To examine the effect of strain on O₂ dissociation at this Co-doped site, we first applied incremental biaxial tensile strain to the pristine MoS₂ unit cell which has the in-plane lattice parameters a = 3.183 Å and b = 5.514 Å, and calculated the corresponding in-plane stress. The resulting biaxial stress–strain relationship is summarized in Fig. 1(b): as the tensile strain ε increases, the in-plane stress σ initially rises nonlinearly, reaches a peak, and then decreases. At ε = 0.22, the stress attains the maximum value σ_max = 24.54 GPa, representing the theoretical biaxial stress derived from the homogeneous σ–ε curve without considering lattice instabilities. However, a detailed evaluation of dynamical stability based on the phonon dispersions presented in Fig. 1(c) reveals that the pristine MoS₂ monolayer becomes unstable at a lower strain. No imaginary phonon frequencies are observed at ε = 0.18; however, at ε = 0.19, the out-of-plane acoustic (ZA) branch exhibits imaginary frequencies at the K point, indicating the onset of phonon instability. These results suggest that the critical strain for the onset of dynamical instability lies between ε = 0.18 and 0.19, corresponding to a stress of approximately 23.42 GPa. Therefore, we adopted this value as the ideal biaxial strength of the pristine MoS₂ monolayer limited by the onset of phonon instability, corresponding to a 2D stress of approximately 14.4 N m⁻¹. This estimate is in good agreement with the results of previous studies. Li et al.²⁷ reported an ideal biaxial strength of 23.8 GPa at ε ≈ 0.20 and identified the onset of phonon instability at the K point, while Cooper et al.⁵¹ constructed a biaxial stress–strain curve in quantitative agreement with the results reported by Li et al. The overall shape of the stress–strain curve obtained here, an initial nonlinear rise followed by a peak and subsequent softening, is consistent with both of these prior reports, and the phonon-instability threshold identified at the K point likewise agrees with the instability mode reported by Li et al., confirming the reliability of the present computational setup for describing the mechanical response of monolayer MoS₂ under biaxial tension. Furthermore, experimental studies have demonstrated that suspended monolayer MoS₂ membranes can sustain biaxial tensile strains of over 5% before rupture and exhibit breaking strains of 6–11% during nanoindentation.^29,52 These experimental findings motivated us to explore strains up to 10%, and our phonon-stability analysis confirmed that the biaxial strain window of 4.5–10% examined in the present work remained dynamically stable, providing a substantial safety margin relative to the ideal phonon-instability threshold. Additionally, phonon dispersion calculations conducted for monolayer MoS₂ containing a single sulfur vacancy and the corresponding Co-substituted structure at a 10% biaxial tensile strain revealed the absence of imaginary modes along the Γ–M–K–Γ path (Fig. S1). These results confirmed that the defected and Co-doped monolayers remained dynamically stable even at the highest strain applied in this work.

Using the strained pristine lattices obtained earlier as the starting points, Co-doped MoS₂ supercells were constructed next. For each strain, the relaxed primitive cell was expanded into a 5 × 3 supercell, where a single S atom was replaced with a Co atom to create a substitutional Co atom at a sulfur site (Co@V_S). The resulting Co-doped structure with an O₂ molecule adsorbed at the Co@V_S site is shown in Fig. 1(a); at each strain, all internal atomic coordinates were fully relaxed. The adsorption of O₂ at this Co site in 1H-MoS₂ resulted in two distinct local minima, which corresponded to the SP and SR configurations (Fig. 2(a)). At biaxial tensile strains up to approximately 6%, structural relaxations can converge to either the SP or SR local minimum. However, beyond a strain of approximately 6%, the system consistently converges to the SR configuration, and the SP structure can no longer be maintained. The most notable difference between these two configurations lies in the Co–S bond length. To compare them more directly, we focus on the 4.5% strain case in Table 1: the SP configuration retains a largely intact lower S sublayer and exhibits a relatively long Co–S separation of approximately 3.23 Å. In this state, the Co atom is not deeply embedded into the MoS₂ lattice but instead protrudes above the surface in an adatom-like geometry. In contrast, the SR configuration at the same strain demonstrates an upward displacement of the nearest-neighbor S atom and corresponding sinking of the Co atom into the MoS₂ monolayer. This results in a shortened Co–S bond length of approximately 2.19 Å with the Co atom more deeply embedded into the lattice. This substrate-mediated reconstruction significantly alters the ligand field and local electronic environment around the Co center, as further supported by the results of PDOS and COHP analyses discussed in the next paragraph. At biaxial tensile strains up to 4%, the SP structure is consistently more stable than the SR structure; however, the energy difference between them gradually decreases and reverses at a 4.5% strain, when the SR configuration becomes energetically favored (Fig. 2(b)). To illustrate this trend and the related transition point, Table 1 summarizes the geometric parameters of the lowest-energy configurations computed at all strains. At a 4.5% strain, the metastable SP configuration is also listed to highlight structural changes across the SP–SR transition. The O–O bond lengths reported in Table 1 already exceed 1.30 Å at both reference strains—0% (SP) and 4.5% (SR)—and the corresponding O–O stretching frequencies in Table S1 are strongly red-shifted from the experimental gas-phase value (1.207 Å and 1556 cm⁻¹⁵³) to approximately 1050 cm⁻¹ for the SP configuration and 980–990 cm⁻¹ for the SR configuration. This indicates that in both cases, the adsorbed O₂ molecules exhibit the characteristics of superoxo-like O₂ species. In transition-metal–dioxygen complexes, the simultaneous elongation of the O–O bond and reduction of the stretching frequency are characteristic signatures of superoxo-like O₂ species, which are commonly associated with the partial filling of the O–O π* antibonding orbitals—that is, a partially activated O₂ species.⁵⁴ According to Table 1, the O–O bond length within the SP branch (0–4% strain) increases gradually with biaxial tension and becomes further elongated in the SR configuration at a 4.5% strain, reflecting a systematic enhancement in the activation of the adsorbed O₂ species. When the strain increases from 4.0% to 4.5%, the additional elongation and activation of the O–O bond originate not only from the increased strain itself but, more importantly, from the accompanying transition from the SP to SR configuration. Thus, ε ≈ 4.5% marks a qualitative transition from a gently strain-activated SP branch to a reconstructed SR branch with distinctly stronger O₂ activation, providing a natural dividing line for the subsequent electronic structure and reaction-pathway analyses.

Fig. 2 (a) Side views of the SP and SR configurations of the O₂ molecule adsorbed at the Co substitutional site (Co@V_S) in monolayer MoS₂. (b) Relative energy difference between the SR and SP configurations (E_SR − E_SP), calculated per 5 × 3 supercell, as a function of biaxial strain. The energy crossover at ε = 4.5% indicates the strain at which the SR structure becomes thermodynamically favored. In (a), the purple, yellow, blue, and red spheres represent Mo, S, Co, and O atoms, respectively.

Table 1 Structural distances (in Å) for O₂ adsorbed at a substitutional Co site in monolayer MoS₂ under the biaxial tensile strain. The Co–O bond lengths, O–O separation, and Co–S distance between the Co atom and S atom located directly beneath it in the lower S sublayer are provided for each strain. Up to a 4% strain, the data correspond to the SP configuration, while at a 4.5% strain and above, they represent the SR configuration, which becomes energetically favorable at a higher strain. Additionally, the metastable SP structure at a 4.5% strain is included to illustrate the structural changes occurring during the SP–SR transition

Configuration	Strain	Co–O1	Co–O2	O1–O2	Co–S
SP	0.000	1.895	1.893	1.357	3.469
	0.020	1.878	1.883	1.360	3.367
	0.040	1.862	1.878	1.364	3.258
	0.045	1.857	1.877	1.364	3.230
SR	0.045	1.828	1.828	1.380	2.194
	0.075	1.826	1.826	1.374	2.148
	0.100	1.825	1.825	1.373	2.117

---

To further clarify these strain-induced structural changes, we compared the spin-resolved PDOS of the substitutional Co center, adsorbed O₂ molecule, and S atom located directly beneath the Co site, along with the corresponding −pCOHP spectra for the two local minima coexisting at the SP–SR crossover strain ε = 4.5% (Fig. 3). In the SP configuration at ε = 4.5% (Fig. 3(a) and (c)), the Co d manifold remains strongly spin-polarized, exhibiting large majority-spin peaks with energies several eV below the Fermi level and comparatively weak minority-spin intensity near the Fermi level. The p states of O1 and O2 exhibit only modest weights around the Fermi level, while the S p contribution is negligible, which is consistent with a superoxo-like O₂ species coupled to a high-spin Co center with a total magnetic moment of approximately 3 µ_B. From a ligand-field perspective, this SP configuration corresponds to a relatively weak ligand field around the Co center that stabilizes the high-spin state. The −pCOHP curves depicted in Fig. 3(c) support this interpretation: the O1–O2 bond shows relatively large positive (bonding) contributions that are well below the Fermi level along with the pronounced antibonding features near the Fermi level. Meanwhile, the Co–O1 and Co–O2 bonds are moderately bonding, and the Co–S interaction is nearly negligible in a good agreement with the large −ICOHP value obtained for O–O and much smaller values determined for Co–O and Co–S (Table 2). In the SR configuration at the same strain (Fig. 3(b) and (d)), the Co d states become more symmetric and broaden around the Fermi level, with reduced overall spin asymmetry. Additionally, both the O p and S p states gain intensity within the same energy range. This evolution indicates stronger Co–O and Co–S covalent interactions and partial delocalization of the Co local moment onto the ligands, which are consistent with the reduced total magnetic moment of approximately 1 µ_B. At the atomic level, the local magnetic moments on the two oxygen atoms are 0.372 and 0.343 µ_B for O1 and O2, respectively, in the SP configuration at ε = 4.5%, confirming that both atoms retain significant spin polarization predominantly in the p orbitals. In the SR configuration at the same strain, these values decrease to 0.257 and 0.257 µ_B, consistent with the progressive quenching of the open-shell character during the SP–SR transition. In other words, the shorter and more covalent Co–O and Co–S bonds in the SR configuration effectively increase coordination and strengthen the ligand field, driving a high spin-to-low spin transition at the Co site. Correspondingly, the −pCOHP curves presented in Fig. 3(d) reveal enhanced bonding contributions for Co–O and especially Co–S bonds, which are accompanied by a modest reduction in the integrated bonding strength of the O1–O2 bond. These observations agree with the configuration-dependent −ICOHP values listed in Table 2. Taken together, these results suggest the following causal chain connecting strain to the regeneration bottleneck. Biaxial tensile strain drives the SP-to-SR structural transition, in which the subsurface S atom approaches the Co center (Table 1: Co–S decreases from approximately 3.23 Å to 2.19 Å) and forms a new covalent Co–S bond (Table 2: −ICOHP increases from 0.097 to 2.038 eV). This increased coordination strengthens the ligand field around the Co site, as supported by the broadening and increased spin symmetry of the Co d states near E_F, the enhanced −ICOHP values for the Co–O bonds, and the reduction of the total magnetic moment from approximately 3 µ_B to 1 µ_B, consistent with partial spin-state reduction. The resulting more covalent and more tightly coordinated Co center simultaneously strengthens the Co–O interaction and deepens the post-dissociation Co–O–Mo product well. This strengthened bonding provides an electronic-structure basis for the regeneration penalty discussed later: the same ligand-field enhancement that promotes O₂ activation also anchors the dissociated oxygen more tightly to the Co center. According to Tables 1, 2 and S1, the observed electronic-structure trends indicate the occurrence of a charge transfer into the antibonding O–O states at the Co@V_S site and produce a coherent picture of the stronger O₂ activation in the SR configuration as compared with the SP state at ε = 4.5%. These results are consistent with the widely accepted mechanism of O₂ activation via the back-donation into antibonding O–O π* orbitals. Although our PDOS and COHP analysis results do not explicitly resolve the orbital symmetry, the observed trends strongly support this charge-transfer model.

Fig. 3 Spin-resolved PDOS and −pCOHP analyses of the O₂ molecule adsorbed on Co@V_S under a biaxial tensile strain of ε = 4.5%. The spin-resolved PDOS calculated for the (a) SP and (b) SR configurations, with the positive and negative values indicating the spin-up and spin-down components, respectively. The Fermi level is set at zero energy. The −pCOHP curves obtained for selected Co–O, O–O, and Co–S bonds in the (c) SP and (d) SR configurations; here, the positive and negative values correspond to the bonding and antibonding contributions, respectively.

Table 2 Integrated −ICOHP (eV per bond) values for selected bonds surrounding adsorbed O₂ in the SP and SR configurations at ε = 4.5%. The higher –ICOHP values correspond to the stronger covalent interactions

Configuration	Co–O1	Co–O2	O1–O2	Co–S
SP	2.145	1.889	11.354	0.097
SR	2.368	2.373	10.776	2.038

---

Before analyzing the strain dependence of O–O dissociation, we first benchmarked our calculations against the reaction pathway previously reported by Zhao et al. for the O₂ molecule adsorbed on a Co-doped MoS₂ monolayer without strain.¹⁵ Along this direct dissociation pathway, the O–O bond is cleaved, and the two oxygen atoms ultimately occupy separate adsorption sites: one atom is bound to the substitutional Co atom, and the other atom is attached to a neighboring S atom in close agreement with the structural motif proposed by Zhao et al. As shown in Fig. 4(a), our CI–NEB calculation yields an activation barrier of approximately 1.30 eV, which is nearly identical to the literature value of ∼1.33 eV, confirming that our computational setup reliably reproduces the established dissociation channel. However, upon the application of biaxial tensile strain, the minimum-energy path progressively evolves toward a distinct, Mo-assisted dissociation pathway. In this Mo-assisted dissociation pathway, the O–O bond is cleaved while both oxygen atoms remain coordinated to the Co center; meanwhile, one oxygen atom forms an additional bridge to the underlying Mo atom, resulting in a substantially more exothermic Co–O–Mo termination (≈0.6 eV more stable than the direct Co/S termination at ε = 0%), as illustrated in Fig. 4(b). At ε = 0%, the Mo-assisted pathway exhibits an activation barrier that is essentially identical to that of the direct dissociation pathway (∼1.30 eV for both pathways); however, it leads to a much more exothermic Co–O–Mo final state (ΔE ≈ −0.93 eV compared with −0.35 eV for the direct pathway; see Fig. 4). Consistently, the total magnetic moment along the Mo-assisted pathway decreases from ∼3 µ_B in the initial superoxo-like state to approximately 1 µ_B in the Co–O–Mo final state, whereas the direct pathway retains a larger moment (∼3 µ_B) in the Co/S termination, indicating a more strongly quenched, oxide-like Co–O–Mo configuration in the Mo-assisted case. It is worth noting that the gas-phase O₂ molecule was computed in its spin-polarized triplet configuration, consistent with its experimental ground state. The total magnetic moment of approximately 3 µ_B in the initial adsorbed state is consistent with a superoxo-like species retaining substantial open-shell character and ferromagnetic alignment with the Co center. The progressive reduction to approximately 1 µ_B in the Co–O–Mo product state reflects the gradual quenching of this open-shell character as the O–O bond is cleaved, and the oxygen atoms become incorporated into the lattice. Spin-polarized calculations were employed throughout all NEB pathways, allowing the magnetic moment to evolve self-consistently along the reaction coordinate. Because this Mo-assisted final state also becomes the preferred product under the finite strain and is closely tied to the strain-induced reconstruction of the Co@V_S site, we focus on this new pathway in the following analysis of strain-dependent barriers and the electronic structure. From a catalytic perspective, this alternative Mo-assisted channel is especially important because it features an activation barrier similar to that of the direct pathway, while producing a significantly deeper Co–O–Mo product well. This indicates that under the applied strain, the Mo-assisted pathway becomes the thermodynamically favored channel, effectively promoting the dissociation process and strongly stabilizing the oxygen species formed on the catalyst surface.

Fig. 4 Minimum-energy pathways determined for the O₂ dissociation on Co-doped MoS₂ at a zero strain via CI–NEB calculations. (a) “Direct dissociation” and (b) “Mo-assisted dissociation” pathways.

The subsequent regeneration analysis focuses on the 4.5–5.5% strain range, as this window encompasses three qualitatively distinct post-dissociation outcomes observed with increasing strain: (i) at ε ≤ 4% (SP branch), O₂ dissociation proceeds along a well-defined pathway with a progressively decreasing barrier (1.30 to 1.04 eV; Fig. 5(a)); (ii) at ε = 4.5–5.0% (SR branch), the dissociation barrier decreases further to approximately 0.68 eV, but the post-dissociation oxygen species penetrate subsurface sites, displacing the Co dopant and causing irreversible structural disruption (Fig. S2), with an energetic penalty of approximately +1.15 to +1.23 eV required to escape this penetrated trap (Fig. S3); and (iii) at ε = 5.5%, a distinct metastable Co-top O intermediate becomes stable for the first time, enabling the construction of a well-defined regeneration pathway. This 4.5–5.5% window, therefore, represents the narrowest range within which the activation–regeneration trade-off transitions from site destruction to a quantifiable kinetic bottleneck, making it the most informative strain window for the subsequent analysis.

Fig. 5 (a) Activation energy obtained for O₂ dissociation along the Mo-assisted pathway as a function of biaxial tensile strain for the SP and SR configurations. The vertical dashed line at ε = 4.5% indicates the SP–SR crossover point. (b) Reaction energy ΔE(ε) = E_FS(ε) − E_IS(ε) per 5 × 3 MoS₂ supercell, where E_IS and E_FS represent the total energies of the initial adsorbed O₂ state and final Co–O–Mo termination at each strain, respectively. The negative values signify the exothermic formation of the Co–O–Mo state. In (b), the inset shows the total magnetic moment (µ_B) of the initial adsorbed O₂ state as a function of strain; the magnetization decrease at the higher strain is consistent with a strain-induced spin reorganization that may contribute to the kink in ΔE(ε) near ε ≈ 7%. A more pronounced change in the total magnetic moment is also observed for the final Co–O–Mo termination beyond ε ≈ 7% (not shown here), which is consistent with the same kink in ΔE(ε).

To quantify how biaxial tensile strain influences the activation of adsorbed O₂ molecules, we calculated minimum-energy paths for O–O dissociation using the climbing-image NEB method. Along the SP branch (ε ≤ 4%), the dissociation barrier decreases almost linearly from 1.297 eV at ε = 0% to 1.036 eV at ε = 4%, demonstrating that strain facilitates the O–O bond cleavage (Fig. 5(a)). This reduction coincides with a gradual elongation of the O–O bond from 1.357 to 1.364 Å and slight shortening of the Co–O distances (Table 1), consistent with the gradual weakening and activation of the intramolecular O–O bond under tensile strain, as commonly observed for transition metal–dioxygen complexes.^53–55 A qualitatively different regime emerges at ε = 4.5%, where the system transitions from the SP to SR configuration. In the SR configuration (ε = 4.5%), the O₂ dissociation barrier decreases to 0.682 eV and remains at a sub-0.7 eV level, reaching ∼0.67 eV at ε = 5.5% (Fig. 5(a)). This reduction is associated with the (i) formation of shorter Co–O bonds with lengths of approximately 1.83 Å and (ii) further elongation of the O–O bond to 1.380 Å as compared with that in the SP configuration. Bader charge analysis (Table 3) shows that the net charge on the O₂ moiety remains around ∼0.6–0.7 e across the studied strain range. However, the SR configuration at ε = 4.5% exhibits pronounced intramolecular charge asymmetry (Δq(O1) ≈ 0.54 e and Δq(O2) ≈ 0.18 e), indicating an electronically asymmetric superoxo-like O₂ species with the charge preferentially localized on one oxygen atom. Along the SP branch (ε ≤ 4%), the Co center is gradually reduced as Δq(Co) increases from −0.37 e to −0.17 e. Upon transitioning to the SR state at ε = 4.5%, Δq(Co) becomes more negative again (≈−0.26 e), indicating a partial reoxidation of Co accompanied by the pronounced intramolecular charge asymmetry within O₂. Although the total Bader charge on O₂ increases moderately from approximately 0.58 e at ε = 4% (SP) to 0.72 e at ε = 4.5% (SR), this change is small compared with the large charge imbalance between O1 and O2. This suggests that the SP–SR transition primarily redistributes charges between Co and the two O atoms rather than significantly altering the overall reduction level of O₂. These findings reveal two distinct activation regimes: at ε ≤ 4%, biaxial tension induces continuous mechanochemical activation within the SP configuration, whereas at ε = 4.5%, a substrate-mediated reconstruction to the SR configuration creates a new initial state, in which cooperative structural and electronic rearrangements synergistically lower the O₂ dissociation barrier. In addition to the barrier heights, we examined the strain dependence of the reaction energy, ΔE(ε) = E_FS(ε) − E_IS(ε), per 5 × 3 supercell (Fig. 5(b)), where E_IS and E_FS denote the total energies of the initial adsorbed O₂ state and final Co–O–Mo termination along the same Mo-assisted pathway at each strain, respectively. Along the SR branch, the reaction energy first gradually increases (becomes less exothermic) between ε = 4.5% and 6.5% and then exhibits a small kink around ε ≈ 7% (Fig. 5(b)), after which it decreases again toward more exothermic values at higher strains. At this point, the total magnetic moment of the final Co–O–Mo state changes from approximately 1 µ_B at lower strains to approximately 0.4 µ_B at ε ≈ 7%. This spin-state rearrangement slightly modifies the total energy of the product state, producing the kink depicted in Fig. 5(b). Nevertheless, O₂ dissociation remains strongly exothermic (ΔE ≈ −1.6 to −1.8 eV) along the entire SR branch, reinforcing the conclusion that the Mo-assisted pathway leads to a deeply stabilized Co–O–Mo termination under strain. From a catalytic perspective, such a deep and strongly exothermic Co–O–Mo termination suggests that the Co@V_S site is likely thermodynamically trapped in this O-rich state under strain. An analysis of the regeneration pathway (see below) shows that this excessive stabilization imposes a significant energy penalty for removing the oxygen species, thereby hindering the reactivation process. Overall, these findings reveal two distinct regimes of O2 activation at Co@VS: a strain-driven mechanochemical activation along the SP branch (ε ≤ 4%), and a reconstruction-driven regime for ε ≥ 4.5%, in which the SR configuration and resulting Co–O–Mo termination work synergistically to lower the dissociation barrier while strongly stabilizing the dissociated oxygen species through the formation of robust Co–O and Mo–O bonds.

Table 3 Strain-dependent Bader charge variations (Δq) obtained for the substitutional Co center and two O atoms in the adsorbed O₂ molecule on Co-doped MoS₂. Here, Δq is defined relative to the neutral atomic valence charges (Co: 9e, O: 6e) with the positive values indicating electron accumulation and negative values indicating electron depletion. For strains ε ≤ 4%, the data correspond to the SP configuration, while for ε ≥ 4.5%, they represent the SR configuration. At ε = 4.5%, both the metastable SP and thermodynamically favored SR structure are included for comparison

Configuration	Strain	Δq(Co)	Δq(O1)	Δq(O2)
SP	0.000	−0.368	0.331	0.325
	0.020	−0.370	0.324	0.322
	0.040	−0.173	0.295	0.282
	0.045	−0.061	0.294	0.271
SR	0.045	−0.263	0.541	0.180
	0.075	−0.337	0.542	0.043
	0.100	−0.348	0.552	0.068

---

It was found that the Mo-assisted dissociation pathway led to a highly stabilized Co–O–Mo termination. Notably, additional structural optimizations conducted at lower strains (ε = 4.5–5.0%) demonstrated that the remaining oxygen atom tended to diffuse into the MoS₂ lattice, displacing the Co dopant and resulting in the formation of a disordered oxide-like structure (Fig. S2). Furthermore, escaping from this penetrated Co–O–Mo trap to representative “recovery” endpoints is a strongly endothermic process (ΔE ≈ +1.15 eV for Co-top; see Fig. S3), indicating a deep thermodynamic trap and the existence of substantial kinetic barriers. Because such structural disruption complicates the recovery of the original single-atom site, we focused our regeneration analysis on the higher-strain regime, in which the active site remained intact. Experimentally, biaxial tensile strains up to approximately 5.6% have been achieved in suspended monolayer MoS₂ membranes through pressure-induced bulging, establishing a practical feasibility limit for strain engineering in free-standing geometries.²⁹ Consequently, we concentrated on regeneration kinetics and energetics within this feasible range (ε ≈ 5–6%). Here, we defined active-site regeneration as the removal or relocation of dissociation-derived O species from the Co center via lateral lattice-oxygen migration, which serves as a proxy for the intrinsic difficulty of recovering the Co coordination environment. Although this proxy does not explicitly account for the electrode potential, solvent effects, or proton–electron transfer, it provides a consistent measure of the intrinsic kinetic barrier to reopening the Co coordination environment after O₂ dissociation—an essential prerequisite for the sustained catalytic turnover. To evaluate this proxy, we constructed a well-defined minimum-energy pathway from the results of CI–NEB calculations, consisting of Steps I–III (Fig. 6(a)) and Steps IV–VI (Fig. 6(b)). We selected ε = 5.5% (SR) as a representative condition, where the active site remains intact and O₂ dissociation is strongly facilitated (E_a for O₂ dissociation ≈ 0.67 eV). At lower strains (ε = 4.5–5.0%), the Step III endpoint (Co-top O intermediate) cannot be stabilized, preventing the formation of a well-defined intermediate that bridges the subsequent lateral-migration steps. Therefore, we focus on ε = 5.5%, at which the Step III intermediate becomes a distinct metastable minimum, allowing the regeneration pathway to be consistently defined. Accordingly, we constructed the regeneration pathway starting from the fully dissociated Co–O–Mo product state at ε = 5.5% in the SR configuration (Fig. 6). Initially, the two oxygen atoms form an O–Co–O–Mo configuration: one O atom bridges Co and a neighboring S atom, while the other O atom binds to both Co and an underlying Mo atom. This minimum-energy path involves six elementary oxygen-migration steps (Steps I–VI). The final state of each step is used as the initial state of the subsequent step. In the first step, the outer O atom breaks its Co–O bond and migrates to a nearby S adsorption site, creating a terminal O atom bound to a surface S site and leaving a single Co–O–Mo linkage. In the second step, the remaining O atom detaches from the Mo site and moves upward toward the Co center, forming a Co–O species oriented closer to the surface normal at the Co-top site. In the third step, the Co-bound oxygen undergoes further rearrangement along the minimum-energy path, maintaining the Co-top O motif before the lateral diffusion. This Step III intermediate represents a distinct strain-dependent metastable minimum that becomes stable only at ε ≥ 5.5%. At lower strains (ε = 4.5–5.0%), it relaxes into a different oxide-like configuration, consistent with the active-site disruption behavior discussed above. Therefore, Step III acts as a preparatory rearrangement that preserves the Co-top O motif and establishes a local geometry for the subsequent lateral migration steps (Steps IV–VI). In the fourth step, the Co-bound oxygen moves (hops) from the Co-top position toward a neighboring S site, forming an intermediate configuration in which the O atom is laterally displaced away from the Co center but continues to interact with it. During the fifth step, the Co-bound oxygen undergoes the rate-limiting lateral migration away from the Co center toward the neighboring S site, substantially weakening—but not completely eliminating—the residual Co–O interaction. Finally, in the sixth step, this remaining weak Co–O interaction is fully released, resulting in a clean Co@VS site coordinated solely by the MoS₂ lattice, with two terminal O-on-S species nearby. The energy profile presented in Fig. 6 indicates that the initial migration steps (Steps I–III) have moderate activation barriers of 0.39, 0.41, and 1.09 eV, respectively. The subsequent oxygen-removal stage (Steps IV–VI) features a moderate barrier in Step IV (0.42 eV), high barrier in Step V (1.63 eV), and small final barrier in Step VI (0.17 eV). Importantly, in Step II, the transition and final states are nearly isoenergetic, resulting in a small reverse barrier (E_TS − E_FS) and thus a non-negligible probability of back-migration under thermal fluctuations. Additionally, the rate-determining barrier in Step V (1.63 eV) is more than two times greater than that of the unstrained system (∼0.77 eV; Fig. S4), indicating that the strain-induced stabilization of the Co–O bond substantially hinders oxygen diffusion. This finding is in good agreement with the larger (more bonding) −ICOHP values for the Co–O bond listed in Table S2. Furthermore, the overall regeneration process is strongly endothermic, with the final state ∼2.20 eV higher than the initial Co–O–Mo state, reinforcing the thermodynamic penalty that deepens the kinetic bottleneck in Step V. Hence, in the lattice-oxygen migration proxy pathway considered herein, tensile strain strengthens the Co–O interaction and elevates the barrier for the lateral oxygen diffusion, thereby hindering the regeneration of the clean Co@V_S site even after the strain-enhanced O₂ activation. Consistent with the deep thermodynamic trapping illustrated in Fig. 5, this phenomenon may imply persistent oxygen coverage and highlights a critical trade-off between facile O₂ dissociation and challenging site reactivation in strain-engineered single-atom catalysts.

Fig. 6 Minimum-energy pathway obtained for regenerating the clean Co@V_S site from the Co–O–Mo-terminated intermediate under a biaxial tensile strain of ε = 5.5% (SR) by performing CI–NEB calculations. (a) Co-site rearrangement following O₂ dissociation (Steps I–III). (b) Oxygen hopping from the Co site to a neighboring S site (Steps IV–VI). Energies are referenced to the initial state of Step I (E = 0) with the final state of each step serving as the initial state for the subsequent step. The activation barriers determined for Steps I–VI via the formula E_a = E_TS − E_IS are 0.39, 0.41, 1.09, 0.42, 1.63, and 0.17 eV, respectively. Representative initial, transition, and final structures are shown for each step.

4. Conclusions

In this study, we performed spin-polarized DFT calculations to explore how strain affected the O₂ activation and regeneration mechanisms in Co-doped MoS₂ (Co@V_S). The obtained results revealed that applying biaxial tensile strain triggered a structural transition from the SP to an SR configuration at a strain of approximately 4.5%, serving as a critical switch that enhanced O₂ activation. Although the strain-induced SR configuration significantly lowered the dissociation barrier via a Mo-assisted pathway, catalyst regeneration remains a significant challenge owing to the concurrent strengthening of the Co–O interaction under the applied strain. Specifically, we identified two strain-dependent regeneration failure modes: at ε = 4.5–5.0%, the dissociation-derived oxygen species penetrate subsurface sites, displacing Co atoms and hindering the site recovery; at ε = 5.5%, the Co@V_S structure remains intact, but the oxygen removal process is limited by the lateral diffusion barrier of 1.63 eV (Step V), which is more than twice as high as the barrier in the unstrained case (∼0.77 eV; Fig. S4). This diffusion bottleneck is followed by the small final barrier of 0.17 eV (Step VI). These findings highlight a notable trade-off: while tensile strain effectively promotes O₂ dissociation, it also stabilizes the oxide-like product states and increases the kinetic difficulty of the oxygen removal process, making active-site regeneration a likely bottleneck within the regeneration (lattice-oxygen migration) pathway studied in this work. More broadly, this activation–regeneration trade-off may apply to other strain-engineered single-atom catalysts, for which O₂ activation correlates with stronger metal–oxygen bonding. Accordingly, our study suggests that future catalyst designs should carefully balance activation kinetics with the kinetic feasibility of regeneration to overcome this intrinsic limitation. It should be noted that the lattice-oxygen migration proxy adopted herein does not incorporate electrode potential, solvation, or proton–electron transfer steps; therefore, the absolute regeneration barriers may shift under realistic electrochemical conditions. PBE + U sensitivity tests (U_eff = 3.3 eV; Table S3) indicate that the SP–SR crossover strain shifts to higher values with the Hubbard correction, while the SR metastable minimum and the monotonic strain-dependent trend are preserved, supporting the qualitative robustness of the conclusions. Nevertheless, the qualitative activation–regeneration trade-off identified in this work may persist, as it originates from the fundamental strain-induced strengthening of Co–O and Mo–O bonds. Under applied electrochemical potential, the absolute regeneration barriers may shift in either direction depending on the specific potential, solvation environment, and proton-coupled reaction steps; thus, the barrier values reported here should be interpreted as baseline structural references rather than operando predictions. Coupling the present structural framework with explicit electrochemical models, such as the computational hydrogen electrode or implicit solvation approaches, represents an important direction for future work.

Finally, achieving uniform biaxial tensile strain under catalytic operating conditions remains experimentally challenging, although strains of approximately 5.6% have been demonstrated in suspended monolayer MoS₂ membranes²⁹ and breaking strains of 6–11% have been reported during nanoindentation.⁵² The present study, therefore, serves as a computational reference framework that identifies the mechanistic transitions governing the activation–regeneration balance; the insights may remain qualitatively relevant for any strain modality that produces comparable lattice expansion near the active site. The electronic-structure analysis relied on PDOS, COHP, and Bader charge methods, which provide atom- and bond-resolved insight into localized defect states; complementary band structure calculations could further clarify the dispersion of these states and represent a valuable direction for future investigation.

Author contributions

Soon-Dong Park: conceptualization, methodology, software, investigation, visualization, writing—original draft. Sung Youb Kim: writing—review & editing, supervision, project administration, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this study—including NEB pathways, energies, and numerical data for plots—are provided in the supplementary information (SI). Atomic structure files, additional input files, and analysis scripts can be obtained from the corresponding author upon reasonable request. References cited in the SI are included in the article's reference list. Supplementary information: Fig. S1–S4, Tables S1–S3. See DOI: https://doi.org/10.1039/d6ta01659f.

Acknowledgements

We gratefully acknowledge support from the Basic Research Laboratory program of the National Research Foundation of Korea (Grant No. 2021RA4A1033224) and the UNIST Supercomputing Center for providing access to supercomputing resources. This work was also supported by the InnoCORE program of the Ministry of Science and ICT (N10250154). The authors reviewed and edited the content thoroughly and took full responsibility for the final manuscript.

References

1. X. Chia, A. Y. S. Eng, A. Ambrosi, S. M. Tan and M. Pumera, Chem. Rev., 2015, 115, 11941–11966.
2. S. Zhao, K. Wang, X. Zou, L. Gan, H. Du, C. Xu, F. Kang, W. Duan and J. Li, Nano Res., 2019, 12, 925–930.
3. N. Karmodak, L. Bursi and O. Andreussi, J. Phys. Chem. Lett., 2022, 13, 58–65.
4. Y. Yan, S. Liang, X. Wang, M. Zhang, S. M. Hao, X. Cui, Z. Li and Z. Lin, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2110036118.
5. M. Wang, X. Huang, Z. Yu, P. Zhang, C. Zhai, H. Song, J. Xu and K. Chen, Nanomaterials, 2022, 12, 4069.
6. S. Xin, Z. Liu, L. Ma, Y. Sun, C. Xiao, F. Li and Y. Du, Nano Res., 2016, 9, 3795–3811.
7. S. Pan, Z. Cai, Y. Duan, L. Yang, B. Tang, B. Jing, Y. Dai, X. Xu and J. Zou, Appl. Catal., B, 2017, 219, 18–29.
8. J. Tang, C. Wang, H. Zhang and J. Guo, J. Alloys Compd., 2022, 911, 164991.
9. J. H. Kim, J. G. Lee and M. J. Choi, Materials, 2024, 17, 4277.
10. L. Fan, X. Dai, F. Li, X. Li, Z. Liu, Q. Guo, C. Zhang, Z. Kang and D. Sun, J. Mater. Sci. Technol., 2025, 226, 1–11.
11. L. Zhang, Y. Dong, L. Li, Y. Shi, Y. Zhang, L. Wei, C. L. Dong, Z. Lin and J. Su, Nano-Micro Lett., 2025, 17, 88.
12. D. Xu, J. Jiang, D. Liu, X. Wei, J. Yang, Z. Zhuang, Y. Dou, X. Liu, J. Duan and D. Wang, Microstructures, 2026, 6, 2026020.
13. C. Ataca and S. Ciraci, J. Phys. Chem. C, 2011, 115, 13303–13311.
14. D. Ma, W. Ju, T. Li, X. Zhang, C. He, B. Ma, Y. Tang, Z. Lu and Z. Yang, Appl. Surf. Sci., 2016, 364, 181–189.
15. B. Zhao, L. L. Liu, G. D. Cheng, T. Li, N. Qi, Z. Q. Chen and Z. Tang, Mater. Des., 2017, 113, 1–8.
16. W. Zhou, X. Zou, S. Najmaei, Z. Liu, Y. Shi, J. Kong, J. Lou, P. M. Ajayan, B. I. Yakobson and J. C. Idrobo, Nano Lett., 2013, 13, 2615–2622.
17. S. Najmaei, Z. Liu, W. Zhou, X. Zou, G. Shi, S. Lei, B. I. Yakobson, J. C. Idrobo, P. M. Ajayan and J. Lou, Nat. Mater., 2013, 12, 754–759.
18. S. Kc, R. C. Longo, R. M. Wallace and K. Cho, J. Appl. Phys., 2015, 117, 135301.
19. A. Wu, Q. Song and H. Liu, Comput. Theor. Chem., 2020, 1187, 112906.
20. H. Lu, A. Kummel and J. Robertson, APL Mater., 2018, 6, 066104.
21. X. Zhang, J. Xu, A. Zhi, J. Wang, Y. Wang, W. Zhu, X. Han and K. Wang, Adv. Sci., 2024, 11, 2408640.
22. K. Dolui, I. Rungger, C. Das Pemmaraju and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 075420.
23. H. P. Komsa and A. V. Krasheninnikov, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 125304.
24. M. Li, T. Li and Y. Jing, RSC Adv., 2022, 12, 31525–31534.
25. G. Liu, A. W. Robertson, M. M. Li, W. C. H. Kuo, M. T. Darby, M. H. Muhieddine, Y. Lin, K. Suenaga, M. Stamatakis and J. H. Warner, Nat. Chem., 2017, 9, 810–816.
26. W. T. Yein, D. Kim and Q. Wang, Chem. Proc., 2024, 15, 3.
27. T. Li, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 235407.
28. H. J. Conley, B. Wang, J. I. Ziegler, R. F. Haglund, S. T. Pantelides and K. I. Bolotin, Nano Lett., 2013, 13, 3626–3630.
29. D. Lloyd, X. Liu, J. W. Christopher, L. Cantley, A. Wadehra, B. L. Kim, B. B. Goldberg, A. K. Swan and J. S. Bunch, Nano Lett., 2016, 16, 5836–5841.
30. R. Frisenda, M. Drüppel, R. Schmidt, S. Michaelis de Vasconcellos, D. Perez de Lara, R. Bratschitsch, M. Rohlfing and A. Castellanos-Gomez, npj 2D Mater. Appl., 2017, 1, 10.
31. W. S. Yun, S. W. Han, S. C. Hong, I. G. Kim and J. D. Lee, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 033305.
32. S.-D. Park and S. Y. Kim, Coupled Syst. Mech., 2016, 5, 305–314.
33. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
34. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
35. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
36. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
38. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
39. R. Dronskowski and P. E. Bloechl, J. Phys. Chem., 1993, 97, 8617–8624.
40. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
41. A. Togo, L. Chaput, T. Tadano and I. Tanaka, J. Phys. Condens. Matter, 2023, 35, 353001.
42. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
43. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, J. Comput. Chem., 2007, 28, 899–908.
44. W. Tang, E. Sanville and G. Henkelman, J. Phys. Condens. Matter, 2009, 21, 084204.
45. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37, 1030–1035.
46. H. Jónsson, G. Mills and K. W. Jacobsen, in Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific, 1998, pp. 385–404.
47. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
48. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
49. N. K. Oh, J. Seo, S. Lee, H. J. Kim, U. Kim, J. Lee, Y. K. Han and H. Park, Nat. Commun., 2021, 12, 4606.
50. A. Grimaud, O. Diaz-Morales, B. Han, W. T. Hong, Y. L. Lee, L. Giordano, K. A. Stoerzinger, M. T. M. Koper and Y. Shao-Horn, Nat. Chem., 2017, 9, 457–465.
51. R. C. Cooper, J. W. Kysar and C. A. Marianetti, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 167401.
52. S. Bertolazzi, J. Brivio and A. Kis, ACS Nano, 2011, 5, 9703–9709.
53. P. H. Krupenie, J. Phys. Chem. Ref. Data, 1972, 1, 423–534.
54. C. J. Cramer, W. B. Tolman, K. H. Theopold and A. L. Rheingold, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 3635–3640.
55. M. M. Montemore, M. A. Van Spronsen, R. J. Madix and C. M. Friend, Chem. Rev., 2018, 118, 2816–2862.

---