High-Throughput Computational-Experimental Screening Protocol for the Discovery of Bimetallic Catalysts

For decades of catalysis research, the d-band center theory that correlates the d-band center and the adsorbate binding energy has successfully enabled the accelerated discovery of novel catalyst materials. Recent studies indicate that, on top of the d-band center value, the full consideration of the d-band shapes describing higher moments of the d-band as well as sp-band properties can help better capturing surface reactivity. However, the density-of-states (DOS) patterns themselves have never been used as a descriptor in combined computational-experimental studies. Here, we propose the full DOS patterns as a key descriptor in high-throughput screening protocols, and prove its effectiveness. For the hydrogen peroxide (H2O2) synthesis as our demo catalytic reaction, the present study focuses on discovering bimetallic catalysts that can replace the prototypic palladium (Pd) one. Through a series of screening processes based on DOS pattern similarities (evaluated using first-principles calculations) and synthetic feasibility, 9 candidates are finally proposed out of 4,350 bimetallic alloy, which then are expected to have a catalytic performance comparable to that of Pd. The subsequent experimental tests demonstrate that 4 bimetallic catalysts (Ni61Pt39, Au51Pd49, Pt52Pd48, Pd52Ni48) indeed exhibit the catalytic properties comparable to those of Pd. Moreover, we discovered a novel bimetallic (Ni-Pt) catalyst, which has not yet been reported for H2O2 direct synthesis. In particular, Ni61Pt39 outperforms the prototypical Pd catalyst for the chemical reaction and exhibits a 9.5-fold enhancement in cost-normalized productivity. This protocol provides a new opportunity for the catalyst discovery for the replacement or reduction in use of the platinum-group metals.


Introduction
One of the important roles of computer simulations in materials science is to predict and/or design novel materials prior to experiments. In particular, first-principles calculations using density functional theory (DFT) have played a vital role in the field of catalysis. [1][2][3] The computational framework enables the prediction or understanding of reaction pathways on catalyst surfaces, which is undoubtedly useful for experimentalists. However, the accurate estimation of some crucial properties such as the reaction barriers on catalytic surfaces is extremely time-consuming in first-principles calculation scheme. In this regard, obtaining a full understanding of the catalytic reaction mechanism via first-principles calculations is often considered rather inefficient. In our experiences, though computing-power-dependent, it requires time costs of several months even for computing a full energetics picture of one metallic catalyst system. The efficiency has therefore been questioned for the first-principlescalculation-guided search of novel catalysts, because carrying out experimental tests without the aid of computer simulations could often be even faster.
To maximize the efficiency of novel catalyst discovery using the combined computational and experimental study, the descriptors or features that bridge these two area should be wisely chosen. To be specific, a simple but physically reasonable descriptor or feature enabling the representation of catalytic properties is critical. A representative example of such a descriptor is the d-band center theory that correlates between the d band center, 4 i.e., the average energy of electronic d states projected onto a surface atom of the catalyst, and the gas adsorption energy, which has been widely accepted in the field of metallic catalysis.
Coupled with the volcano relationship between catalytic activity and adsorption energy, 5 it is theoretically possible to predict catalytic activity from the d-band center without exploration of the full reaction mechanism. Afterward, along with the d-band center theory, consideration of the d-band shapes describing higher moments of the d-band as well as the sp-band properties describing local Pauli electronegativity were also suggested to capture the surface reactivity of transition metal alloys. [6][7][8] These studies indicate that the electronic density of states (DOS) patterns themselves projected onto surface atoms of catalysts can serve as an improved descriptor for the development of novel metal catalysts because they include more comprehensive information of not only d states but also sp states (both their values and shapes). Nonetheless, the full DOS patterns themselves have never been used as a descriptor in combined computational-experimental screening process.
Since electronic structures are key to determining the physical/chemical properties of materials, materials with similar electronic structures tend to exhibit the similar properties.
Indeed, a DFT study revealed that a solid-solution Ir 50 Au 50 alloy can exhibit catalytic properties for H 2 dissociation reaction with a comparable activity of Pt, and this result is attributed to their similar electronic structures (DOS patterns) between the Ir 50 Au 50 alloy and Pt. 9 In another example of Rh 50 Ag 50 versus Pd for hydrogen storage, their similar electronic structures were also experimentally verified. 10 The similar electronic structural characteristics enables the Rh 50 Ag 50 alloy to exhibit hydrogen storage properties as superior as Pd, 11 although neither Rh nor Ag alone shows hydrogen storage properties. 9 Thus, bimetallic alloys with a DOS pattern similar to that of Pd, which is well-known as a representative metal catalyst, would exhibit catalytic performance similar to that of Pd, which is a good hypothesis for the high-throughput development of bimetallic catalysts.
Here, we use the full DOS pattern as a key descriptor in the high-throughput computationalexperimental screening protocols. As an exemplary demonstration in the high-throughput search of bimetallic catalysts, we select palladium (Pd), the prototypical catalyst for hydrogen peroxide (H 2 O 2 ) synthesis from hydrogen (H 2 ) and oxygen (O 2 ) gases, 12-14 as our reference material. Using DFT calculations, we first screened 4,350 crystal structures of bimetallic alloys to investigate their thermodynamic stabilities (only closest-packed surfaces considered). Next, the similarities in the DOS patterns between the alloys and Pd were quantified, and subsequently their synthetic feasibility was also evaluated. Through these processes, we finally proposed 9 candidates with the high DOS similarities with Pd, which then are expected to show catalytic properties comparable to those of Pd. These 9 screened alloys were experimentally synthesized and tested for H 2  found much to outperform the prototypic Pd with a 9.5-fold enhancements in cost-normalized productivity owing to the high content of inexpensive Ni.

Results
High-throughput computational screening for novel bimetallic catalysts There are three screening steps to identify catalyst candidates from 4,350 candidates. The top (cyan color) section indicates the first step showing 30 elements and 10 crystal structures to collect binary alloys with an 1:1 composition. The middle (magenta color) and bottom (violet color) sections indicate the second and the last steps relevant to the thermodynamic stability and the DOS similarity, respectively. Fig. 1 shows a schematic diagram of our high-throughput screening protocol for the discovery of bimetallic catalysts. Based on 30 transition metals in periods IV, V, and VI, we considered a total of 435 ( 30 C 2 ) binary systems with a 1:1 (50:50) composition. For each alloy combination, 10 ordered phases that are available for the 1:1 composition were investigated, namely, B1, B2, B3, B4, B11, B19, B27, B33, L1 0 , and L1 1 , leading to a screening of 4,350 (435×10) crystal structures. Using DFT calculations, we calculated the formation energy (∆E f ) of each phase and determined which crystal structure is most stable in the bulk phase. Here, a negative formation energy implies that the phase is thermodynamically favorable (or miscible element combination), and vice versa. Because nonequilibrium alloyed phases can be stabilized by the nanosize effect, although the elemental combinations are immiscible in the bulk phases, 11,15,16 a margin of ∆E f < 0.1 eV was considered when screening the thermodynamic stabilities of the bimetallic systems. Eventually, of a total of 435 binary systems, 249 alloys were filtered at the thermodynamic screening step.
For the 249 thermodynamically screened alloys, we calculated the DOS pattern projected on the close-packed surface for each structure and compared it with that of the Pd(111) surface.
To quantitatively compare the similarity of two DOS patterns, we defined the following measure: Here, due to the rarity of Hf and Ta, the CoHf and CoTa alloys were excluded from the candidate list. Additionally, we excluded NiCu and CrRh from the list because the large reduction potential difference between the elements in the binary system leads to unfavorable formation of the alloyed structures in which two elements are homogeneously mixed via a wet-chemical method for nanoparticle (NP) synthesis. 17 Figure 2. Map of the thermodynamic stabilities and DOS similarities of bimetallic alloy systems. The diagonal components denote the most thermodynamically stable crystal structures (e.g., FCC, BCC, HCP, and SC) of pure metals. The occupied components under the diagonal denote the most thermodynamically stable crystal structures (e.g., B2, B11, B19, B27, B33, L1 0 , and L1 1 ) of bimetallic alloys. The background colors above the diagonal describe the formation energies, and the occupied components above the diagonal denote the DOS similarity of each alloy compared with Pd (111).

Experimental verification for the DFT-guided bimetallic catalysts
For the remaining 13 candidates, we tried to synthesize their alloyed nanoparticles (NPs) and confirm whether their crystal structures are identical to those (Fig. 2) predicted by DFT calculations. For NP synthesis, we used a butyllithium reduction method. 18 Details of the NP synthesis are found in the Methods section. Because it was assumed in Fig. 2  CoFe, NiIr, NiMo, and NiNb were not alloyed because the reducing power of butyllithium used as a reducing agent was likely insufficient. To confirm whether the alloyed NPs have the same crystal structures as predicted by DFT, X-ray powder diffraction (XRD) patterns were investigated, as shown in Fig. 3. All of the samples showed alloyed diffraction patterns without signals from phases of the pure elements that compose each alloy. We also simulated XRD patterns of the alloys on the basis of the DFT crystal structures and found that the simulated patterns are matched well with the experimental patterns. This clearly reveals that we successfully synthesized alloyed NPs with the crystal structures predicted by DFT calculations for 9 bimetallic systems. The lattice parameters for the alloyed samples are presented in Table. 1 in Supporting Information. We further characterized the NP samples by a high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM)   although the system has been well known as an oxygen reduction reaction (ORR) catalyst in fuel cells. [27][28][29] As the price of Pd rapidly increases, 30 our Pd-free Ni-Pt alloy can serve as a more cost-effective option. To identify the cost-effectiveness of the Ni-Pt catalyst, we list the cost-normalized productivity (CNP; unit, mmol H2O2 ·$metal -1 ·h -1 ) 31 in  Table 1.

On the other hand, Wilson and Flaherty proposed another mechanism for H 2 O 2 formation on
Pd NPs called a heterolytic reaction mechanism, 35 similar to the two-electron oxygen reduction reaction (ORR), in which an electron is generated during the oxidation process of H 2 on the NP surfaces. We also investigated the mechanism on Ni 50 Pt 50 (111) and Pd(111) by DFT calculations where the computational hydrogen electrode method was used along with consideration of implicit water solvents (Fig. 4b). On both Pd (111) 25 Pt-Pd, 26 and Pt-Ag 36 bimetallic catalysts have been reported. The successful story of our protocol basically results from the validity of our descriptor, the similarity of electronic DOS patterns, for catalyst screening. In this regard, this protocol will immediately provide new opportunities for catalyst discovery for the replacement or reduction in use of the platinum-group metals and then readily apply them into other chemical reactions.

Methods
Computational details for screening bimetallic catalysts. We generated 10 ordered phases for each bimetallic alloy with a 1:1 composition. Specifically, the crystal structure prototypes were obtained by NaCl (B1), CsCl (B2), zinc blende (B3), wurtzite (B4), TiCu (B11), CdAu (B19), FeB (B27), CrB (B33), AuCu (L1 0 ), and CuPt (L1 1 ) compounds from the Material Project. 37 Here, the initial lattice parameters were replaced by values commensurate with covalent radii of two elements in the target compounds. By means of the DFT optimization process, we collected the most stable bulk phases of 4,350 binary alloys. Furthermore, the formation energies (ΔE f ) of all the bulk phases were calculated as follows: To investigate thermodynamic stabilities in bulk phases and the DOS patterns of surface atoms, spin-polarized DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP) 39 with projector-augmented-wave (PAW) pseudopotentials 40 and the revised Perdew-Burke-Ernzerhof (RPBE) exchange-correlation functional. 41 Then, we used a plane-wave kinetic energy cutoff of 400 eV and Monkhorst-Pack k-point meshes of 8×8×8 for the bulk systems and 4×4×1 for the slab systems. The bulk of the 1×1×1 unit cell and the slab of the 2×2 unit cell with 4 layers were used. In the slab systems, a vacuum spacing of 15 Å was used to avoid interactions between slabs. All relaxations were performed until the energy change was less than 1 × 10 −6 eV/cell and the forces on the individual atoms were less than 0.01 eV Å -1 .

H 2 O 2 formation mechanism.
The energy diagrams for H 2 O 2 production on the Ni 50 Pt 50 (111) surface and Pd(111) surfaces were also calculated by the DFT method. The simulations were performed using the VASP 39 using the PAW 40 pseudopotentials to describe the potential from the ionic core. For the exchange and correlation terms, we employed the RPBE functional 41 with Grimme's D3 dispersion correction. 42 An energy cut-off of 400 eV and Monkhorst-Pack k-point meshes of 3×3×1 for the slab calculations were used. A large vacuum spacing of 15 Å was used to prevent interslab interactions. The top two layers and the reactant molecules were optimized until the energy change was less than 1 × 10 −4 eV/cell and the force on each atom was less than 0.02 eV Å −1 . In the case of the Langmuir-Hinshelwood mechanism, the climbing image nudged elastic band (CI-NEB) 43 method was employed to investigate the transition states for the main reaction. For the electron-proton transfer mechanism, the computational hydrogen electrode method 44 with U = 0 V in addition to the water solvation energy calculated with the implicit solvent model implemented within the VASP sol package 45 was used to construct the energy diagram on the following reactions: Catalyst synthesis. All alloy metal catalysts and Pd catalysts were synthesized using the butyllithium reduction method. To synthesize the bimetallic XY catalysts, 0.017 mmol of X precursor and Y precursor (details are described in Table 3 of Supporting Information) were dissolved in 10 mL of dioctyl ether with 2 mL of oleylamine at 50 °C. The solution was then injected into a butyllithium solution containing 15 mL of dioctyl ether and 1.2 mL of 2.0 M butyllithium in cyclohexane at 50 °C. The colloids were stirred (500 rpm) for 20 min and then heated to 120 °C for 1.5 h in an Ar atmosphere. The reaction mixture was further heated to 260 °C for 1 h. The mixture was cooled to room temperature, and 1.25 mL of trioctylphosphine was injected to protect the colloids. The resulting NPs were washed with ethanol 3 times.
Catalyst characterization. Transmission electron microscopy (TEM) images and energydispersive X-ray spectra (EDS) were acquired using a transmission electron microscope (FEI Talos F200X) equipped with scanning transmission electron microscopy-energy dispersive X-ray spectroscopy (Bruker Super-X EDS system). The crystal structure was examined by X-ray diffraction (Rigaku Dmax 2500) with Cu Kα radiation (λ = 1.5406 Å). The samples for the TEM images and X-ray diffraction patterns were purified by centrifugation three times to remove the surfactants and/or excess reactants. Then, the NPs were dispersed in ethanol. The resulting solution was dropped on a copper grid coated with an amorphous carbon film for the TEM images and a soda-lime-silica glass for the X-ray diffraction pattern. The actual compositions of each alloy catalyst were measured via inductively coupled plasma atomic emission spectrometry (ICP-OES) using iCAP6500 Duo (Thermo).
Catalyst performance measurements.