Atomistic Modeling of Spinel Oxide Particle Shapes and Reshaping under OER Conditions

: The surface configurations of the low-index facets of a set of spinel oxides are investigated using DFT+U calculations to derive surface energies and predict equilibrium nanoparticle shapes via the Wulff construction. Two very different conditions are investigated, corresponding to application either in heterogeneous catalysis or in electrocatalysis. First, the bare stoichiometric surfaces of NiFe 2 O 4 , CoFe 2 O 4 , NiCo 2 O 4 , and ZnCo 2 O 4 spinels are studied to model their use as high-temperature oxidation catalysts. Second, focusing attention on the electrochemical oxygen evolution reaction (OER) and on the CoFe 2 O 4 inverse spinel as the most promising OER catalyst, we generate surface configurations by adsorbing OER intermediates and, in an innovative study, we recalculate surface energies taking into account adsorption and environmental conditions, i.e., applied electrode potential and O 2 pressure. We predict that under OER operating conditions, (111) facets are dominant in CoFe 2 O 4 nanoparticle shapes, in fair agreement with microscopy measurements. Importantly, in the OER case, we predict a strong dependence of nanoparticle shape upon O 2 pressure. Increasing O 2 pressure increases the size of the higher-index (111) and (110) facets at the expense of the (001) more catalytically active facet, whereas the opposite occurs at low O 2 pressure. These predictions should be experimentally verifiable and help define the optimal OER operative conditions.


Introduction
Nanomaterials, such as metal or metal oxide nanoparticle systems, find applications in many areas of technology, such as biomedical, sensors, optics, storage, and energy harvesting materials as well as in heterogeneous catalysis, where they play a dominant role [1].Their function is closely linked to their structure in terms of nanoparticle size and shape [2].Therefore, controlling the shape of metal or metal oxide nanoparticles is one of the major challenges in applications as this has a dramatic effect on material properties, e.g., determining the number and type and performance of active sites in catalysis.
The shapes of nanoparticles can be determined via the Wulff construction, proposed by the mineralogist George Wulff at the beginning of the twentieth century [3].This method was used in geology, mineralogy, and crystallography and was then rediscovered in materials science to predict the shapes of nanoparticles.The Wulff construction allows one to predict the nanoparticle shape through knowledge of their surface energies, data which can be and have been calculated for many systems via, e.g., theoretical simulations based on first-principles quantum mechanical methods [1].Some characteristic equilibrium shapes often found in materials with cubic symmetry are shown in Figure 1.Following the fundamental work by Gibbs [4] and introducing the principle of minimum surface free energy to determine the shape of a crystal at equilibrium, Wulff proposed [3] that the shape of a crystalline material be constructed via the following recipe: where γhkl is the energy required to create a surface of unit area normal to the (hkl) vector and rhkl is the distance from the center of the crystal to the (hkl) plane.This process is repeated for all sets of Miller indexes, hkl, uniquely identifying a surface.The ratio of the surface energy, γhkl, to the distance from the center of the crystal, rhkl, to any plane is constant in the Wulff construction.
In heterogeneous catalysis by metal nanoparticles, the catalytic functionality is strongly linked to their shape, especially for structure-sensitive reactions, which makes predictions of their morphology via, e.g., computational methods, an important goal.For example, Nørskov et al. obtained the equilibrium polyhedron shape of Ru nanoparticles, estimated the catalytic rates by summing contributions of the various facets, and compared with experimentally measurements of transmission electron microscopy (TEM) and catalytic rates, finding reasonable agreement with experiment [5].Vilé et al. investigated supported Ag catalyst nanoparticles (2-20 nm in diameter) to determine active sites for the olefin reaction using the atomistic Wulff construction and reported that the nanoparticles have the highest density of B5 sites, which are the active centers in the reaction [6].Larger nanoparticles may expose higher-index facets: Wulff shapes of large Au nanoparticles were predicted to include (332), (001), (111), (211), and (322) facets, whereas nanoparticles smaller than 16 nm in diameter expose only (111) and (001) facets [7].However, the equilibrium shape rarely contains facets with indexes higher than (001), (111), and (110): high-index surfaces have a higher surface tension and, even if a highindex surface and a low-index surface have equal surface tensions, the low-index face will have a larger area as the high-index face will be steeper and be hidden in the Wulff construction [1].
Bare metal oxide facets have also been considered.In a systematic study, the surface structure, energies, and morphology of bare cubic cobalt spinel nanocrystals with the formula M[CoM′]O4 (M = Mg, Zn, Fe, and Co and M′ = Ni, Al, Mn, and Co) were investigated via density-functional theory plus Hubbard correction (DFT+U) calculations and compared directly with TEM and STEM experiments [8].The authors predicted that the abundance of the (110) faces is always low (below 9%) and the surface energies typically increase in the sequence γ100 < γ111 < γ110 for cobalt spinel nanocrystals.
Computational studies of bare metal or metal oxide nanoparticles are appropriate for modeling high-temperature and low-coverage conditions.At high temperatures and low pressures, there is an entropic driving force for adsorbate species to leave the surface for the gas phase, giving rise to a low-coverage régime.These are the conditions of thermal catalysis, such as those occurring in oxidation processes: chemical looping, water gas shift reaction [9], methanol [10] and propanol [11] oxidation, H2 production [12], etc.Following the fundamental work by Gibbs [4] and introducing the principle of minimum surface free energy to determine the shape of a crystal at equilibrium, Wulff proposed [3] that the shape of a crystalline material be constructed via the following recipe: where γ hkl is the energy required to create a surface of unit area normal to the (hkl) vector and r hkl is the distance from the center of the crystal to the (hkl) plane.This process is repeated for all sets of Miller indexes, hkl, uniquely identifying a surface.The ratio of the surface energy, γ hkl , to the distance from the center of the crystal, r hkl , to any plane is constant in the Wulff construction.
In heterogeneous catalysis by metal nanoparticles, the catalytic functionality is strongly linked to their shape, especially for structure-sensitive reactions, which makes predictions of their morphology via, e.g., computational methods, an important goal.For example, Nørskov et al. obtained the equilibrium polyhedron shape of Ru nanoparticles, estimated the catalytic rates by summing contributions of the various facets, and compared with experimentally measurements of transmission electron microscopy (TEM) and catalytic rates, finding reasonable agreement with experiment [5].Vilé et al. investigated supported Ag catalyst nanoparticles (2-20 nm in diameter) to determine active sites for the olefin reaction using the atomistic Wulff construction and reported that the nanoparticles have the highest density of B5 sites, which are the active centers in the reaction [6].Larger nanoparticles may expose higher-index facets: Wulff shapes of large Au nanoparticles were predicted to include (332), (001), (111), (211), and (322) facets, whereas nanoparticles smaller than 16 nm in diameter expose only (111) and (001) facets [7].However, the equilibrium shape rarely contains facets with indexes higher than (001), (111), and (110): high-index surfaces have a higher surface tension and, even if a high-index surface and a low-index surface have equal surface tensions, the low-index face will have a larger area as the high-index face will be steeper and be hidden in the Wulff construction [1].
Bare metal oxide facets have also been considered.In a systematic study, the surface structure, energies, and morphology of bare cubic cobalt spinel nanocrystals with the formula M[CoM ′ ]O 4 (M = Mg, Zn, Fe, and Co and M ′ = Ni, Al, Mn, and Co) were investigated via density-functional theory plus Hubbard correction (DFT+U) calculations and compared directly with TEM and STEM experiments [8].The authors predicted that the abundance of the (110) faces is always low (below 9%) and the surface energies typically increase in the sequence γ100 < γ111 < γ110 for cobalt spinel nanocrystals.
Computational studies of bare metal or metal oxide nanoparticles are appropriate for modeling high-temperature and low-coverage conditions.At high temperatures and low pressures, there is an entropic driving force for adsorbate species to leave the surface for the gas phase, giving rise to a low-coverage régime.These are the conditions of thermal catalysis, such as those occurring in oxidation processes: chemical looping, water gas shift reaction [9], methanol [10] and propanol [11] oxidation, H 2 production [12], etc.
For nanoparticles under lower temperatures and therefore higher-coverage conditions, instead, the influence of the chemical environment should be included in the Wulff con-struction by rescaling surface energies with adsorption energies [13].For example, Au nanoparticles exposed to CO were found to be more spherical and more reactive compared to Au nanoparticles in non-interacting environments, in agreement with experimental data [7] and also in agreement with previous studies investigating the shape and chemical ordering upon CO and H adsorption of mono-and bi-metallic explicit atomistic nanoparticle models, thus going beyond the ideal Wulff construction [14].Similarly, the surface free energies of MgO low-index crystal facets were investigated as a function of the temperature and water pressure through density functional theory, also comparing predictions with experimental observations [15].It was predicted that (110) is a transient facet and that the equilibrium Wulff shape of MgO crystals only involves (001) and (111) facets.Additionally, the experiment and theory agreed on the pattern of the hydroxylated surface energies as γ111 < γ100 < γ110, supporting the partial dissociation of water on MgO(001).
At the next level of complexity, the important question is how Wulff shape changes under reaction conditions and how to account for this phenomenon at the theoretical level [16,17].There are a few examples of this research in thermal catalysis.In a pioneering study, the size-and corresponding morphology-dependent theoretical activity and selectivity of silver nanoparticles between 9 and 23 nm in diameter for the partial oxidation of propylene and the evolution of these quantities along the progress of the reaction were experimentally measured and compared with theoretical predictions based on the Wulff construction and explicit reaction mechanisms, finding an excellent agreement with experimental observations of evolution of the aspect ratio [16].More recently, the concept of dynamical Wulff construction was explicitly introduced to account for particle reshaping under reaction conditions and employed to predict the shape of Si-doped iron nanoparticles during ammonia synthesis [17].
All this previous work refers to thermal catalysis.In contrast, the application of Wulffbased theoretical protocols in the field of electro-catalysis is, to the best of our knowledge, lacking, despite the ever growing interest of the scientific and technological communities in this field [18,19].In particular, the oxygen evolution reaction (OER) in anion exchange membrane water electrolysis (AEMWE or water splitting) represents a critical step in the realization of a sustainable hydrogen production whence a global hydrogen economy would arise, so that discovering and optimizing active and non-critical-material OER catalysts is one of the major challenges to XXI century chemistry.
Here, we address and fill this lack of knowledge in electro-catalysis by investigating time spinel-type oxides via theoretical atomistic Wulff construction in a wide range of OER electro-catalytic conditions for the first.Notably, spinel-like iron-based spinels are promising OER electro-catalysts [20] due to their significant electrical conductivity, structural stability, and catalytic performance, stemming from multiple valence of the cations and the ability to switch among different oxidation states [21,22].We predict, via first-principles DFT+U, the surface energies of the most abundant (001), (110), and (111) facets of selected spinel oxide structures.To provide an exhaustive ensemble of data, we derive the equilibrium atomistic Wulff shapes of the corresponding nanoparticles under two very different sets of conditions: and ZnCo 2 O 4 (inverse or normal) spinels as bare stoichiometric but non-symmetric surfaces to be compared with previous studies [8] and with characterization experiments and catalytic investigations under high-temperature low-coverage conditions, such as those occurring in heterogeneous thermal catalysis [20]; (ii) The CoFe 2 O 4 inverse spinel, the most promising OER catalyst [23], for which we consider an extensive set of surface configurations with various coverage of adsorbate intermediates relevant during the electrochemical OER (H, OH, H 2 O, O, and O 2 ).
For the OER case, we also explore a wide set of environmental variables (applied electrode potential and O 2 pressure) so as to determine Wulff reshaping under different OER conditions and then finally compare our predictions with available experimental characterization.
In case (i), in agreement with the experiment and a previous systematic study [8], we predict that underhigh-temperature low-coverage conditions (100) typically dominate.
In case (ii), also in agreement with the experiment, we find that the most abundant CoFe 2 O 4 facets change from (100) for bare particles to (111) under OER conditions.Moreover, importantly, in the CoFe 2 O 4 OER case, we find a strong phenomenon of nanoparticle reshaping upon O 2 pressure: increasing O 2 pressure reshapes the catalyst, increasing the size of the higher-index facets, i.e., (111) and (110), at the expense of the (001) facet, whereas the opposite occurs at low O 2 pressure.Given that the OER catalytic activity of CoFe 2 O 4 spinels depends upon the exposure of different facets [24][25][26], these predictions should be experimentally verifiable and help defining the optimal operative conditions for OER catalysis (and also oxidative catalysis), with potential implications, e.g., on the design of OER electro-catalysts fulfilling societal requirement on activity, selectivity, efficiency, and cost.

Materials and Methods
Spin-polarized density functional theory plus Hubbard correction (DFT+U [27]) calculations were performed using a plane wave and ultrasoft pseudopotential framework [28] as implemented in the Quantum ESPRESSO [29] suite of codes.The Perdew-Burke-Ernzerhof (PBE) [30] exchange-correlation function was used as DFT, augmented with Hubbard U parameters chosen as 3.3, 5.5, 3.0, and 2.0 eV for Fe, Ni, Co, and Zn, respectively.It should be noted that we use an U eff value of 3.0 eV for all types of Co atoms (i.e., octahedrally or/and tetrahedrally coordinated, depending on spinel type) to be consistent with a previous systematic theoretical study by Selloni et al. [31][32][33].In our previous study [23], we used a value for U eff of 4.5 eV as optimal for describing iron-based inverse spinels, in which cobalt only occupies an octahedral coordination.However, we verified that there is no significant difference between the two U eff values in terms of configuration energetics.Kinetic energy cutoffs of 40 and 200 Ry were chosen for describing the wave function and the charge density, respectively.A 3 × 3 × 1 Monkhorst−Pack k-point mesh was utilized for energy and structural calculations.The WulffPack-1.2 package was used to predict theoretical equilibrium shapes [34].
Spinels have the general chemical formula AB 2 X 4 where A II is a divalent cation like Mg, Cr, Mn, Fe, Co, Ni, Cu, Zn, Cd, or Sn; B III is a trivalent cation like Al, Ga, In, Ti, V, Cr, Mn, Fe, Fe, Co, or Ni; and X is O, S, Se, etc.A normal spinel can be represented as (A II ) tet (B III ) 2 oct O 4 , whereas an inverse spinel can be represented as (B III ) tet (A II ) oct (B III ) oct O 4 , where the superscript "tet" and "oct" refer to tetrahedral and octahedral sites.The selected catalysts with an inverse spinel structure (NiFe 2 O 4 , CoFe 2 O 4 , and NiCo 2 O 4 ) are ferrimagnetic [34], i.e., they exhibit a spin arrangement: ↑↓↑ for A(Oh)/B(Td)/B(Oh), using the notation (A = Ni, Co; B = Fe, Co), spin up (↑), spin down (↓), and Oh and Td refer to octahedral and tetrahedral sites, respectively.At variance, ZnCo 2 O 4 only exists in the normal spinel structure and is ferromagnetic [35,36], with spin arrangement ↑↑ for A(Td)/B(Oh), where (A = Zn; B = Co).We verified that our DFT+U predictions are consistent with these expectations.
For the surface energy calculations under vacuum, corresponding to low-coverage thermal catalysis, the stoichiometric/non-symmetric slabs utilized for the three considered facets [(001), (110), and (111)], including dipole correction, for inverse and normal spinel slab models are shown in Figures 2 and 3.For the (001) orientation, our model corresponds to a total of 56 atoms per unit cell for a thickness of 9 layers.For the (110) orientation, our slab with a thickness of 9 layers contains 126 atoms per unit cell, while for the (111) orientation our 9-layer thick unit cell contains 112 atoms.Notably, there are two types of possible terminations for the (110) spinel orientation, which are usually labeled A and B in the literature [8,37].As the 110-A termination was reported to be more stable than the 110-B one in previous work [38], the latter was excluded from our analysis.Therefore, the (110) label only refers to the most stable (110)-A termination in this study.
We use Equation (2) to calculate surface energies from optimized energies of bare slab models: The lattice parameters a, b, and c of each stoichiometric slab have been derived after implementing variable-cell (vc) relaxation on the bulk phases of the investigated spinels and are shown in Table S1 of the Supporting Information (SI).For the (001) orientation, our model corresponds to a total of 56 atoms per unit cell for a thickness of 9 layers.For the (110) orientation, our slab with a thickness of 9 layers contains 126 atoms per unit cell, while for the (111) orientation our 9-layer thick unit cell contains 112 atoms.Notably, there are two types of possible terminations for the (110) spinel orientation, which are usually labeled A and B in the literature [8,37].As the 110-A termination was reported to be more stable than the 110-B one in previous work [38], the latter was excluded from our analysis.Therefore, the (110) label only refers to the most stable (110)-A termination in this study.
We use Equation ( 2) to calculate surface energies from optimized energies of bare slab models: The lattice parameters a, b, and c of each stoichiometric slab have been derived after implementing variable-cell (vc) relaxation on the bulk phases of the investigated spinels and are shown in Table S1 of the Supporting Information (SI).For the (001) orientation, our model corresponds to a total of 56 atoms per unit cell for a thickness of 9 layers.For the (110) orientation, our slab with a thickness of 9 layers contains 126 atoms per unit cell, while for the (111) orientation our 9-layer thick unit cell contains 112 atoms.Notably, there are two types of possible terminations for the (110) spinel orientation, which are usually labeled A and B in the literature [8,37].As the 110-A termination was reported to be more stable than the 110-B one in previous work [38], the latter was excluded from our analysis.Therefore, the (110) label only refers to the most stable (110)-A termination in this study.
We use Equation ( 2) to calculate surface energies from optimized energies of bare slab models: The lattice parameters a, b, and c of each stoichiometric slab have been derived after implementing variable-cell (vc) relaxation on the bulk phases of the investigated spinels and are shown in Table S1 of the Supporting Information (SI).
Physchem 2024, 4, FOR PEER REVIEW 6 We also use the bulk optimized lattice parameters to build non-stoichiometric/symmetric inverse spinel models for the high-coverage, i.e., electrochemical conditions, to investigate surfaces with various coverage of adsorbed OER intermediates on CoFe2O4.Non-stoichiometric but symmetric slab unit cells with dimensions of 5.99 Å × 5.99 Å × 25 Å for (001), 5.99 Å × 8.48 Å × 25 Å for (110), and 5.99 Å × 5.19 Å × 25 Å for (111) of inverse spinel slabs were built as shown in Figure 4.The total number of atoms per unit cell for (001), (110), and (111) slabs used in bare and adsorbed surface calculations is summarized in Table S2 of the Supplementary Materials; brief details about the comparison of those slabs are shown in Figure S1.
To calculate surface energies for surfaces with adsorbates, we re-scale them as shown in Equation ( 3), where adsorption energies are taken from optimization calculations of non-stoichiometric slab models [13,39]: In the surfaces covered with oxygen, hydrogen, and water, it is necessary to include entropy terms in the chemical potentials of these species, to obtain the system free energy (G).Free energies are calculated using the same scheme utilized in previous studies [40,41] and explained in detail in our previous work [23].The total energy of each configuration was obtained directly via DFT.The contributions (∆H, ∆ZPE, and T∆S) to the free energies (∆G) of small molecules (H2O, hydrogen, and oxygen) were added empirically as ΔGi = ΔEi + ΔHi + ΔZP-i − TΔSi [41], where standard thermodynamic data [42] are used to obtain the T and P contributions to the G values of aqueous H2O and gaseous H2.
The standard hydrogen electrode electrochemistry model [43] was finally used to calculate how the configuration energies depend upon the applied bias U via the protonelectron (H + + e − ) transfer steps.In this approach, the standard hydrogen electrode (SHE) is used as a reference, so the proton (G[H + ]) and electron (G[e − ]) free energies are replaced by ½G[H2]-|e|U, where G[H2] is the free energy of H2 and U is the electrode potential versus SHE.To get rid of the dependence of the applied potential upon pH, a reversible hydrogen electrode (RHE) reference is employed, so that the equilibrium potential re- To calculate surface energies for surfaces with adsorbates, we re-scale them as shown in Equation (3), where adsorption energies are taken from optimization calculations of non-stoichiometric slab models [13,39]: In the surfaces covered with oxygen, hydrogen, and water, it is necessary to include entropy terms in the chemical potentials of these species, to obtain the system free energy (G).Free energies are calculated using the same scheme utilized in previous studies [40,41] and explained in detail in our previous work [23].The total energy of each configuration was obtained directly via DFT.The contributions (∆H, ∆ZPE, and T∆S) to the free energies (∆G) of small molecules (H 2 O, hydrogen, and oxygen) were added empirically as [41], where standard thermodynamic data [42] are used to obtain the T and P contributions to the G values of aqueous H 2 O and gaseous H 2 .
The standard hydrogen electrode electrochemistry model [43] was finally used to calculate how the configuration energies depend upon the applied bias U via the protonelectron (H + + e − ) transfer steps.In this approach, the standard hydrogen electrode (SHE) is used as a reference, so the proton (G[H + ]) and electron (G[e − ]) free energies are replaced by ½G[H 2 ]-|e|U, where G[H 2 ] is the free energy of H 2 and U is the electrode potential versus SHE.To get rid of the dependence of the applied potential upon pH, a reversible hydrogen electrode (RHE) reference is employed, so that the equilibrium potential required for OER is always 1.23 V vs. RHE at all pH values [44].In other words, in calculating the over-potential, the pH term is canceled by taking the difference between predicted onset potentials and the equilibrium OER potential [45].

Results and Discussion
As anticipated in the Introduction, we divide our study into two parts.First, we focus on NiFe 2 O 4 , CoFe 2 O 4 , NiCo 2 O 4 , and ZnCo 2 O 4 systems as bare stoichiometric surfaces: this corresponds to experiments at high temperature and low coverage, i.e., under the conditions of thermal heterogeneous catalysis [9][10][11][12].This will also allow us to compare our predictions with and validate against previous studies [8] (Section 3.1).Then, we switch to the most innovative part of our investigation, i.e., the Wulff construction under electrochemical OER conditions; in Section 3.2, we consider the CoFe 2 O 4 system and study configurations with various coverage of (H, OH, H 2 O, O, and O 2 ) adsorbate species under electrochemical OER and varying environmental variables (applied electrode potential and O 2 pressure), so as to be able to make predictions on the optimal experimental conditions for OER [20].After imposing the lattice parameters derived from variable-cell (vc) relaxation of the bulk, slab models were constructed, Cartesian coordinates were relaxed, and surface energies (in J/m 2 ) of each facet were calculated using Equation ( 2) and reported in Figure 5 (top).

Wulff Shapes of Bare Spinel Surfaces as Heterogenous Catalysts
quired for OER is always 1.23 V vs. RHE at all pH values [44].In other words, in calculating the over-potential, the pH term is canceled by taking the difference between predicted onset potentials and the equilibrium OER potential [45].

Results and Discussion
As anticipated in the Introduction, we divide our study into two parts.First, we focus on NiFe2O4, CoFe2O4, NiCo2O4, and ZnCo2O4 systems as bare stoichiometric surfaces: this corresponds to experiments at high temperature and low coverage, i.e., under the conditions of thermal heterogeneous catalysis [9][10][11][12].This will also allow us to compare our predictions with and validate against previous studies [8] (Section 3.1).Then, we switch to the most innovative part of our investigation, i.e., the Wulff construction under electrochemical OER conditions; in Section 3.2, we consider the CoFe2O4 system and study configurations with various coverage of (H, OH, H2O, O, and O2) adsorbate species under electrochemical OER and varying environmental variables (applied electrode potential and O2 pressure), so as to be able to make predictions on the optimal experimental conditions for OER [20].

Wulff Shapes of Bare Spinel Surfaces as Heterogenous Catalysts
For each investigated spinel, the NiFe2O4, CoFe2O4, NiCo2O4, and ZnCo2O4 oxides stoichiometric (AB2O4) but non-symmetric surface slab models were constructed for all considered terminations as shown in Figures 2 and 3.In the case of inverse spinels, NiFe2O4, CoFe2O4, and NiCo2O4, A ++ (M = Ni, Co) cations occupy octahedral sites and B 3+ (M = Fe, Co) cations occupy both octahedral and tetrahedral sites.For the normal spinel ZnCo2O4, A ++ (M = Zn) cations occupy tetrahedral sites and B 3+ cations (M = Co) occupy octahedral sites.Protruding cations on tetrahedral sides (B T and A T ) on the top and bottom surface in the (001) plane were retained to maintain the stoichiometry in both inverse and normal spinels.
After imposing the lattice parameters derived from variable-cell (vc) relaxation of the bulk, slab models were constructed, Cartesian coordinates were relaxed, and surface energies (in J/m 2 ) of each facet were calculated using Equation ( 2) and reported in Figure 5 (top).As is apparent in Figure 5 (top), the calculated surface energies increase in the sequence γ001 < γ111 < γ110 for the three inverse spinels (NiFe 2 O 4 , CoFe 2 O 4 , and NiCo 2 O 4 ).We thus expect the (001) facets to have a larger extension since they exhibit the lowest surface energy in all cases.For the normal spinel ZnCo 2 O 4 , the sequence changes to γ100 < γ110 < γ111.However, the (001) facet is still expected to be dominant, indeed up to the point of being unique, whence a predicted cubic equilibrium shape for ZnCo 2 O 4 .
By using surface energies as calculated above, the predicted equilibrium Wulff configurations of nanoparticle structures were modeled and schematically depicted in Figure 5 (bottom).Fe-based spinels, NiFe 2 O 4 and CoFe 2 O 4 , are shown in Figure 5a,b, respectively, while Co-based spinels, NiCo 2 O 4 and ZnCo 2 O 4 , are shown in Figure 5c,d, respectively.For NiFe 2 O 4 , the predicted Wulff shape consists of (001) facets with a fraction of 81.7% and (111) facets with a fraction of 18.3%, while for CoFe 2 O 4 we find a little fraction of (110) (1.1%) as well, together with 19.1% fraction of (111) and a bigger fraction of (001) (79.8%).For Co-based spinels, the Wulff shape of NiCo 2 O 4 exhibits all three facets with a fraction of 71.9%, 22.7%, and 5.4% for (001), (111), and (110) facets, respectively, whereas in contrast, the Wulff shape of the normal spinel ZnCo 2 O 4 has a pure cube shape without any inclusion of (111) or (110) facets.
As discussed in the introduction, the modeling of spinel surfaces as bare so far refers to low-pressure(vacuum)/low-coverage/high-temperature conditions.In the study of Zasada et al. [8], the authors employed an analogous approach based on periodic DFT+U calculations for a different set of spinel oxide materials and compared their predictions with the results of TEM/STEM experiments on bare nanoparticles.Similar to the present study, the three most stable, (001), (110), and (111), planes exposed by mixed cobalt spinel nanocrystals were investigated [8].In their findings, the abundance of the (110) faces is always low (below 9%) regardless of the nature of the secondary metal in the mixed cobalt spinels.These predictions match our findings very well, predicting a very low abundance of (110) facets in the Wulff nanoparticle shapes as well as predicting the fact that the (001) facet dominates all investigated spinels with a high abundance (>70%).Our predictions should set the ground for a detailed study of the activity of these different materials and their surfaces in heterogeneous catalysis [9][10][11][12]20].

Wulff Shapes of Spinel Surfaces under Electrochemical Conditions
In Section 3.1 we have not included the effects of reaction conditions and of the chemical environment on spinel oxide nanoparticle shapes.Rescaling of surface energies must be considered to obtain Wulff nanoparticle shapes under operative electrochemical conditions, on which we focus in this subsection.Therefore, we investigate OH, H 2 O, and O 2 adsorption on (001), (110), and (111) facets of CoFe 2 O 4 , we predict Wulff reshaping of CoFe 2 O 4 nanoparticles under OER electrochemical conditions, and we compare them with experimental microscopy observations.CoFe 2 O 4 is selected to investigate the change in nanoparticle shape under OER conditions because of its higher OER performance [23].For example, Li et al. prepared various inverse spinel MFe 2 O 4 (M = Co, Ni, Cu, and Mn) nanofibers by the electrospinning technique [46] and found that the OER performance increases in the following order: MnFe 2 O 4 < NiFe 2 O 4 < CuFe 2 O 4 < CoFe 2 O 4 , with CoFe 2 O 4 exhibiting the highest catalytic activity (overpotential, η = 408 mV at 5 mA cm −2 ).Moreover, CoFe 2 O 4 is the system on which there are more experimental characterization results to compare with, as we will discuss in Section 3.2.1.
According to the Wulff-Kaishev theorem [47], for particles deposited on a substrate, the surface energy in the Wulff construction has to be rescaled by the adhesion energy.As adhesion stabilizes the system, the rescaled surface energy will be lower and the corresponding facet will increase in size.Analogously, in the case of particles immersed in a chemical environment, the adsorption energies of ligand species must be included when rescaling surface energies [13,39].Clearly, it is necessary to explore a thorough set of adsorbate configurations to determine the lowest energy state of the system (corresponding to the resting state under reaction conditions).A variety of adsorption coverages and modes of OH, H 2 O, and O 2 on CoFe 2 O 4 slab models was therefore investigated to predict Wulff reshaping.Non-stoichiometric but symmetric slabs were built for the three planes [(001), (110), and (111)] of CoFe 2 O 4 inverse spinel, as shown in Figure 4.
One important initial observation is that for (110) and ( 111), the bare-surface termination ended up with significant distortion of the tetrahedral (under-coordinated) surface Fe atoms upon DFT relaxation.Adsorption of O 2 molecular species was found to solve this issue.An O 2 molecule can indeed be put on (110) as a chelate ligand in bridge adsorption mode in two different positions: (a) on tetrahedral Fe-Fe atoms and (b) on octahedral Co-Fe atoms (see Figure 6a,b).Notably, in Figure 6 for reasons of clarity, the four surface lattice oxygens are labeled as O1, O2, O3, and O4 along with octahedrally coordinated Co (Co O ) and Fe (Fe O ) atoms and tetrahedrally coordinated Fe (Fe T ) atoms.The addition of O 2 in a bridge absorption mode brought about a surface stabilization of −0.28 eV when O 2 was put on (a) Fe-Fe tetrahedral atoms, whereas in the (b) case a destabilization of +0.07 eV was observed.O 2 was thus kept on tetrahedral Fe sites for (110) surface adsorption calculations.On (111), O 2 can be put in bridge between one tetrahedral Fe and one octahedral Co site (see Figure 6c).As a result of the surface cut, one Fe valence coordination and three Co valence coordinations are missing.The addition of O 2 in a bridge absorption mode has the effect of partially completing the coordination of the upper layer; indeed, we found that O 2 adsorption stabilizes the surface by −0.72 eV.Note, however, that while tetrahedral Fe has now completed its four-fold coordination, two coordinated species are still missing for the octahedral Co.
to the resting state under reaction conditions).A variety of adsorption coverages and modes of OH, H2O, and O2 on CoFe2O4 slab models was therefore investigated to predict Wulff reshaping.Non-stoichiometric but symmetric slabs were built for the three planes [(001), (110), and (111)] of CoFe2O4 inverse spinel, as shown in Figure 4.
One important initial observation is that for (110) and ( 111), the bare-surface termination ended up with significant distortion of the tetrahedral (under-coordinated) surface Fe atoms upon DFT relaxation.Adsorption of O2 molecular species was found to solve this issue.An O2 molecule can indeed be put on (110) as a chelate ligand in bridge adsorption mode in two different positions: (a) on tetrahedral Fe-Fe atoms and (b) on octahedral Co-Fe atoms (see Figure 6a,b).Notably, in Figure 6 for reasons of clarity, the four surface lattice oxygens are labeled as O1, O2, O3, and O4 along with octahedrally coordinated Co (Co O ) and Fe (Fe O ) atoms and tetrahedrally coordinated Fe (Fe T ) atoms.The addition of O2 in a bridge absorption mode brought about a surface stabilization of −0.28 eV when O2 was put on (a) Fe-Fe tetrahedral atoms, whereas in the (b) case a destabilization of +0.07 eV was observed.O2 was thus kept on tetrahedral Fe sites for (110) surface adsorption calculations.On (111), O2 can be put in bridge between one tetrahedral Fe and one octahedral Co site (see Figure 6c).As a result of the surface cut, one Fe valence coordination and three Co valence coordinations are missing.The addition of O2 in a bridge absorption mode has the effect of partially completing the coordination of the upper layer; indeed, we found that O2 adsorption stabilizes the surface by −0.72 eV.Note, however, that while tetrahedral Fe has now completed its four-fold coordination, two coordinated species are still missing for the octahedral Co.Then, as in our previous investigation of the OER catalysis on the (100) facet of CoFe2O4 [23], undissociated water molecules were used to fill in the missing coordination of metal sites on each surface [(001), (110), and (111)].In a further step, the transformation of waters into hydroxyls and hydrogen adatoms via dissociation (*H2O → *OH + *H) was considered and the affinity of *H2O, *OH, and *H on the metal sites were predicted by calculating the relative energy of the configurations associated with each coverage pattern.The sampling of several coverage patterns on the catalyst surface allowed us to determine the lowest energy state (i.e., the resting state under reaction conditions) on each of the three facets of CoFe2O4.Notably, along with stoichiometric coverage (i.e., dissociated and/or un-dissociated water coverage), off stoichiometric patterns with excess hydrogens on surface oxygens as well as deprotonated surfaces were also considered.The complete set of studied structures is shown and discussed in Figures S2-S6 of the SI.Here, we only show the configurations which we predict as resting states under selected applied potentials, i.e., U = 1.23 V, 1.48 V, and 1.63 V vs. SHE.After adsorption of reaction intermediates, the surface energies were recalculated by adding adsorption energies (E ads ) as formalized in Equation ( 5).Notably, free-energy contributions of reaction species must be added to the DFT energetics and included into adsorption energy as shown in Equation ( 4).Finally, E ads terms including free energy contributions must be divided by the area of adsorbed surfaces for normalization.
The lowest energetic structures (i.e., the resting states) resulting from these calculations under different applied potential (U = 1.23 V, 1.48 V, 1.63 V vs. SHE) spanning a realistic set of reaction conditions for OER are shown in Figure 7A-C.
tributions of reaction species must be added to the DFT energetics and included into adsorption energy as shown in Equation ( 4).Finally, E terms including free energy contributions must be divided by the area of adsorbed surfaces for normalization.
The lowest energetic structures (i.e., the resting states) resulting from these calculations under different applied potential (U = 1.23 V, 1.48 V, 1.63 V vs. SHE) spanning a realistic set of reaction conditions for OER are shown in Figure 7A-C.In Figure 7B, four selected resting states for (001), (110), and (111) under bias U = 1.48 V (η = 0.25) are shown.For the (001) surface, there are two resting states corresponding to minimum energy (i.e., dissociated waters on Fe and Co sites, respectively, as shown in Figure S2b,c of the SI) that are shown in Figure 7B (configurations d1, d2).As for (110), the deprotonated pattern with respect to two waters gave the lowest energy at U > 1.24 V among adsorbed patterns.For (111), the same resting state as at the U = 1.23 V case (two water coordinated to Co cations) was predicted.As a next step, the energies of resting state configurations under U = 1.23 V, 1.48 V, and 1.63 V vs. SHE (shown in Figure 7) were used to recalculate the surface energies based on Equation ( 5) and thus to estimate the Wulff nanoparticle shapes of CoFe 2 O 4 under the given bias.
In the table in Figure 8A, a comparison of surface energies calculated from nonadsorbed (i.e., bare) surfaces with the surface energies recalculated for adsorbed surfaces under applied potentials (U = 1.23 V, 1.48 V, and 1.63 V vs. SHE) for CoFe 2 O 4 inverse spinel is reported.As is apparent in this table, the adsorption process reduces the value of the surface energy for each surface, as expected.The recalculated surface energies for adsorbed facets now increase in the sequence γ111 < γ001 < γ110, which is different from the sequence for bare surfaces: γ001 < γ111 < γ110.We thus predict an inversion from low coverage to OER conditions, i.e., we expect to have a larger extension of (111) in the Wulff nanoparticle shape under electrochemical conditions.Notably, we did not include the reaction energy to transform the stoichiometric bare facets to non-stoichiometric covered facets in our estimates; this reaction involves desorption of Fe cations from the stoichiometric facets, whose energetics is difficult to estimate as it will depend upon environmental conditions (the chemical potential of Fe cations).However, in both free and OER conditions, the (110) facets have the highest surface energy (therefore the least stable plane) for spinel surfaces.This is consistent with experimental results, see the discussion below.

Comparison with the Experiment
In order to validate our theoretical predictions, we now proceed to a comparison available information from previously published experimental work on spinel oxid noparticles in which microscopic characterization of the systems under OER cond was provided and which therefore lend themselves best to be compared with our W construction predictions.To this purpose, it would be ideal to have microscopy mea ments on CoFe2O4 nanoparticles detected under OER operando conditions.However, best of our knowledge, there are no examples in the literature of such operando ex ments.Therefore, we will compare four experimental studies from the literature in w TEM and HRTEM images of as-prepared CoFe2O4 nanoparticles or of CoFe2O4 nano Using Equation ( 5) and the WulffPack-1.2python package [29], the surface energy values in Figure 8A were used to obtain Wulff nanoparticle shapes of CoFe 2 O 4 under applied potentials and the results are shown in Figure 8B.We find that under OER applied potentials, (111) dominates in all cases with a fraction higher than 90%, while (001) facets occupy less than 10% and (110) facets do not contribute to the Wulff shape.Note that, as the U value increases, (111) occupies a larger and larger area: at U = 1.23 V, (111) has a fraction of 90.2%, then passing to a 92.8% fraction at U = 1.48 V, and ending up at 97.1% at U = 1.63 V.However, these changes, as a function of the bias U, are not drastic and may be difficult to observe experimentally.

Comparison with the Experiment
In order to validate our theoretical predictions, we now proceed to a comparison with available information from previously published experimental work on spinel oxide nanoparticles in which microscopic characterization of the systems under OER conditions was provided and which therefore lend themselves best to be compared with our Wulff construction predictions.To this purpose, it would be ideal to have microscopy measurements on CoFe 2 O 4 nanoparticles detected under OER operando conditions.However, to the best of our knowledge, there are no examples in the literature of such operando experiments.Therefore, we will compare four experimental studies from the literature in which TEM and HRTEM images of as-prepared CoFe 2 O 4 nanoparticles or of CoFe 2 O 4 nanoparticles after some cycles of electrochemical treatment have been reported.
Gebrelase et al. observed a transformation of the CoFe 2 O 4 structure by varying and optimizing the ratio of CoFe 2 O 4 and dopamine contents in a peculiar synthesis recipe [48], reporting a final over-potential of 440 mV at 10 mA/cm 2 for pristine CoFe 2 O 4 .These CoFe 2 O 4 samples were prepared via the synthesis of spinel oxide nanoparticles by the hydrothermal method, addition and polymerization of dopamine in a nanoparticle suspension, and finally followed by carbonization via calcination, as described in the original paper [48].Regarding the CoFe 2 O 4 nanoparticle shape (their size was 3 to 45 nm in diameter), using HRTEM imaging, Gebrelase et al. reported an octahedron-like structure for pristine CoFe 2 O 4 before the OER.The experimental TEM and HRTEM images from their work are reported in Figure 9a,b and clearly show dominance of the (111) surface, matching the predicted Wulff shapes in Figure 8b-d.
Physchem 2024, 4, FOR PEER REVIEW 13 of 432 mV at 10 mA/cm 2 (1.66 V vs. RHE) and reported TEM images of CoFe2O4 nanoparticles that are reported in Figure 9c.Comparing microscopy images of pristine CoFe2O4 nanoparticles in Ref. [49] with our calculated nanoparticles at U = 1.63 V in Figure 8c, one notice that the TEM images often show a more spherical shape than that predicted by the present study.However, in the experiments one can find more faceted nanoparticles as well; for example, the nanoparticle enclosed in a yellow circle in Figure 9c exhibits (111) planes due to its prismatic/sharp edges and compares very favorably with our predictions.Arrassi et al. analyzed 4 nm-sized CoFe2O4 spinel nanoparticles and showed that these particles retain their size and crystal structure after OER as observed via area electron diffraction measurements (SAED) [51].CoFe2O4 was tested in its intrinsic catalytic response without binders or additives and OER was recorded at a potential of 1.86 V vs. RHE.They also detected (111) faceting in nanoparticles with SAED analysis before and after electrochemistry experiments and observed no structural changes after OER, in agreement with theoretical calculations under a wide range of bias in the present study.
Finally, in the study of Kargar et al. [52], the sample morphology and material composition of CoFe2O4 nanoparticles on carbon fiber papers (NPs-on-CFPs) were examined after long-term stability testing under OER.After more than 15 h of long cycle testing (>1000 cycles), the sample was thoroughly examined using SEM imaging, XRD, and elemental mapping analyses-we refer to the original work [52] for the details for the experimental analysis.Similarly to previous studies, they concluded that the morphology of CoFe2O4 NPs-on-CFP did not change significantly; the samples showed long term stability without any morphological or compositional modifications.

Effect of Oxygen Pressure on Nanoparticle Shape
We noted above that O2 adsorption has a role in stabilizing the (110) and (111) surfaces of CoFe2O4.The physical reason of this unusually high adsorption energy lies in the charge state of the O2 molecule.We have performed a Bader analysis of the configuration depicted in Figure 6c, which exhibits the highest O2 desorption energy, and found that   [50].Also, they reported an OER activity of pristine CoFe 2 O 4 with an overpotential of 432 mV at 10 mA/cm 2 (1.66 V vs. RHE) and reported TEM images of CoFe 2 O 4 nanoparticles that are reported in Figure 9c.Comparing microscopy images of pristine CoFe 2 O 4 nanoparticles in Ref. [49] with our calculated nanoparticles at U = 1.63 V in Figure 8c, one notice that the TEM images often show a more spherical shape than that predicted by the present study.However, in the experiments one can find more faceted nanoparticles as well; for example, the nanoparticle enclosed in a yellow circle in Figure 9c exhibits (111) planes due to its prismatic/sharp edges and compares very favorably with our predictions.
Arrassi et al. analyzed 4 nm-sized CoFe 2 O 4 spinel nanoparticles and showed that these particles retain their size and crystal structure after OER as observed via area electron diffraction measurements (SAED) [51].CoFe 2 O 4 was tested in its intrinsic catalytic response without binders or additives and OER was recorded at a potential of 1.86 V vs. RHE.They also detected (111) faceting in nanoparticles with SAED analysis before and after electrochemistry experiments and observed no structural changes after OER, in agreement with theoretical calculations under a wide range of bias in the present study.
Finally, in the study of Kargar et al. [52], the sample morphology and material composition of CoFe 2 O 4 nanoparticles on carbon fiber papers (NPs-on-CFPs) were examined after long-term stability testing under OER.After more than 15 h of long cycle testing (>1000 cycles), the sample was thoroughly examined using SEM imaging, XRD, and elemental mapping analyses-we refer to the original work [52] for the details for the experimental analysis.Similarly to previous studies, they concluded that the morphology of CoFe 2 O 4 NPs-on-CFP did not change significantly; the samples showed long term stability without any morphological or compositional modifications.

Effect of Oxygen Pressure on Nanoparticle Shape
We noted above that O 2 adsorption has a role in stabilizing the (110) and (111) surfaces of CoFe 2 O 4 .The physical reason of this unusually high adsorption energy lies in the charge state of the O 2 molecule.We have performed a Bader analysis of the configuration depicted in Figure 6c, which exhibits the highest O 2 desorption energy, and found that each oxygen atom bears a charge of −0.4 e, for a total of a charge of −0.8 e for the O 2 molecule.For reference, the oxygen anions in the bulk, which formally correspond to doubly negatively charged species, bear a Bader charge of −1.Given this sizeable adsorption energy, in this section, we will explore whether it is possible to change the nanoparticle shape by controlling the O 2 pressure.We use standard thermodynamics to calculate the Gibbs free energy of O 2 in actual conditions (∆G), which is related to the Gibbs free energy in the standard state (∆G • ) by the following relationship: Here, Q r is the reaction quotient which allows one to estimate the changes (i.e., temperature, concentration, pressure, etc.) under non-standard conditions.The effect of pressure was therefore included by adding the last term in Equation ( 6) to the adsorption energy (E ads ) calculations, see Equation ( 5) for ( 110) and (111) facets.We focused on a bias of U = 1.48 V (also considering 1.63 V in the SI) and investigated both increasing oxygen pressure up to 30 atm (1 atm, 5 atm, 15 atm, and 30 atm) and also decreasing oxygen pressure down to 10 −4 atm.Notably, the (001) surface energy was kept at 0.75 J/m 2 and 0.71 J/m 2 at U = 1.48 V and 1.63 V, respectively, as this facet does not carry O 2 in its resting state.Rescaled surface energies and Wulff constructions are reported in Figure 10 at U = 1.48 V, as this bias is presently the optimal target of only 0.25 V over-potential and is considered as the one realistically closest to the minimum voltage necessary for water electrolysis (1.23 V).Results at U =1.63 V are additionally shown in Figure S8 of the Supplementary Materials.Figure 10 clearly demonstrates that changing the oxygen pressure has a dramatic effect on the surface energies and the Wulff shape, with the (111) facet favored at high O 2 pressure and the (100) facet favored at low O 2 pressure.In contrast, it can be noted that, despite these dramatic changes, the (110) surface energies are still too high to make the contribution of this facet significant compared to other surfaces of CoFe 2 O 4 .Under 15 to 30 atm of O 2 pressure at U = 1.43-1.63V, the surface energies of (110) were computed lower than (001); nonetheless, the (111) surface energy in parallel also decreases significantly, so that, eventually, we predict a full fraction 100% of (111) in a purely octahedron shape (see Figure 10e).At the opposite, when the O 2 pressure is decreased down to 0.5 atm, the surface energy of (111) approaches the value of (001) and the (001) area on the nanoparticles gradually increases.Dramatically, when the O 2 pressure is reduced to very low values of 10 −4 atm, the nanoparticle surfaces become a full fraction of 100% of (001) with a purely cube shape, as shown in Figure 10a.In light of the above findings, we conclude that increasing O2 pressure might poison the catalyst in the OER process.Indeed, increasing O2 pressure increases the size of the higher-index facets, i.e., (111) and (110), at the expense of the low-index (001) facet, whereas decreasing O2 pressure has the opposite effect.Now, although the (111) facets is catalytically active [24][25][26], its activity can be poisoned by the adsorption of the O2 reaction product, whereas the catalytic activity of the (001) facet may be higher and should not be affected by O2 pressure [23].These findings thus suggest working experimentally at the lowest O2 pressure as possible in order to maximize the OER catalytic activity of CoFe2O4, generally spinel oxide, materials.According to our DFT modeling, in fact, the (001) facet should completely dominate the nanoparticle shape at an O2 pressure of 10 −4 atm.

Conclusions
Our aim in this work is to predict and elucidate nanoparticle shape under diverse environments, for both validation and production purposes, to provide ground for operando monitoring of these promising catalysts and compare them with available experimental data.We are especially interested in predicting from first-principle simulation na- In light of the above findings, we conclude that increasing O 2 pressure might poison the catalyst in the OER process.Indeed, increasing O 2 pressure increases the size of the higher-index facets, i.e., (111) and (110), at the expense of the low-index (001) facet, whereas decreasing O 2 pressure has the opposite effect.Now, although the (111) facets is catalytically active [24][25][26], its activity can be poisoned by the adsorption of the O 2 reaction product, whereas the catalytic activity of the (001) facet may be higher and should not be affected by O 2 pressure [23].These findings thus suggest working experimentally at the lowest O 2 pressure as possible in order to maximize the OER catalytic activity of CoFe 2 O 4 , generally spinel oxide, materials.According to our DFT modeling, in fact, the (001) facet should completely dominate the nanoparticle shape at an O 2 pressure of 10 −4 atm.

Conclusions
Our aim in this work is to predict and elucidate nanoparticle shape under diverse environments, for both validation and production purposes, to provide ground for operando monitoring of these promising catalysts and compare them with available experimental data.We are especially interested in predicting from first-principle simulation nanoparticle shapes under OER electrochemical conditions, which, to the best of our knowledge, has not been attempted thus far.To achieve this goal, we systematically investigated the surface structures of the low-index facets [(001), (110), and (111)] of a set of spinel oxides (NiFe 2 O 4 , CoFe 2 O 4 , NiCo 2 O 4 , and ZnCo 2 O 4 ) using periodic DFT+U calculations under two different sets of conditions, corresponding to the two possible catalytic applications of these systems: (i) as bare surfaces under vacuum or low-coverage high-temperature conditions typical of chemical looping or oxidation reactions [9][10][11][12] and (ii) as adsorbate-covered facets under OER (limiting the CoFe 2 O 4 system as a paradigmatic example) [20][21][22][23].The calculated surface energies were then used in the Wulff construction to predict the equilibrium nanoparticle shapes of spinel oxides under different environmental conditions.
For bare surfaces under vacuum, regardless of the nature of the alien metal in the mixed iron and cobalt spinels, we found that the abundance of the (110) faces is always low and that (001) facets dominate the nanoparticle shapes, in agreement with the experiment and a previous systematic study [8].
CoFe 2 O 4 , a promising non-critical-material candidate for OER/AEMWE catalysis, was then selected to comprehensively study the adsorbate-covered resting state surface structures of low index surfaces [(001), (110), and (111)] under OER conditions, exploiting results of a previous study on OER catalysis [23].
First, we found that under OER, the ordering of surface energies changes with respect to bare systems from γ001 < γ111 < γ110 to γ111 < γ001 < γ110 under a wide range of realistic applied potentials.Then, building the optimal nanoparticle shape via the Wulff construction, we predict that CoFe 2 O 4 exhibit dominant (111) facets up to 1.63 V bias vs. SHE, in fair agreement with experimental TEM and HRTEM observations.Second, we investigated the effect of oxygen pressure to the best of our knowledge for the first time herein.We predict that, upon O 2 adsorption, a stabilization takes place for the (110) and especially the (111) surfaces of CoFe 2 O 4 .We rationalize this prediction in terms of the formation of a strongly adsorbed superoxide-like species on (110) and (111).Dramatic changes in the Wulff shapes are predicted: at high O 2 pressure, the (111) facet will dominate in a purely octahedron nanoparticle shape, whereas, at low O 2 pressure, the (001) will eventually become dominant, to achieve a purely cubic nanoparticle shape.Whereas (001) is a catalytically active facet whose activity does not depend on O 2 pressure [23], (111) is catalytically active [24][25][26] but may be poisoned by the adsorption of the reaction product; this implies that O 2 can poison the catalysts and ideal OER conditions may correspond to the lowest O 2 pressure possible.
In perspective, we underline the importance of investigating computationally nanoparticle reshaping of spinel oxides under operating conditions together with experimental operando characterization of these systems.This should open the way to cross-validation between theory and experiment, which should lead to a deeper understanding and then, eventually, to the rational design of much-needed more efficient and sustainable optimized catalysts in both oxidation and OER catalysis.In perspective, it will also be very interesting to further extend our approach by including defects such as oxygen vacancies in the surface energetics [53] or more complex composite architectures [54,55].

Figure 1 .
Figure 1.Examples of Wulff constructions for shapes of high cubic symmetry.

Figure 1 .
Figure 1.Examples of Wulff constructions for shapes of high cubic symmetry.

Figure 2 .
Figure 2. Stoichiometric but non-symmetric surface slab models for the considered surface terminations of inverse spinel (AB2O4) systems in the ball and stick representation: (100) (a), (110) (b), and (111) (c).O and T denote octahedral and tetrahedral coordination of A and B metals, respectively.Conditions of thermal catalysis.

Figure 3 .
Figure 3. Stoichiometric but non-symmetric surface slab models for the considered surface terminations of normal spinel (AB2O4) systems in the ball and stick representation: (100) (a), (110) (b), and (111) (c).T and O denote tetrahedral and octahedral coordination of A and B metals, respectively.Conditions of thermal catalysis.

Figure 2 .
Figure 2. Stoichiometric but non-symmetric surface slab models for the considered surface terminations of inverse spinel (AB 2 O 4 ) systems in the ball and stick representation: (100) (a), (110) (b), and (111) (c).O and T denote octahedral and tetrahedral coordination of A and B metals, respectively.Conditions of thermal catalysis.

Figure 2 .
Figure 2. Stoichiometric but non-symmetric surface slab models for the considered surface terminations of inverse spinel (AB2O4) systems in the ball and stick representation: (100) (a), (110) (b), and (111) (c).O and T denote octahedral and tetrahedral coordination of A and B metals, respectively.Conditions of thermal catalysis.

Figure 3 .
Figure 3. Stoichiometric but non-symmetric surface slab models for the considered surface terminations of normal spinel (AB2O4) systems in the ball and stick representation: (100) (a), (110) (b), and (111) (c).T and O denote tetrahedral and octahedral coordination of A and B metals, respectively.Conditions of thermal catalysis.

Figure 3 .
Figure 3. Stoichiometric but non-symmetric surface slab models for the considered surface terminations of normal spinel (AB 2 O 4 ) systems in the ball and stick representation: (100) (a), (110) (b), and (111) (c).T and O denote tetrahedral and octahedral coordination of A and B metals, respectively.Conditions of thermal catalysis.

Figure 4 .
Figure 4. Non-stoichiometric/symmetric slab models for the considered surface terminations of inverse spinel systems in the following ball and stick representations: (a) (001), (b) (110), and (c) (111).O and T denote octahedral and tetrahedral coordination of A(Co) and B(Fe) metals, respectively.Conditions of electro-catalysis.

Figure 4 .
Figure 4. Non-stoichiometric/symmetric slab models for the considered surface terminations of inverse spinel systems in the following ball and stick representations: (a) (001), (b) (110), and (c) (111).O and T denote octahedral and tetrahedral coordination of A(Co) and B(Fe) metals, respectively.Conditions of electro-catalysis.The total number of atoms per unit cell for (001), (110), and (111) slabs used in bare and adsorbed surface calculations is summarized in Table S2 of the Supplementary Materials; brief details about the comparison of those slabs are shown in Figure S1.To calculate surface energies for surfaces with adsorbates, we re-scale them as shown in Equation(3), where adsorption energies are taken from optimization calculations of non-stoichiometric slab models[13,39]: For each investigated spinel, the NiFe 2 O 4 , CoFe 2 O 4 , NiCo 2 O 4 , and ZnCo 2 O 4 oxides stoichiometric (AB 2 O 4 ) but non-symmetric surface slab models were constructed for all considered terminations as shown in Figures 2 and 3.In the case of inverse spinels, NiFe 2 O 4 , CoFe 2 O 4 , and NiCo 2 O 4 , A ++ (M = Ni, Co) cations occupy octahedral sites and B 3+ (M = Fe, Co) cations occupy both octahedral and tetrahedral sites.For the normal spinel ZnCo 2 O 4 , A ++ (M = Zn) cations occupy tetrahedral sites and B 3+ cations (M = Co) occupy octahedral sites.Protruding cations on tetrahedral sides (B T and A T ) on the top and bottom surface in the (001) plane were retained to maintain the stoichiometry in both inverse and normal spinels.

Figure 6 .
Figure 6.O 2 -bridged (110) CoFe 2 O 4 surfaces sited on (a) tetrahedral Fe-Fe atoms and (b) octahedral Co-Fe atoms.(c) O 2 -bridged (111) CoFe 2 O 4 surface sited on octahedral Co-tetrahedral Fe atoms.Then, as in our previous investigation of the OER catalysis on the (100) facet of CoFe 2 O 4[23], undissociated water molecules were used to fill in the missing coordination of metal sites on each surface [(001), (110), and (111)].In a further step, the transformation of waters into hydroxyls and hydrogen adatoms via dissociation (*H 2 O → *OH + *H) was considered and the affinity of *H 2 O, *OH, and *H on the metal sites were predicted by calculating the relative energy of the configurations associated with each coverage pattern.The sampling of several coverage patterns on the catalyst surface allowed us to determine the lowest energy state (i.e., the resting state under reaction conditions) on each of the three facets of CoFe 2 O 4. Notably, along with stoichiometric coverage (i.e., dissociated and/or un-dissociated water coverage), off stoichiometric patterns with excess hydrogens on surface oxygens as well as deprotonated surfaces were also considered.The complete set of studied structures is shown and discussed in FiguresS2-S6of the SI.Here, we only show the configurations which we predict as resting states under selected applied potentials, i.e., U = 1.23 V, 1.48 V, and 1.63 V vs. SHE.After adsorption of reaction intermediates, the surface energies were recalculated by adding adsorption energies (E ads ) as formalized in Equation (5).Notably, free-energy contributions of reaction species must be added to the DFT energetics and included into adsorption energy as shown in Equation (4).Finally, E ads terms including free energy contributions must be divided by the area of adsorbed surfaces for normalization.

Figure 7 .
Figure 7. Resting-state configurations of (001), (110), and (111) surfaces of CoFe2O4 at (A) U = 1.23 V, at (B) U = 1.48 V, and (C) U = 1.63 V vs. SHE for the (001) = (a,d,g); (110) = (b,e,h), and (111) = (c,f,l) facets, respectively.Oxygen, hydrogen, iron, and cobalt atoms are colored red, white, violet, and indigo-blue, respectively.The three selected resting states for (001), (110), and (111) under bias U = 1.23 V (η = 0.0) are shown in panel Figure 7A.As for the (001) surface, water adsorption on metal sites along with excess hydrogens on lattice oxygens gave the lowest energy at U = 1.23 V.For (110), the pattern (b) which has an *OH adsorbed on a bridge region between octahedral Co and Fe atoms along with two hydrogens on lattice oxygens (O1 and O4) gave the lowest energy.For the (111) facet, the state in which two water adsorbates complete the surface Co s missing two octahedral coordination was predicted as the resting state under U= 1.23 V, shown as (c) in Figure 7A (see also panel (a) in Figure S7 of the Supplementary Materials).

Figure 7 .
Figure 7. Resting-state configurations of (001), (110), and (111) surfaces of CoFe 2 O 4 at (A) U = 1.23 V, at (B) U = 1.48 V, and (C) U = 1.63 V vs. SHE for the (001) = (a,d,g); (110) = (b,e,h), and (111) = (c,f,i) facets, respectively.Oxygen, hydrogen, iron, and cobalt atoms are colored red, white, violet, and indigo-blue, respectively.The three selected resting states for (001), (110), and (111) under bias U = 1.23 V (η = 0.0) are shown in panel Figure 7A.As for the (001) surface, water adsorption on metal sites along with excess hydrogens on lattice oxygens gave the lowest energy at U = 1.23 V.For (110), the pattern (b) which has an *OH adsorbed on a bridge region between octahedral Co and Fe atoms along with two hydrogens on lattice oxygens (O1 and O4) gave the lowest energy.For the (111) facet, the state in which two water adsorbates complete the surface Co's missing two octahedral coordination was predicted as the resting state under U= 1.23 V, shown as (c) in Figure 7A (see also panel (a) in Figure S7 of the Supplementary Materials).In Figure7B, four selected resting states for (001), (110), and (111) under bias U = 1.48 V (η = 0.25) are shown.For the (001) surface, there are two resting states corresponding to minimum energy (i.e., dissociated waters on Fe and Co sites, respectively, as shown in FigureS2b,c of the SI) that are shown in Figure7B(configurations d1, d2).As for (110), the deprotonated pattern with respect to two waters gave the lowest energy at U > 1.24 V among adsorbed patterns.For (111), the same resting state as at the U = 1.23 V case (two water coordinated to Co cations) was predicted.

Figure 9 .
Figure 9. (a) TEM image of pristine CoFe2O4 nanoparticles.(b) HRTEM images with and lattice fringe analyses of pristine CoFe2O4.(a,b) Figures adapted from Ref[48].Copyright© Elsevier 2022 (c) TEM image of the pristine CoFe2O4 nanoparticles in the work of Xiang et al. Figure adapted from Ref [49].Copyright© Springer Nature 2022.
Figure 9. (a) TEM image of pristine CoFe2O4 nanoparticles.(b) HRTEM images with and lattice fringe analyses of pristine CoFe2O4.(a,b) Figures adapted from Ref[48].Copyright© Elsevier 2022 (c) TEM image of the pristine CoFe2O4 nanoparticles in the work of Xiang et al. Figure adapted from Ref [49].Copyright© Springer Nature 2022.

Figure 9 .
Figure 9. (a) TEM image of pristine CoFe 2 O 4 nanoparticles.(b) HRTEM images with and lattice fringe analyses of pristine CoFe 2 O 4 .(a,b) Figures adapted from Ref[48].Copyright© Elsevier 2022 (c) TEM image of the pristine CoFe 2 O 4 nanoparticles in the work of Xiang et al. Figure adapted from Ref [49].Copyright© Springer Nature 2022.
Figure 9. (a) TEM image of pristine CoFe 2 O 4 nanoparticles.(b) HRTEM images with and lattice fringe analyses of pristine CoFe 2 O 4 .(a,b) Figures adapted from Ref[48].Copyright© Elsevier 2022 (c) TEM image of the pristine CoFe 2 O 4 nanoparticles in the work of Xiang et al. Figure adapted from Ref [49].Copyright© Springer Nature 2022.Xiang et al. used atom probe tomography (APT) to elucidate the evolution of the 3D structure of 10 nm-sized Co 2 FeO 4 and CoFe 2 O 4 nanoparticles during OER [49].They 2 e.The adsorbed O 2 thus corresponds essentially to a superoxide anion, O 2 -.A perusal of bond distances confirms this analysis: the O-O distance in the configuration depicted in Figure 6c is 1.36 Å, to be compared with a distance of 1.22 Å for the (neutral) O 2 molecule, 1.35 Å in the (singly negative) superoxide LiO 2 molecule, and 1.58 Å in the (doubly negative) peroxide Li 2 O 2 , respectively.
), (111)] planes of CoFe 2 O 4 , and Wulff shapes as a function of oxygen pressure.Table S1.Calculated lattice parameters (a,b) of investigated spinels in Angstrom (Å).c = 20 Å for all facets.Table S2.The number of atoms in the unit cell used in bare and adsorbed surface calculations of CoFe 2 O 4 .Figure S1.Stoichiometric (a,c,e) and Nonstoichiometric slabs (b,d,f) of inverse spinel (AB 2 O 4 ) used for bare and covered surfaces of CoFe 2 O 4 , respectively.Symmetry was protected for all slab models.Figure S2.Water coverage patterns on CoFe 2 O 4 (001).(a) Surface covered with adsorbed undissociated water molecules on each metal sites and nonprotonated lattice oxygens.(b) Adsorbed OH (*OH) on the Co site, while water is still adsorbed (*H2O) on the Fe site, and one lattice surface oxygen is protonated.(c) Adsorbed OH (*OH) on the Fe site, while water is still adsorbed (*H2O) on the Co site, and one lattice surface oxygen is protonated.(d) Adsorbed *OH on both Fe and Co metal sites, while two lattice surface oxygens are protonated.Oxygen, hydrogen, iron, and cobalt atoms are colored red, white, violet, and indigo blue, respectively.Figure S3.Offstoichiometric coverage patterns on CoFe 2 O 4 (001).(a-c) Surface covered with excess hydrogens on lattice oxygens.(d,e) Missing hydrogens with respect to two water coverage.Free energies referenced to the resting states in Figure S2.(b,c).Energetics with respect to applied potentials (U = 1.23 V, 1.48 V and 1.63 V) has been reported each converged structure.The energetics are with respect to Figure S2 (b,c) configurations.Oxygen, hydrogen, iron, and cobalt atoms are colored red, white, violet, and indigo blue, respectively.Figure S4.Dissociated water Coverage patterns on CoFe 2 O 4 (110).Oxygen, hydrogen, iron, and cobalt atoms are colored red, white, violet, and indigo blue, respectively.Figure S5.Dissociated water Coverage patterns on CoFe 2 O 4 (110).Oxygen, hydrogen, iron, and cobalt atoms are colored red, white, violet, and indigo blue, respectively Figure S6.Off-stoichiometric coverage patterns on CoFe 2 O 4 (110) with (a-e) excess hydrogens on lattice oxygens, and (f) with missing hydrogen.Figure S7.Water coverage (a-c), and off-stoichiometric coverage patterns on CoFe 2 O 4 (111) with (d) excess hydrogens on lattice oxygens and (e) missing hydrogen.Figure S8.Calculated DFT Surface Energies with oxygen pressure (10−4 to 30 atm) on (001), (110) and (111) surfaces of CoFe 2 O 4 under applied potential U = 1.63 V vs SHE (above).Wullf NP shapes of CoFe 2 O 4 under U= 1.63 V, ranged from 10−4 atm (a) to 30 atm (e).Color codes: blue for (111), yellow for (001) facets.