Computationally Accelerated Discovery and Experimental Demonstration of Gd0.5La0.5Co0.5Fe0.5O3 for Solar Thermochemical Hydrogen Production

Solar thermochemical hydrogen (STCH) production is a promising method to generate carbon neutral fuels by splitting water utilizing metal oxide materials and concentrated solar energy. The discovery of materials with enhanced water-splitting performance is critical for STCH to play a major role in the emerging renewable energy portfolio. While perovskite materials have been the focus of many recent efforts, materials screening can be time consuming due to the myriad chemical compositions possible. This can be greatly accelerated through computationally screening materials parameters including oxygen vacancy formation energy, phase stability, and electron effective mass. In this work, the perovskite Gd0.5La0.5Co0.5Fe0.5O3 (GLCF), was computationally determined to be a potential water splitter, and its activity was experimentally demonstrated. During water splitting tests with a thermal reduction temperature of 1,350°C, hydrogen yields of 101 μmol/g and 141 μmol/g were obtained at re-oxidation temperatures of 850 and 1,000°C, respectively, with increasing production observed during subsequent cycles. This is a significant improvement from similar compounds studied before (La0.6Sr0.4Co0.2Fe0.8O3 and LaFe0.75Co0.25O3) that suffer from performance degradation with subsequent cycles. Confirmed with high temperature x-ray diffraction (HT-XRD) patterns under inert and oxidizing atmosphere, the GLCF mainly maintained its phase while some decomposition to Gd2-xLaxO3 was observed.


INTRODUCTION
Solar thermochemical hydrogen (STCH) production has been studied as a potential path to produce alternative fuels (Miller et al., 2014). This reaction is generally a two-step process that generates hydrogen (H 2 ) gas by splitting water using metal oxide materials and concentrated solar energy (Steinfeld, 2005;Smestad and Steinfeld, 2012). In the first step of a typical STCH cycle, the metal oxide is thermally reduced under inert atmosphere at high temperature (>1,200°C, achieved using concentrated solarthermal flux), creating oxygen vacancies in the metal oxide and releasing oxygen gas (Scheffe and Steinfeld, 2014). In the second re-oxidation step at a lower temperature under steam, the oxygen-deficient metal oxide splits water to produce H 2 gas while regenerating the metal oxide for consecutive water splitting cycles. While a combined photovoltaic/electrolytic system can be an alternative method for generating H 2 from water (Ivy, 2004), the efficiency of STCH systems were predicted to potentially exceed the efficiency of the combined system (Siegel et al., 2013).
With the vast composition space of inorganic materials, preliminary computational materials screening has become an important tool for accelerating materials discovery. Computational methods were previously applied to materials screening for the solid-state hydrogen storage reaction and chemical looping process (Clary et al., 2020;Singstock et al., 2020). For STCH, ternary (ABO 3 ) and quaternary (AA'BO 3 ) perovskites were explored to computationally evaluate their oxygen vacancy formation energies, electronic properties, and thermodynamic stabilities, with several promising candidate materials predicted based on these results (Emery et al., 2016;Sai Gautam et al., 2020). In this work, we report Gd 0.5 La 0.5 Co 0.5 Fe 0.5 O 3 (GLCF) as a new STCH material that was discovered using a computational screening approach and experimentally demonstrated as a STCH producing material.

Computational Screening Framework
Potential STCH compounds from the A 2 BB'O 6 , AA'B 2 O 6 , and AA'BB'O 6 compositional spaces were first screened for stability as perovskites using the machine learned descriptor τ (Bartel et al., 2019). τ classifies potential perovskite compositions as perovskite or non-perovskite using the Shannon radii (rA, rB, rX) of the A, B, and X site ions and the formal oxidation state of the A site cation as inputs. This descriptor exhibits 92% accuracy for predicting theoretical perovskite synthesizability for ABX 3 compositions and 91% accuracy for A 2 BB'X 6 compositions, where fractional weighting of the B site cation radii is used (Bartel et al., 2019). Of the compounds identified by τ as synthesizable as perovskites, the STCH relevant properties of >1,000 Gd-containing compositions were evaluated in a highthroughput optimization scheme (Bare et al., 2021a). The compositions La 2 CoFeO 6 , GdLaCoFeO 6 , and Gd 2 CoFeO 6 are predicted by τ to be stable as perovskites (with formal oxidation states of +3 used for Gd, La, Co, and Fe and -2 used for O) and are predicted by DFT to have favorable STCH properties (see Results and Discussion).
PySPuDS (Bare et al., 2021b), a custom high-throughput python wrapper for the bond valence method (BVM) based Structure Prediction and Diagnostic Software (SPuDS), was used to generate initial perovskite geometries for DFT optimization (Lufaso and Woodward, 2001). The BVM Global Instability Index (GII), which SPuDS minimizes to predict the phase and magnitude of perovskite octahedral tilting, is strongly correlated with DFT energy in perovskite oxides (Morelock et al., 2021). This enables SPuDS to accurately predict perovskite ground state polymorph structures consistent with DFT, thereby substantially reducing the computational expense associated with high-throughput DFT investigations (Bare et al., 2021b). Initial geometries for DFT optimizations of the La 2 CoFeO 6 , GdLaCoFeO 6 , and Gd 2 CoFeO 6 perovskites were generated in the a-b+a-Glazer mode, as this is the ground state tilting mode most frequently predicted by DFT for experimentally observed ABO 3 perovskite oxides (Bare et al., 2021b). Atomic configurations for cation alloying on the B site were generated using rock salt site ordering, while configurations for alloying on the A sites were generated that minimize the Ewald sum.
The specific pseudopotentials and Hubbard +U parameters used for GGA+U DFT optimizations are compatible with the Materials Project (MP) database , which tabulates the structures and energies of >130,000 inorganic materials. Calculations were performed using the Vienna Ab initio Simulation Program (VASP 5.4.1) (Kresse and Hafner, 1993;Kresse and Hafner, 1994;Kresse and Furthmüller, 1996a;Kresse and Furthmüller, 1996b) with periodic boundary conditions utilizing projector augmented wave (PAW) pseudopotentials (Kresse and Joubert, 1999) and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional (Perdew et al., 1996). +U parameters of 3.32 and 5.3 were used for Co 3+ and Fe 3+ , respectively, consistent with pymatgen's MPRelaxSet . The electronic wave functions were expanded in a plane wave basis with an energy cutoff of 520 eV. The Brillouin zones were sampled during geometry optimizations using an automatically generated Γ-centered Monkhorst-Pack k-point mesh with a grid density of at least 1,000/(atoms/unit cell). Oxide-specific corrections to DFT total energies were included to maintain compatibility with the MP .
We explicitly considered the effects of magnetism by first performing two consecutive spin-polarized relaxations initialized in a high-spin ferromagnetic configuration with species-specific initial magnetic moments dictated by the default MP spin parameters. Then, magnetic sampling of the computed structures was performed for up to 20 different magnetic symmetries using pymatgen's MagneticStructureEnumerator . Finally, the internal coordinates of the DFT structures-with lattice vectors and initial magnetic moments fixed from previous optimizations-were optimized such that total energies were converged to within 10 -6 eV, and forces converged to within 0.01 eV/Å. The effects of spin configuration on the electronic density of states (DOS) were also explicitly described in this manner. Calculations to determine oxygen vacancy energies were performed using the aforementioned convergence criteria for all symmetrically unique vacancies at a defect concentration of C d 0.0833.
Synthesis of Gd 0.5 La 0.5 Co 0.5 Fe 0.5 O 3 Synthesis was conducted via a solid state synthesis route. For a 2 g scale reaction, stoichiometric amounts of gadolinium (III) oxide (0.7151 g), lanthanum (III) oxide (0.6427 g), cobalt (II, III) oxide (0.3167 g), and iron (III) oxide (0.3150 g) were ground by hand in an agate mortar and pestle for ∼10 min. The resulting powder mixture was calcined in air at 600°C (5°C/min, 12 h dwell), then sintered in air at 1,300°C (9°C/min, 12 h dwell) with intermediate grinding.

Characterization
Powder X-ray diffraction (XRD) was collected on Bruker D2 Phaser X-ray Diffractometer with Cu Kα radiation. Profile fitting of diffraction patterns were performed with GSAS-II (Toby and Von Dreele, 2013). Thermogravimetric analysis (TGA) was performed on a Netzsch STA 449 F1 Jupiter thermal analyzer.
Thermochemical cycling experiments were conducted under gas flow rates of 100 ml/min. For thermal reduction under Ar (100 ml/min), the sample (∼50 mg) was first heated to 1,250°C (10°C/min), held isothermally for 30 min (thermal reduction), cooled to 400°C (25°C/min), then held isothermally for 30 min. For re-oxidation, the gas was switched to a mixture of air (80 ml/ min) and Ar (20 ml/min), and the sample was heated to 1,100°C (10°C/min), isothermally held for 30 min, then cooled to 200°C (25°C/min). For repeated consecutive analyses, the sample was reweighed between runs. TGA baseline correction was performed with an empty crucible. All thermograms shown here are corrected.
Water splitting experiments were conducted in a stagnation flow reactor (SFR) equipped with a laser-based sample heater and a mass spectrometer; the experimental details are described in previous papers (Scheffe et al., 2011;Arifin et al., 2012;Scheffe et al., 2013). In brief, the powder sample was placed in a tube furnace to maintain the oxidation temperature, and then irradiated by an IR laser through an optical access window to achieve the desired reduction temperature. The amount of oxygen gas evolved during reduction under Ar and hydrogen gas evolved during oxidation under 40 vol% steam were measured using a mass spectrometer.
High temperature XRD (HT-XRD) was performed on a Scintag PAD X-ray diffractometer, equipped with a Buehler hot-stage with Pt/Rh heating strip and surround heater. The hot stage resides within a sealed chamber with an X-raytransparent beryllium window, and gas flow (either helium or air) was maintained at 200 ml/min. The data was analyzed using MDI Jade 8.2 software, and the plot was constructed with square root of intensity to easily observe low intensity peaks.

RESULTS AND DISCUSSION
To evaluate compound stability relative to decomposition, the DFT computed energies of the La 2 CoFeO 6 (L2CF), GdLaCoFeO 6 (GLCF), and Gd 2 CoFeO 6 (G2CF) perovskites optimized in the monoclinic, triclinic, and monoclinic space groups, respectively, (Supplementary Figure S1) were compared to the DFT computed energies of their potential decomposition products tabulated in the MP database. The energy of a material relative to the convex hull, E hull , quantifies a material's stability relative to its decomposition products. Monoclinic L2CF, triclinic GLCF, and monoclinic G2CF exhibit E hull of 0, 6.8, and 0 meV/atom, respectively, predicting that all three STCH candidates can be successfully synthesized as perovskites relative to their competing phases. The combination of the τ-predicted stabilities and E hull values of L2CF, GLCF, and G2CF suggests that Gd-La-Co-Fe-O is a promising compositional space for experimentally synthesizable perovskite oxides. We thus selected L2CF, GLCF, and G2CF from our high-throughput screening and computed additional STCH-relevant properties for these compounds, including DOS and charge neutral oxygen vacancy formation enthalpies (ΔH Ovac ).
The distributions of ΔH Ovac computed by DFT for the symmetrically unique sites of L2CF, GLCF, and G2CF are shown in Figure 1A. The ΔH Ovac distribution of GLCF is bounded by the ΔH Ovac distributions of L2CF (lower bound) and G2CF (upper bound). Herein, we use the DFT computed ΔH Ovac of CeO 2 (3.95 eV/atom) (Abanades and Flamant, 2006;, the gold standard STCH redox mediator (Muhich et al., 2016), as the upper bound for the STCH active range and a liberal lower bound of 2 eV/atom to account for uncertainty in DFT energetics (Naghavi et al., 2020). Materials approaching the minimum ΔH Ovac for STCH activity exhibit slow oxidation kinetics and/or degradation during redox cycling, Frontiers in Energy Research | www.frontiersin.org October 2021 | Volume 9 | Article 750600 whereas materials approaching the maximum ΔH Ovac for STCH activity suffer from reduced H 2 capacity (Muhich et al., 2016). GLCF exhibits both O vacancies with low and high ΔH Ovac , where vacancies with low ΔH Ovac participate in H 2 production and those with high ΔH Ovac are less likely to form under STCH operating conditions and therefore enable preservation of the perovskite lattice during redox cycling. This ΔH Ovac distribution predicts enhanced cyclability of GLCF relative to L2CF and increased H 2 production capacity relative to G2CF. Additionally, GLCF exhibits a large DOS effective mass, m e , that arises from its large concentration of accessible electronic states near the conduction band minimum (CBM). Lany showed that large m e corresponds with large electronic contributions to the entropy of reduction S red in STCH processes (Lany, 2018) that benefits STCH performance at high temperatures (Meredig and Wolverton, 2009). Figure 1B shows the DFT computed DOS for L2CF, GLCF, and G2CF. GLCF exhibits a larger DOS m e at 1000 K (m e 8.356) than both L2CF (m e 7.304) and G2CF (m e 3.209). The introduction of Gd into the nominally La-occupied A site of L2CF results in unoccupied states that lie closer to the CBM than those of La alone, which increases m e . However, complete substitution of Gd for La increases the splitting of the unoccupied Gd-f states that results in lowering m e relative to GLCF. Due to its favorable ΔH Ovac distribution and larger m e relative to L2CF and G2CF, GLCF was recommended for experimental synthesis and characterization of STCH performance.
The redox behavior of GLCF was examined through thermogravimetric analysis (TGA), which monitors the mass change with respect to temperature and atmosphere (Ar or air). In the case of these materials, the mass change corresponds to gain/loss of oxygen. The mass change was converted to the extent of reduction (δ), assuming Gd 0.5 La 0.5 Co 0.5 Fe 0.5 O 3-δ during redox cycling. The as-prepared powder was subjected to redox cycles, comprising a reduction step at 1,250°C under Ar followed by a re-oxidation step at 1,100°C under air, corresponding to 0.16 atm O 2 . The redox cycle was repeated twice with the two redox cycles showing similar behavior. The second redox cycle is shown in Figure 3. During thermal reduction, a shallow weight loss was observed at ∼750°C, followed by a steeper weight loss onset at ∼1,030°C. After the isotherm at 1,250°C, the extent of reduction (δ) was ∼0.13. During re-oxidation, the mass sharply increased when the gas was switched from Ar to air at 400°C, followed by a more gradual increase as the temperature was raised. The dip in the extent of reduction seen between 900 and 1,100°C is attributed to a phase change in the material, that is, formation of a phase at around 900°C that ejects oxygen to achieve stability. This change in δ is reproducible between cycles. The sample mass returns to the starting mass after the re-oxidation step.
The water splitting capability of GLCF was verified using a stagnation flow reactor (SFR). Figure 4 shows the water splitting results performed with a thermal reduction temperature (T R ) of 1,350°C for 330 s and re-oxidation temperatures (T O ) of 850 and 1,000°C for 1,200 s under 40 vol% H 2 O. Water splitting was observed under both sets of conditions, with the amount of H 2 produced increasing with each cycle. For T O 850°C ( Figure 4A), the amount of H 2 produced was 67, 90, and 101 μmol/g for each consecutive cycle. Moreover, significantly less O 2 was released in the second reduction cycle compared to the first cycle, likely due to the water splitting step being kinetically limited. This is also evident from the long tails of H 2 gas evolved during the re-oxidation steps. For T O 1,000°C ( Figure 4B), the H 2 capacity increased (127 and 141 μmol/g), evident of an improvement in the rate of water splitting at higher temperature. Similar to the T O 850°C experiment, more H 2 was produced as cycle-number increased. In an ideal situation, once steady state has been achieved, the amount of H 2 should be twice of the amount of O 2 produced. From the SFR experiments with repeated cycles, GLCF is evolving and approaching this ideal steady state, though additional studies are needed to understand this behavior.
While LaFeO 3 was previously reported to have negligible solar thermochemical H 2 O and CO 2 conversion behavior (Jiang et al., 2014;Chen et al., 2017), similar perovskites to GLCF, La 0.6 Sr 0.4 Co 0.2 Fe 0.8 O 3 and LaFe 0.75 Co 0.25 O 3 , were shown to be active for solar thermochemical CO 2 conversion (STCH activity is unknown) (Nair and Abanades, 2018), indicating that the mixing of Co and Fe may have contributed to the solar thermochemical conversion activity. For La 0.6 Sr 0.4 Co 0.2 Fe 0.8 O 3 and LaFe 0.75 Co 0.25 O 3 , however, CO production decreased substantially during subsequent cycles (Nair and Abanades, 2018). The substitution of Gd for La may have contributed to minimizing performance degradation as predicted by DFT calculations.
In terms of its water splitting ability, at T R 1,350°C and T O 850°C, GLCF produced more H 2 (101 μmol/g) than CeO 2 (50 μmol/g) (Barcellos et al., 2018). However, GLCF produced less H 2 compared to the previously studied perovskite materials BaCe 0.25 Mn 0.75 O 3 (BCM) and Sr 0.4 La 0.6 Mn 0.6 Al 0.4 O 3 (SLMA4664) Barcellos et al., 2018). BCM and SLMA4664 produced 140 μmol/g (T R 1,350°C, T O 850°C) and 307 μmol/g (T R 1,350°C, T O 1,000°C) of H 2 , respectively Barcellos et al., 2018). However, direct comparisons of performance reported for different conditions (temperature, atmosphere, and time) for various materials that have different optimized conditions for STCH can lead to incorrect conclusions about the H 2 production capabilities of candidate materials. In the present case, the water splitting experimental conditions implemented for GLCF have not yet been optimized. Nevertheless, computational screening greatly accelerated the discovery of GLCF as a water splitting material, which would otherwise have been experimentally time consuming due to the vast chemical space of perovskite materials.
To understand the phase stability in GLCF during redox cycling, HT-XRD patterns were collected first from room temperature to 1,250°C to 50°C under He for thermal reduction, then with an identical temperature profile under air for re-oxidation. The GLCF sample was redox-cycled (T R 1,250°C under Ar, T O 1,100°C under air) before collecting the HT-XRD patterns. The HT-XRD pattens are shown in Figure 5. During the HT-XRD experiment, GLCF appears to maintain its perovskite phase as the major phase, but the major peaks broaden and exhibit some peak splitting occurring during reduction. In addition to the changes in the major GLCF phase, during cool down at ∼1,200°C under He, additional peaks appear corresponding to Gd 2 O 3 phase. Due to the shift in peaks, it may have different chemistry (i.e. Gd 2-x La x O 3 ). The additional Gd 2 O 3 phase remains in the pattern until heated under air at ∼700-800°C during re-oxidation. The disappearance of the Gd 2 O 3 phase may correspond to the increase in the extent of reduction observed in the TG experiment at ∼900°C during reoxidation. Moreover, a spinel phase appears during cool down under air. Overall, during the redox cycle, the perovskite phase persists, agreeing with the stable STCH activity observed in the SFR results and TGA cycling tests.

CONCLUSION
Solar thermochemical hydrogen production (STCH) significantly contributes to the renewable energy portfolio.
However, materials discovery with high efficiency is needed, and computational screening can greatly accelerate this process. This was demonstrated with Gd 0.5 La 0.5 Co 0.5 Fe 0.5 O 3 (GLCF), a perovskite oxide, that was computationally determined first, then experimentally demonstrated. When compared to LaCo 0.5 Fe 0.5 O 3 , it was predicted that GLCF would have higher phase stability due to the incorporation of Gd on the La site. The synthesis of a single phase GLCF sample was achieved. With redox activity confirmed through thermogravimetric analysis, stable water splitting behavior over multiple cycles was also observed. When compared to previously reported La 0.6 Sr 0.4 Co 0.2 Fe 0.8 O 3 and LaFe 0.75 Co 0.25 O 3 , which suffered from performance degradation with subsequent cycles, the Gd substitution of La plays a significant role to maintain water splitting performance. Based on high temperature x-ray diffraction experiments, the GLCF perovskite phase persists which potentially contributes to the stable water splitting performance. This work demonstrates that computational materials screening can greatly accelerate the discovery of new water splitting materials. Through computational screening, potential water splitting materials were narrowed down to a promising STCH material from the vast chemical space of perovskite materials, and a unique strategy of incorporating a rare earth element that improved the stability was demonstrated.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author. October 2021 | Volume 9 | Article 750600 6