Strain-controlled spin regulation in Fe–N–C catalysts for enhanced oxygen reduction reaction activity†

Mingyuan Yu^abc, Jiaxiang Wu^abc, Yashi Chen^d, Yongping Du^abc, Ang Li^abc, Erjun Kan*^abc and Cheng Zhan*^abc
aSchool of Physics, Nanjing University of Science and Technology, Nanjing 210094, China
^bMIIT Key Laboratory of Semiconductor Microstructure and Quantum Sensing, Nanjing University of Science and Technology, Nanjing 210094, China
^cEngineering Research Center of Semiconductor Device Optoelectronic Hybrid Integration in Jiangsu Province, Nanjing 210094, China
^dSchool of Materials Science and Engineering, Nanjing University of Science and Technology, Nanjing 210094, China

First published on 7th August 2024

---

---

Introduction

Efficient energy conversion is the key principle in the design of electrocatalysts for renewable energy applications.^1–4 Single-atom-catalysts (SACs) with isolated single metal sites dispersed on substrates have become a new frontier in catalysis science, owing to their high efficiency of atomic utilization and activity.^5–10 Among different SACs, Fe-embedded N-doped graphene (Fe–N–C) is the most representative case that has shown excellent ORR activity in fuel cell catalysis.^11–15 However, partial occupation of d orbitals in Fe could bring diverse spin states in electrocatalysis, which make mechanistic understanding more difficult and challenging.^16–18

In recent years, the spin effects in SACs have shown significant impact on the catalytic mechanism and activity, which has been both experimentally and theoretically studied.^19–21 For instance, Yang et al. introduced Mn single atoms around Fe sites in the Fe–N–C system by pre-polymerization and pyrolysis processes, which caused the Fe³⁺ spin state transition from the low spin (LS) to intermediate spin (IS) state, leading to the half-wave potential of 0.93 V at 0.1 M KOH.²² Zhang et al. also observed the spin state transition of Fe from LS to IS in Fe–N–C by introducing Fe atomic clusters, leading to the spin electron filling in the Fe–O σ* orbital and favors OH⁻ desorption.²³ However, Cheng et al. and Hou et al. found that Fe could change from HS to IS by introducing Fe and Ru clusters and they proposed that spin transition enhances the interaction with the π* antibonding orbitals of oxygen.^24,25 Zhong et al. reported a nearly linear relationship between magnetism and the oxygen reduction reaction (ORR) activity of an Fe single atom supported on C₂N (C₂N–Fe), suggesting that the electronic spin is a promising descriptor for Fe single-atom-based electrocatalysts.²⁶ Hu et al. found that the ORR catalytic performance of Fe–N–C is related to the spin state, which can be controlled by external fields.²⁷ Our previous work revealed the impact of spin multiplicity on ORR/OER catalysis and proposed the spin-switching pathway with the consideration of potential-dependent spin crossover.²⁸ This work also emphasized that the spin state preference is jointly determined by the valence (electron number) and ligand (coordination number), but the latter one is more dominative in spin manipulation.

According to the crystal field theory, the most preferable spin state of a complex is determined by its coordination environment and metal–ligand (M–L) interaction.^18,29 Thus, to manipulate the spin state, extra ligand introduction and strain engineering are the commonly used strategies. Xin et al. investigated the electronic and magnetic properties of 3d transition metals in TMN₄ SACs using density functional theory calculations and found two different spin states that can be switched over via external tensile strain.³⁰ Wang et al. demonstrated that the introduction of XO₂ groups with different polarity (SeO₂ > TeO₂ > SO₂) into the FeN₄ site leads to the spin electron transfer and spin state change due to the symmetry breaking.³¹ However, early research mainly focuses on the spin state transition and the effect of magnetic moment on activity, while the origin of diverse spin states was rarely discussed in early studies. In particular, the correlation between spin multiplicity, orbital hybridization, coordination geometry and frontier orbital information in SACs and reaction intermediates is not well understood at the current stage.

To reveal the origin of spin multiplicity in single atom electrocatalysis, in this work we take Fe–N–C as an example and employ density functional theory (DFT) to calculate the high-spin (HS) and intermediate-spin (IS) structures of all ORR intermediates. We analyze the bonding modes of adsorbates, the orbital occupation of spin states through the projected density of state (PDOS), and the bonding information by Crystal Orbital Hamilton Population (COHP) to explain the observed magnetic moment of each ORR intermediate and its corresponding spin state. In addition, we also investigate how the spin state is modulated by mechanical strain and its impact on the ORR activity, which brings a feasible way to improve the electrocatalytic activity via spin manipulation.

Methods

Spin-polarized DFT calculations were performed by using the Vienna ab initio simulation package (VASP).^32–34 In our DFT calculation, the exchange-correlation effect was treated by the general gradient approximation with the Perdew–Burke–Ernzerhof formalism (GGA-PBE).³⁵ The van der Waals interaction of adsorbates was described by the empirical correction in Grimme's scheme (DFT-D3).³⁶ The DFT + U method was employed to deal the strongly correlated 3d orbital of Fe and the value of U_eff = U–J was set to 4.03 eV according to our previous research.²⁸ The electron–nuclei interaction was described using the Projector-Augmented Wave (PAW) pseudopotential.³² The energy cutoff for the plane-wave basis was set to 550 eV. The k-point sampling in the Brillouin zone was implemented with a gamma-centered 2 × 2 × 1 mesh for geometry relaxation and 10 × 6 × 1 mesh for electronic structure calculation, which has been fully tested for energy convergence in our previous study.²⁸ The Fe–N–C SAC was modeled by a FeN₄ site embedded into a monolayer graphene in an orthorhombic cell with the lattice parameters of a = 9.94 Å and b = 17.04 Å, and the surface was modeled in contact with a 30 Å vacuum layer. The convergence criteria in our DFT calculation were 1 × 10⁻⁶ eV and 0.02 eV Å⁻¹ for electronic and ionic optimization. VASPKIT was used for electronic structure and density of states (DOS) analysis.³⁷ The LOBSTER package was used to analyze the bonding orbital and binding strength between Fe and adsorption via COHP analysis.^38,39 The visualized orbitals were calculated using the wannier90 package.⁴⁰

The adsorption free energy and geometry relaxation at constant electrode potential were computed by VASPsol integrated with the optimization algorithm developed by Xiao et al.^41–44 The implicit solvation based on the linearized Poisson–Boltzmann model (PBM) was adopted to describe the electrochemical response behavior of the aqueous electrolyte with a dielectric constant of 80 and Debye length of 3 Å in our simulation.⁴⁵ The potential-dependent free energy can be expressed as:

E = Evasp − EF × Δq	(1)

where E_vasp is DFT energy determined by VASPsol, E_F is the Fermi energy and Δq is the extra number of electrons. The electrode potential referenced to the standard hydrogen electrode U is:

U = −4.6V − (EF − Esol)	(2)

where E_sol is electrostatic potential in the bulk of solution. Then, the calculated potential-dependent free energies can be fitted using a quadratic function:⁴⁶

E(U) = −1/2 × C × (U − U0)2 + E0	(3)

where C is the interfacial capacitance of the electrode, U₀ is the potential of zero charge (PZC), and E₀ is the energy at the PZC. For each structure, we calculated at charges ranging from 0e to +3e with the increment of +0.5e; corresponding fit parameters are shown in Tables S1–S3.† The energy of the reacting molecules is given in Table S4.† The potential-dependent free energy was used to construct the free energy profile of the ORR and set up the micro-kinetics simulation; all the detailed information is provided in the ESI.†

The ab initio molecular dynamics (AIMD) simulation was performed to validate the stability of Fe spin states in contact with explicit water solvation. Specifically, six H₂O molecules were displaced on each side of FeN^IS₄, and AIMD simulations were carried out using the Nosé-thermostat at 300 K for 1 ps with a time step of 0.5 fs. The out-of-plane displacement of the Fe and spin states were analyzed and compared with the implicit solvent model.

Results and discussion

Spin rearrangement of under external strain

In the FeN₄ model, Fe is coordinated with four N atoms to form a square-planar complex. Four N atoms in sp² hybridization donate six electrons in total and form four σ bonds with the dsp² hybridized Fe center,⁴⁷ as shown in Fig. 1a. Except two electrons from the 4s orbital of Fe that have bonded to N, the five 3d orbitals could be selectively filled with the remaining six 3d electrons, generating three possible spin states: S = 0, 1, or 2. Here we note that, in the case of HS (S = 2), the d_xy orbital filled with spin-up electrons is actually the antibonding orbital σ* of Fe–N bonds, as shown in Fig. 1a. This attribution has been clearly validated by the COHP analysis (Fig. S2†) of the Fe–N bond in FeN₄. We obtained four structures of three spin states, among which S = 2 has two configurations with out-of-plane Fe and in-plane Fe, as marked by FeN^HS₄-1 and FeN^HS₄-2 in Fig. S1 and S3.† Their stability follows the order: FeN^IS₄ (S = 1) > FeN^HS₄-1 (S = 2) > FeN^HS₄-2 (S = 2) > FeN^LS₄ (S = 0). Then, we applied a monoaxial strain from −6% to +8% along the y direction to these four structures as shown in Fig. 1b. The relationship of energy vs. strain in Fig. 1c shows that spin state transition from S = 1 to S = 2 occurs when the applied strain is beyond +6% or below −4%. To understand this strain-induced spin state transition, we compared the structural parameters under different strains (Fig. S4†). The Fe–N bond length increases with strain in FeN^IS₄, FeN^HS₄-1 and FeN^LS₄ due to the expanding N₄ vacancy site (Fig. 1d), which weakens the σ_Fe–N bond and reduces the energy of the antibonding orbital. By comparing PDOS of FeN^IS₄ under different strains, we found that the energy of the d_xy orbital decreases when the strain increases but other d orbitals are less affected by strain (Fig. S5†). Partial spin density was found on N atoms in FeN^HS₄-1 and FeN^HS₄-2 due to the electron filling in the d_xy orbital, which refers to the antibonding orbital , as shown in Fig. S6.† This anomalous electron filing in the antibonding orbital also explains why FeN₄ energetically prefers IS at equilibrium in Fig. 1c. To lower the energy of the d_xy orbital and stabilize the HS state, application of tensile strain is a feasible way. As shown in Fig. 1c and e, under 6% tensile strain, FeN₄ prefers the HS electronic configuration rather than IS. When the strain is between −4% and +6%, the four nonbonding d orbitals (the d_xy orbital forms four σ_Fe–N bonds) are filled with six d electrons. According to Hund's rule, the electrons should occupy as many different orbitals as possible, which lead to the most preferable IS electronic configuration, as shown in Fig. 1e. For FeN^HS₄-1, tensile or compressive strain increases the Fe–N bond length, but large tensile strain could drag the out-of-plane Fe back to the graphene plane, resulting in the transition to FeN^HS₄-2 (Fig. 1e and S4†). When the applied strain is below −2%, due to the bending of the graphene substrate under compressive strain (Fig. 1d), the Fe–N bond in the out-of-plane configuration continues to increase and weakens the σ_Fe–N bond, making the HS state more favorable. We also calculated the strain applied on the y-axis in bare FeN₄. As shown in Fig. S7,† when tensile strain is applied to the y-axis, the same spin state transition characteristics are shown as when strain is applied to the x-axis. Specifically, a spin state transition from S = 1 to S = 2 occurs when the applied strain is beyond +6%. However, under compressive strain, the IS state remains the most stable spin state instead of turning to out-of-plane configuration of the HS state. This can be attributed to the different compressibility of graphene in different directions. As shown in Fig. S8,† under −6% x-axis strain, the compression ratio of the hexagonal ring is only 4.2%, while it increases to 5.5% with the strain along the y axis. This means that the compression ratio of the FeN₄ center under x-axis strain is larger than that of the y-axis strain, resulting in the insensitive response of the spin transition under y-axis strain.

To demonstrate the feasibility of applying strain, we calculated the formation energy of FeN₄ with different spin states under different strains along the x-axis. As shown in Fig. S9,† the formation energy of IS and HS states of FeN₄ with no strain is negative, while the formation energy of LS is slightly positive, indicating good stability of IS and HS FeN₄ under zero strain. Under the tensile strain, the formation energies of all spin states decrease, suggesting that the formation of FeN₄ becomes more favorable in stretched graphene. In contrast to the tensile strain, the formation energy of all spin states increases under compressive strain, except for the HS state with an out-of-plane configuration (FeN^HS₄-1), which can be attributed to the out-of-plane displacement of Fe that weakens the strong interatomic repulsion. These results indicate that it is feasible to tune the spin state of FeN₄ by tensile strain in the experiment from a thermodynamic perspective. Since the tensile strain is experimentally more feasible for graphene-based SACs,^48–50 we will mainly focus on the x-axis tensile strain of 4% and 8% in ORR intermediates.

Strain-induced spin regulation in five-coordinate FeN₄–X (X = *OH/*O/*OOH)

Considering the strain-induced spin state transition in bare FeN₄, the possible spin state transition in ORR intermediates should be taken into consideration. According to early studies, within the operating potential of the ORR, the reaction intermediates can be adsorbed on one side of the Fe–N–C monolayer and form a new active site with five-fold coordination on Fe.^28,51–53 We first constructed three possible active sites in HS/IS states without lattice strain: FeN₄*OH (S = 5/2, 3/2), FeN₄*O (S = 2, 1), and FeN₄*OOH (S = 5/2, 3/2). Their corresponding electronic configuration is determined by PDOS analysis, as shown in Fig. S10.† By comparing the spin density and COHP analysis results of Fe–O/Fe–N bonds, we interestingly found that the Fe center is still in dsp² hybridization, the same as that in the bare FeN₄. This finding is different from the prediction by conventional Valence Shell Electron Pair Repulsion (VSEPR) theory which should predict d²sp³ hybridization in the Fe center based on the number of valence electrons. The electronic configuration of IS/HS FeN₄–X (X: *OH/*O/*OOH) is similar to that of bare FeN₄. Two electrons are either paired in the nonbonding d_x2–y2 orbital for the IS configuration or individually filled in d_xy and d_x2–y2 orbitals for the HS configuration because the adsorbate does not change the bonding type of Fe–N, as shown in Fig. S11.† However, besides the formation of σ_Fe–N, the Fe center also bonds to the absorbed OH and O via d_xz, d_yz, and d_z2 orbitals, leading to the presence of spin density in σ_Fe–O and π_Fe–O, as shown in Fig. 2. The sp² of O and d_z2 of Fe form a σ bond with three electrons, leading to the bond order of 0.5 due to one electron in the orbital.²³ This has been validated by COHP analysis in Fig. 2c, which shows that the –ICOHP value of Fe(d_z2)–O(p_z) for spin-up electrons is close to 0 (Fig. 2c and S12†). Besides, the p_x orbital of O forms a π₂³ bond with the d_xz orbital of Fe with the bond order of 0.5, while the bonding interaction of p_y–d_yz is relatively much weaker, as shown in Fig. 2c and S12.† By scanning the rotational potential energy surface of O–H orientation in Fig. 2d, we found that the OH bond tends to align with the y axis (pointing to a six-membered ring) rather than the x axis (pointing to a five-membered ring) with the energy difference of 0.03 eV. This observation indicated that the coordination of Fe–N₄ is not a rigid D_4h symmetry. An identical bonding type was also found in FeN₄*OOH which includes the σ_Fe–O and π_Fe–O both with the bond order of 0.5 (Fig. S13†) and the O–O bond also tends to align with the y axis (Fig. S14†). Due to the partial filling of electrons in the and orbitals, we found obvious spin density on O atoms (Fig. S15a–d†). In FeN₄*O, the empty O(p_z) orbital bonds to the Fe(d_z2) orbital, generating a σ_Fe–O bond with the bond order of 1. The d_xz and d_yz orbitals in Fe form two π³₂ bonds with the p_x and p_y orbitals of O in side-by-side mode (Fig. 2b and S16†) with obvious spin density observed on O atoms in Fig. S15e–f.† In the HS state of FeN₄*O, partial electron transfer to the antibonding orbital (Fig. S16b†) was observed. The COHP analysis shows that the spin-up –ICHOP for the σ bond in Fe–O decreases from 0.79 (IS) to 0.24 (HS), suggesting that the valence of O is between −2 and −1. This result also agrees with the experimental finding by Frevel et al., who reported the coexistence of μ₂-O and μ₁-O oxyl species in OER intermediates through O K-edge X-ray absorption spectroscopy.⁵⁴ In summary, through the electronic structure and wavefunction analysis, we find that FeN₄ and FeN₄–X possess the same dsp² hybridization mode on Fe, generating four stable σ_Fe–N bonds. The IS–HS transition is dominated by the spin electron transfer from the nonbonding d_x2–y2 orbital to the antibonding (d_xy) orbital. With the adsorption of *OH and *OOH, Fe–O forms one σ_Fe–O bond (d_z2–sp²) and one π³₂ bond (d_xz–p_x) both with the bond order of 0.5. With the *O adsorption, Fe–O forms one σ_Fe–O bond (d_z2–p_z) and two π³₂ bonds (d_xz–p_z and d_yz–p_y) with the bond order of 1 and 0.5 (for each π bond) respectively. Different from FeN₄*OH and FeN₄*OOH, the antibonding (d_z2–p_z) orbital also partially contributes to the spin electron transfer in the IS–HS transition.

As we mentioned before, since the IS–HS transition is mainly dominated by the spin electron transfer from the nonbonding d_x2–y2 orbital to the antibonding (d_xy) orbital, tuning the Fe–N bond length is a feasible way to manipulate the spin state of Fe. The energy and structural parameters of FeN₄–X under different strains are listed in Table 1. Comparing the structural parameters in Table 1, *OH/*O/*OOH adsorption brings an out-of-plane displacement (0.46–0.8 Å) on Fe and increases the Fe–N bond length from 1.90–1.97 Å (FeN₄) to 1.95–2.06 Å (FeN₄–X). This observation explains why the HS of five-coordinate FeN₄–X is energetically more favorable than IS.

---

With the introduction of tensile strain of 4% and 8%, an unexpected spin state transition was observed in FeN₄–X. Under 4% and 8% tensile strain, the spin state of FeN₄*OOH and FeN₄*OH varies from S = 5/2 to S = 2, while the FeN₄*O still prefers HS (S = 2) under strain, as shown in Table 1. Unlike the regular spin state transition between d_xy and d_x2–y2 orbitals, our PDOS analysis in Fig. 3 indicates that the transition from S = 5/2 to S = 2 originates from external spin electron injection in the d_xy orbital, while electronic occupation of other d orbitals is not affected. Bader charge analysis shows that the strain-induced electron filling in the d_xy orbital comes from the neighboring N and C atoms, as shown in Fig. S17 and S18.† In Fig. 3, the peak position of the d_xy (spin-up component) in FeN^IS₄*OOH/*OH is slightly beyond the Fermi level at 0% strain. With the increase of strain, the energy of d_xy gradually drops down below the Fermi level, eventually leading to the electron filling and spin state transition. Besides the charge density analysis, the wavefunction analysis could provide a deeper understanding of this spin state change. We performed the COHP analysis of Fe–N bonds in FeN^IS₄*OH and FeN^HS₄*OOH at different strains, as shown in Fig. S19.† As tensile strain increases, the –ICOHP of the Fe–N bond gradually decreases, indicating that the tensile strain weakens the Fe–N bond. Additionally, the –ICOHP value of the spin-down in the IS state is consistently larger than that of HS until 8% strain. This indicates that the overlapping part of the wavefunction in N and Fe is significantly weakened, leading to the partial bond breaking and backward electron transfer. Different from FeN₄*OH and FeN₄*OOH, the spin state transition was not observed in FeN₄*O due to the relatively shorter Fe–N bond, as shown in Table 1. The presence of two π³₂ bonds (d_xz–p_z and d_yz–p_y) in FeN₄*O could inhibit the out-of-plane shift of Fe, thus shortening the Fe–N bond and preventing the spin state transition.

Since the four- and five-coordinate Fe are both ORR-active, it's important to track how the strain-induced spin state transition affects the partition of different active moieties. To evaluate the potential-dependent evolution of active sites with HS/IS states, we calculated constant-potential free energy of FeN₄–X at different strains, as shown in Fig. S20–S22.† At 0% strain, potential-dependent spin state crossover was observed in the ORR reaction potential window (0–1.23 V vs. SHE) for both bare FeN₄ and FeN₄*O, and the corresponding potential was approximately 1.1 V vs. SHE. This is attributed to the effect of electric double layer capacitance of solution in different spin states. With the increase of strain, the crossover potential increases owing to the larger energy difference at the potential of zero charge. At the same time, FeN₄*OOH also shows spin crossover at 4% and 8% strains. The lower energy of the IS state below the crossover potential can be attributed to electron injection in the d_xy orbital, but increasing potential could lead to electron loss in FeN₄. When the electrode potential reaches the crossover potential, an electron loss in the d_xy orbital is observed, resulting in the spin state turning back to S = 3/2, as illustrated in Fig. S23.† However, the FeN₄*OH under 8% strain exhibits the same energy in HS and IS when U > 0.50 V vs. SHE due to the electron loss in the d_x2–y2 orbital instead of the d_xy orbital (Fig. S24†).

Then, we calculated the coverage of each configuration and spin states at different potentials and different strains based on formation energy, as illustrated in Fig. S25.† Two dominant active sites (FeN^IS₄ and FeN^HS₄*OH) at 0% strain are found in our calculation. As shown in Fig. 4a, the active site is dominated by bare FeN^IS₄ when U < 0.45 V vs. SHE and by FeN^HS₄*OH when U > 0.47 V vs. SHE. As strain increases to 4%, FeN^HS₄-2 becomes the dominant active site with the coverage potential range increasing to U < 0.62 V vs. SHE, while the coverage of FeN^IS₄ at U = 0 V vs. SHE is only 30% and gradually decreases as the potential increases. When U > 0.62 V vs. SHE, the active site is dominated by FeN^HS₄*OH. With further increasing the strain to 8%, the coverage of FeN^IS₄ is very low and can be ignored. The active site is dominated by FeN^HS₄, and the coverage potential range further increases to U < 0.86 V vs. SHE. According to our discussion on spin stability of bare FeN₄ under different strains, FeN^HS₄-2 becomes the main active site in the low potential range at 4% and 8% strain because tensile strain enhanced the stability of the HS state. Furthermore, the coverage range of FeN^HS₄*OH active sites shifts toward a higher potential region with strain, which indicates that tensile strain reduces the adsorption strength between *OH and Fe. We used LOBSTER to analyze the strength of the Fe–O bond in FeN^HS₄*OH under different strains and found that the –ICOHP value decreases as strain increases, suggesting the weakened Fe–O bond as the strain increases (Fig. S26†).

Unusual spin regulation in the six-coordinate X–FeN₄–Y (X/Y = *OH/*O/*OOH)

Subsequently, we calculated all the possible ORR intermediates with a six-coordinate Fe center with the HS/IS state at 0% strain, and the corresponding results are listed in Table 2. Interestingly, the spin state of six-coordinate Fe is either S = 3 or 1 in X–FeN₄–Y (X/Y: *OOH/*OH) and S = 2 or 1 in O*FeN₄–X (X: *OOH/*O/*OH). We calculated the PDOS to identify the orbital occupation of HS/IS intermediates in Fig. S27.† Interestingly, it shows that no matter what ligands are bonded to Fe, the d orbital of Fe always retains at least 4 electrons, which means the highest valence of Fe is +4. In order to reveal the bonding information of six-coordinate ORR intermediates, according to the PDOS and COHP analysis, we provided a schematic of electronic occupation and bonding types in X–FeN₄–Y, as shown in Fig. 5. In the six-coordinate model, the Fe center is bonded with four N atoms and two axial ligands to form an octahedral complex with d²sp³ hybridization.⁴⁷ Four hybridized d²sp³ orbitals of Fe are bonded to the sp² orbital of N in the xy plane and the other two d²sp³ orbitals are bonded to O atoms in the z direction, as shown in Fig. 5a. In X–FeN₄–Y (X/Y: *OH/*OOH), besides the σ_Fe–N and σ_Fe–O formed by d²sp³ hybridized orbitals, the p_x orbital of adsorbates (OH/OOH) and the d_xz orbital form a delocalized π⁵₃ bond, as described in Fig. 5a. Different from the five-coordinate cases, we obtained the S = 1 and S = 3 spin states in the HO*FeN₄*OH complex. In the S = 1 state, the nonbonding d_yz orbital of Fe and the π⁵₃ bond (d_xz orbital) are filled with two unpaired electrons. To convert the S = 1 to S = 3 state, the spin electron transfer from nonbonding d_x2–y2 to and from σ_Fe–O to is observed in Fig. 5a. According to the schematic description of orbital occupations, the bond order of Fe–O is 1 in the S = 1 state, indicating a stronger OH binding in HO*FeN₄*OH than that in FeN₄*OH. This also can be confirmed by the fact that the Fe–O bond length in HO*FeN₄*OH is shorter than that in FeN₄*OH/*OOH (Table S5†). However, in the HS (S = 3) state, PDOS shows a spin electron filling in the d_z2 orbital due to the electron transfer from σ_Fe–O to . According to the COHP analysis of the Fe–O bond in Fig. S28,† the –ICOHP value of the spin up state of Fe(d_z2)–O(p_z) bond shows a significant difference in the HS state (0.44) and IS state (0.72), indicating the typical spin electron transfer in the HS–IS transition. In addition, the spin charge density on O also can be clearly observed in the HS state, as shown in Fig. S29a.† Dumbbell-shaped spin charge density on O in parallel to the y-axis in the IS state is observed in Fig. S29b,† indicating the spin up electron filling in the π⁵₃(O–Fe–O) bond. Bader charge analysis of Fe shows that the HS and IS states hold 6.17|e| and 6.24|e| electronic charges, suggesting that no external charge transfer exists in the HS–IS transition. This conclusion was also supported by the observed increasing Fe–O bond length from 1.85 Å (IS) to 2.00 Å (HS), as shown in Fig. S29.†

---

With the replacement of *OH/*OOH by *O, the σ_Fe–O is formed by the empty O(p_z) orbital and the Fe(d²sp³) orbital. The p_x and p_y orbitals of O could form the π bonds with the d_xz and d_yz orbitals of Fe, as illustrated in Fig. 5b and c. The presence of π_Fe–O bonds can be clearly confirmed by COHP analysis in Fig. S30.† Taking the O*FeN₄*OH as an example, the –ICOHP value of the d_xz–p_z orbital is 0.00 and 0.14 for spin up and spin down states, indicating that the π_Fe–O in Fe–OH is very weak. Compared to the Fe–OH bond, the –ICOHP value of the d_xz–p_x orbital (Fe–O bond) is 0.00 and 0.66 with a relatively higher peak in COHP analysis, indicating that there is a typical π_Fe–O bond formed by the Fe(d_xz) and O(p_x) in Fe–O. Thus, the *OH ligand only provides the p_y orbital to form the π bond with Fe, while the *O ligand could provide both p_x and p_y orbitals to form the π bond. Therefore, there is a π³₂ and a π⁵₃ in the HO*FeN₄*O, but two π⁵₃ in the O*FeN₄*O. This bonding scenario of three-centered delocalized π bonds is similar to that of the HN₃ and N₃⁻ anion. In addition, the COHP analysis in Fig. S28† and Fig. S31† shows that the HS–IS transition doesn't change the electronic occupation of Fe–O bonds in HO*FeN₄*O and O*FeN₄*O. In summary, the six-coordinate X–FeN₄–Y adopts the d²sp³ hybridization in Fe favoring the IS in all intermediates (Table 2) due to the unshifted Fe atom (along the z axis) in the graphene plane. The *OH and *O ligands could provide p_y and p_y(p_x) orbitals to form three-centered π bonds, leading to highly stable O–Fe–O bonding that is not affected by the HS–IS transition in HO*FeN₄*O and O*FeN₄*O.

With the introduction of +4% tensile strain on X–FeN₄–X (X/Y: *OOH/*O/*OH), a strain-induced spin state transition from IS (S = 1) to HS (S = 3) was found in HO*FeN₄*OOH and HO*FeN₄*OH due to the expanded N₄ vacancy as we discussed before, but the IS–HS transition was not observed in HOO*FeN₄*OOH and O*FeN₄–X (X: *OOH/*O/*OH). However, we interestingly found the IS(S = 1)–IS(S = 3/2) transition in HOO*FeN^IS₄*OOH at +4% strain due to an extra electron injection from N atoms, as shown in Fig. S32.† This electron injection-induced S = 1 to S = 3/2 spin state transition was also observed in HO*FeN₄*OOH and HO*FeN₄*OH (Fig. S33†) when the strain increases to +8%. This phenomenon can be explained by the effect of *O/*OH ligands on the energy of the d_xy orbital. Owing to the presence of π, the Fe–O bond in *O-containing adsorbates is stronger than that in *OH-containing adsorbates, leading to the larger orbital splitting energy in the octahedral crystal field and higher energy of the d_xy orbital. Thus, it is energetically more difficult to induce the electron filling of the d_xy orbital in *O-containing intermediates. This analysis also explains why the HOO*FeN₄*O, O*FeN₄*O and HO*FeN₄*O still retain the IS state under +4% strain while HOO*FeN₄*OH and HO*FeN₄*OH converted to HS, as shown in Table 2.

Enhanced ORR kinetics via strain-induced spin regulation

After a detailed discussion on the spin states and electronic configuration, we calculated the energy vs. potential relation of all ORR intermediates under 0%, +4%, and +8% strain, as shown in Fig. S34–S36.† Based on potential-dependent partition of different active sites in Fig. 4, the ORR pathways and corresponding reaction free energy under different strains were obtained as shown in Fig. 6 and S37.† According to the reaction free energy results, the ORR only occurs on the FeN₄ and FeN₄OH sites, and we defined it as path A and path B respectively. At 0% strain, the bare FeN₄ and HO*FeN₄–X (X: *OOH/*O/*OH) prefer IS, while FeN₄–X prefer HS, leading to the ligand-dependent IS/HS switching ORR mechanism (Fig. 6a). Under +4% strain, the bare FeN₄ and HO*FeN₄–Y turn to the HS state, while the HO*FeN₄–Y (Y: *OOH/*OH) goes back to the IS at low overpotential. FeN₄*OOH exhibits the IS state at high overpotential owing to the electron injection in the d_xy orbital and its induced spin state transition, but it turns to HS at low overpotential (Fig. 6c). With the increase of strain to +8%, all the reaction intermediates prefer the HS state at low overpotential (Fig. 6e), leading to the steady HS Fe center along the ORR path. According to potential-dependent reaction free energy in Fig. S37,† we confirmed that whether there is applied strain or not, *OH desorption and *OOH adsorption are the potential-limiting step in path A and path B respectively. In addition, we also found that OH desorption is easier in HO*FeN₄*OH than in FeN₄OH despite the orbital occupation suggesting a higher Fe–O bond order in HO*FeN₄*OH because the orbital is unoccupied. The limiting potential of path A is 0.4 V at 0% strain, 0.7 V at +4% strain and 0.75 V at +8% strain since the strain-weakened Fe–O bond promoted *OH desorption (Fig. S26†), which significantly enhances the ORR activity.⁴⁸ For path B, the limiting potential is 0.8 V at 0% strain, 0.7 V at 4% strain and 1.05 V at 8% strain due to the spin state transition of active sites. According to the polarization current simulation, the overall onset potential is reported as 0.8 V at 0% strain, 0.7 V at 4% strain, and 1.05 V at 8% strain as shown in Fig. 6. These results suggest that the strain-induced spin manipulation is a feasible strategy to improve the ORR performance of Fe–N–C in fuel cell applications.

Finally, to validate the impact of solvation on the Fe spin state in our simulation, we also employed the AIMD simulation with explicit water solvation in contact with FeN₄. Our simulation shows that, within 1 ps, the Fe atoms shift slightly in the vertical direction (less than 0.05 Å), and the spin state of FeN₄ does not change, which indicates that explicit solvation does not significantly affect the position of Fe in the z direction and does not affect the preferable spin state, as shown in Fig. S38.† In addition, we also calculated the free energy of IS and out-of-plane HS of FeN₄ with a H₂O adsorption in the implicit solvent model, as shown in Fig. S39.† After adsorption of H₂O, the IS state is still the most preferable state in FeN₄ and the out-of-plane displacement of Fe is 0.05 Å, which is much smaller than that in FeN^HS₄ (0.53 Å). These results suggest that the spin state of Fe is only dominated by chemical factors, such as valence or ligand change, while the electrostatic interaction from solvation has very little influence.

Summary and conclusion

In our quest to unravel the intricacies of spin multiplicity in ORR intermediates on Fe–N–C, we conducted a systematic exploration of the electronic configuration and bonding paradigms of both high-spin (HS) and intermediate-spin (IS) states across the ORR process. Our comprehensive analysis, based on Projected Density of States (PDOS), Crystal Orbital Hamilton Population (COHP), and spin density visualization, has unveiled the nuanced electronic configurations of HS/IS states in these ORR intermediates. We discovered that within the four-coordinate FeN₄, the nonbonding orbitals d_yz, d_z2, d_xz, and d_x2–y2 are filled with six d electrons, culminating in the most stable IS state according to Hund's rule. A pivotal finding of our work is the dominance of spin electron transfer from the nonbonding d_x2–y2 orbital to the antibonding orbital in governing the IS–HS transition, a process that can be modulated by adjusting the Fe–N bond length. The application of tensile strain has been shown to drive an IS–HS transition, attributed to the decrease of the orbital energy with concurrent elongation of the Fe–N bond. In five-coordinate FeN₄–X species (X = *OH/*OOH), our findings reveal the complex electronic configuration of σ_Fe–O and π_Fe–O bonds, with implications for understanding the IS–HS transition under strain. Notably, the *O adsorption introduces two π bonds between Fe–O formed by Fe(d_xz)–O(p_x) and Fe(d_yz)–O(p_y), leading to the steady bonding of Fe and O with the bond order of 2. In six-coordinate X–FeN₄–Y complexes, we observed a d²sp³ hybridization pattern in Fe, with the highest oxidation state of +4, and the formation of both σ and π bonds in *OH/*OOH and *O adsorbates. Intriguingly, the bonding mode in HO–Fe–O and O–Fe–O complexes parallels that of hydrazoic acid and the azide anion, featuring three-centered π bonds. The tensile strain effect extends to these six-coordinate X–FeN₄–Y complexes, prompting a spin state transition from S = 1 to S = 3/2 due to electron transfer from N/C to Fe. At +8% tensile strain, the spin states of most ORR intermediates turn to S = 3 or S = 2, underscoring the critical role of spin manipulation in tuning catalytic activity. Our potential-dependent formation energy calculations reveal a dynamic conversion between FeN^IS/HS₄ and FeN^HS₄*OH as active sites across different potentials, with FeN^HS₄ gaining dominance under increased strain. By pinpointing the most stable spin states for each intermediate, we demonstrate that tensile strain effectively suppresses ligand-related spin transitions, stabilizing the HS state on the Fe center across the whole ORR process. Finally, the simulated current density indicates that +8% strain significantly improves the ORR activity, pushing the onset potential to 1.05 V vs. SHE. Our study not only elucidates the origins of spin multiplicity in Fe–N–C but also proposes a viable strategy for enhancing ORR activity through spin manipulation. Our work could significantly facilitate the development of spin-related electrocatalysis from both experimental and theoretical perspectives.

Data availability

The authors confirm that the data supporting the findings of this study are available within the article and its ESI.† The code and input parameters of this simulation work are available on request from the authors.

Author contributions

C. Z. designed the project. M. Y. performed the simulation, data analysis, and manuscript writing. M. Y., Y. D., J. W., Y. C., A. L., E. K. and C. Z. participated in the technical discussion and manuscript revision.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was funded by the Funding of NJUST (No. TSXK2022D002), Natural Science Foundation of Jiangsu Province (No. BK20220929), Natural Science Foundation of China (No. 22303039), and Startup Grant of NJUST.

References

1. B. C. Steele and A. Heinzel, Nature, 2001, 414, 345–352 CrossRef CAS PubMed .
2. M. Shao, Q. Chang, J.-P. Dodelet and R. Chenitz, Chem. Rev., 2016, 116, 3594–3657 CrossRef CAS PubMed .
3. M. K. Debe, Nature, 2012, 486, 43–51 CrossRef CAS PubMed .
4. J. Shui, M. Wang, F. Du and L. Dai, Sci. Adv., 2015, 1, e1400129 CrossRef PubMed .
5. Y. Wang, H. Su, Y. He, L. Li, S. Zhu, H. Shen, P. Xie, X. Fu, G. Zhou, C. Feng, D. Zhao, F. Xiao, X. Zhu, Y. Zeng, M. Shao, S. Chen, G. Wu, J. Zeng and C. Wang, Chem. Rev., 2020, 120, 12217–12314 CrossRef CAS PubMed .
6. X. Cui, W. Li, P. Ryabchuk, K. Junge and M. Beller, Nat. Catal., 2018, 1, 385–397 CrossRef CAS .
7. W. H. Lee, Y.-J. Ko, J.-Y. Kim, B. K. Min, Y. J. Hwang and H.-S. Oh, Chem. Commun., 2020, 56, 12687–12697 RSC .
8. H. Fei, J. Dong, D. Chen, T. Hu, X. Duan, I. Shakir, Y. Huang and X. Duan, Chem. Soc. Rev., 2019, 48, 5207–5241 RSC .
9. Y. Wang, J. Mao, X. Meng, L. Yu, D. Deng and X. Bao, Chem. Rev., 2018, 119, 1806–1854 CrossRef PubMed .
10. M. B. Gawande, P. Fornasiero and R. Zbořil, ACS Catal., 2020, 10, 2231–2259 CrossRef CAS .
11. Y. Dai, F. Kong, X. Tai, Y. Zhang, B. Liu, J. Cai, X. Gong, Y. Xia, P. Guo and B. Liu, Electrochem. Energy Rev., 2022, 5, 22 CrossRef CAS .
12. W. Song, C. Xiao, J. Ding, Z. Huang, X. Yang, T. Zhang, D. Mitlin and W. Hu, Adv. Mater., 2024, 36, 2301477 CrossRef CAS PubMed .
13. S. Liu, C. Li, M. J. Zachman, Y. Zeng, H. Yu, B. Li, M. Wang, J. Braaten, J. Liu, H. M. Meyer III, M. Lucero, A. J. Kropf, E. E. Alp, Q. Gong, Q. Shi, Z. Feng, H. Xu, G. Wang, D. J. Myers, J. Xie, D. A. Cullen, S. Litser and G. Wu, Nat. Energy, 2022, 7, 652–663 CrossRef CAS .
14. X. Zhang, L. Truong-Phuoc, T. Asset, S. Pronkin and C. Pham-Huu, ACS Catal., 2022, 12, 13853–13875 CrossRef CAS .
15. H. Shen, T. Thomas, S. A. Rasaki, A. Saad, C. Hu, J. Wang and M. Yang, Electrochem. Energy Rev., 2019, 2, 252–276 CrossRef CAS .
16. K. Song, B. Yang, X. Zou, W. Zhang and W. T. Zheng, Energy Environ. Sci., 2023, 17, 27–48 RSC .
17. J. Li, J. Ma, Z. Ma, E. Zhao, K. Du, J. Guo and T. Ling, Adv. Energy Sustainability Res., 2021, 2, 2100034 CrossRef CAS .
18. Z. Zhang, P. Ma, L. Luo, X. Ding, S. Zhou and J. Zeng, Angew. Chem., Int. Ed., 2023, 62, e202216837 CrossRef CAS PubMed .
19. L. Scarpetta-Pizo, R. Venegas, P. Barrías, K. Muñoz-Becerra, N. Vilches-Labbé, F. Mura, A. M. Méndez-Torres, R. Ramírez-Tagle, A. Toro-Labbé, S. Hevia, J. H. Zagal, R. Oñate, A. Aspée and I. Ponce, Angew. Chem., Int. Ed., 2024, 136, e202315146 CrossRef .
20. Y. Sun, S. Sun, H. Yang, S. Xi, J. Gracia and Z. J. Xu, Adv. Mater., 2020, 32, 2003297 CrossRef CAS PubMed .
21. D. Xia, X. Yang, L. Xie, Y. Wei, W. Jiang, M. Dou, X. Li, J. Li, L. Gan and F. Kang, Adv. Funct. Mater., 2019, 29, 1906174 CrossRef CAS .
22. G. Yang, J. Zhu, P. Yuan, Y. Hu, G. Qu, B.-A. Lu, X. Xue, H. Yin, W. Cheng, J. Cheng, W. Xu, J. Li, J. Hu, S. Wu and J.-N. Zhang, Nat. Commun., 2021, 12, 1734 CrossRef CAS PubMed .
23. H. Zhang, H.-C. Chen, S. Feizpoor, L. Li, X. Zhang, X. Xu, Z. Zhuang, Z. Li, W. Hu, R. Snyders, D. Wang and C. Wang, Adv. Mater., 2024, 36, 2400523 CrossRef CAS PubMed .
24. C. Chen, Y. Wu, X. Li, Y. Ye, Z. Li, Y. Zhou, J. Chen, M. Yang, F. Xie, Y. Jin, C. Jones, N. Wang, H. Meng and S. Chen, Appl. Catal., B, 2024, 342, 123407 CrossRef CAS .
25. J. Hou, Y. Jian, C. Chen, D. Zhang, F. Xie, J. Chen, Y. Jin, N. Wang, X. Yu and H. Meng, Mater. Chem. Front., 2024, 8, 2592–2598 RSC .
26. W. Zhong, Y. Qiu, H. Shen, X. Wang, J. Yuan, C. Jia, S. Bi and J. Jiang, J. Am. Chem. Soc., 2021, 143, 4405–4413 CrossRef CAS PubMed .
27. X. Hu and N. Q. Su, J. Phys. Chem. Lett., 2023, 14, 9872–9882 CrossRef CAS PubMed .
28. M. Yu, A. Li, E. Kan and C. Zhan, ACS Catal., 2024, 14, 6816–6826 CrossRef CAS .
29. D. Benito-Garagorri, I. Lagoja, L. F. Veiros and K. A. Kirchner, Dalton Trans., 2011, 40, 4778–4792 RSC .
30. X. Xin, C. Guo, R. Pang, X. Shi and Y. Zhao, Appl. Surf. Sci., 2021, 570, 151126 CrossRef CAS .
31. R. Wang, L. Zhang, J. Shan, Y. Yang, J. F. Lee, T. Y. Chen, J. Mao, Y. Zhao, L. Yang, Z. Hu and L. Tao, Adv. Sci., 2022, 9, 2203917 CrossRef CAS PubMed .
32. P. E. Blöchl, Phys. Rev. B, 1994, 50, 17953 CrossRef PubMed .
33. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50 CrossRef CAS .
34. G. Kresse and D. Joubert, Phys. Rev. B, 1999, 59, 1758 CrossRef CAS .
35. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865 CrossRef CAS PubMed .
36. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799 CrossRef CAS PubMed .
37. V. Wang, N. Xu, J.-C. Liu, G. Tang and W.-T. Geng, Comput. Phys. Commun., 2021, 267, 108033 CrossRef CAS .
38. V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Phys. Chem. A, 2011, 115, 5461–5466 CrossRef CAS PubMed .
39. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37, 1030–1035 CrossRef CAS PubMed .
40. G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi and J. R. Yates, J. Phys.: Condens. Matter, 2020, 32, 165902 Search PubMed .
41. K. Mathew, V. Kolluru, S. Mula, S. N. Steinmann and R. G. Hennig, J. Chem. Phys., 2019, 151, 234101 CrossRef PubMed .
42. K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. Arias and R. G. Hennig, J. Chem. Phys., 2014, 140, 084106 CrossRef PubMed .
43. J. A. Gauthier, S. Ringe, C. F. Dickens, A. J. Garza, A. T. Bell, M. Head-Gordon, J. K. Nørskov and K. Chan, ACS Catal., 2018, 9, 920–931 CrossRef .
44. Z. Duan and P. Xiao, J. Phys. Chem. C, 2021, 125, 15243–15250 CrossRef CAS .
45. A. Hagopian, M.-L. Doublet, J.-S. Filhol and T. Binninger, J. Chem. Theory Comput., 2022, 18, 1883–1893 CrossRef CAS PubMed .
46. Y. Huang, R. J. Nielsen and W. A. Goddard III, J. Am. Chem. Soc., 2018, 140, 16773–16782 CrossRef CAS PubMed .
47. Y. Cui, C. Ren, Q. Li, C. Ling and J. Wang, J. Am. Chem. Soc., 2024, 15640–15647 CrossRef CAS PubMed .
48. S. Zhu, L. Ding, X. Zhang, K. Wang, X. Wang, F. Yang and G. Han, Angew. Chem., Int. Ed., 2023, 62, e202309545 CrossRef CAS PubMed .
49. L. Li, J. Zhu, F. Kong, Y. Wang, C. Kang, M. Xu, C. Du and G. Yin, Matter, 2024, 7, 1517–1532 CrossRef CAS .
50. Z. Miao, X. Wang, Z. Zhao, W. Zuo, S. Chen, Z. Li, Y. He, J. Liang, F. Ma and H. L. Wang, Adv. Mater., 2021, 33, 2006613 CrossRef CAS PubMed .
51. Z. Huang and Q. Tang, J. Phys. Chem. C, 2022, 126, 21606–21615 CrossRef CAS .
52. F. Sun, F. Li and Q. Tang, J. Phys. Chem. C, 2022, 126, 13168–13181 CrossRef CAS .
53. Y. Wang, Y.-J. Tang and K. Zhou, J. Am. Chem. Soc., 2019, 141, 14115–14119 CrossRef CAS PubMed .
54. L. J. Frevel, R. Mom, J.-J. s. Velasco-Vélez, M. Plodinec, A. Knop-Gericke, R. Schlögl and T. E. Jones, J. Phys. Chem. C, 2019, 123, 9146–9152 CrossRef CAS .

---

Footnote

† Electronic supplementary information (ESI) available: Structural parameters; PDOS and electron occupation; COHP analysis and spin density; chemical reactions and free energy calculation; constant-potential free energy diagrams; occupation calculation. See DOI: https://doi.org/10.1039/d4ta04458d

This journal is © The Royal Society of Chemistry 2024