Computational screening of transition metal-doped phthalocyanine monolayers for oxygen evolution and reduction

Rationally designing efficient, low-cost and stable catalysts toward the oxygen evolution reaction (OER) and the oxygen reduction reaction (ORR) is of significant importance to the development of renewable energy technologies. In this work, we have systematically investigated a series of potentially efficient and stable single late transition metal atom doped phthalocyanines (TM@Pcs, TM = Mn, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ir and Pt) as single-atom catalysts (SACs) for applications toward the OER and ORR through a computational screening approach. Our calculations indicate that TM atoms can tightly bind with Pc monolayers with high diffusion energy barriers to prevent the isolated atoms from clustering. The interaction strength between intermediates and TM@Pc governs the catalytic activities for the OER and ORR. Among all the considered TM@Pc catalysts, Ir@Pc and Rh@Pc were found to be efficient OER electrocatalysts with overpotentials ηOER of 0.41 and 0.44 V, respectively, and for the ORR, Rh@Pc exhibits the lowest overpotential ηORR of 0.44 V followed by Ir@Pc (0.55 V), suggesting that Rh@Pc is an efficient bifunctional catalyst for both the OER and ORR. Moreover, it is worth noting that the Rh@Pc catalyst can remain stable against dissolution under the pH = 0 condition. Ab initio molecular dynamic calculations suggest that Rh@Pc could remain stable at 300 K. Our findings highlight a novel family of two-dimensional (2D) materials as efficient and stable SACs and offer a new strategy for catalyst design.


Introduction
Increasing energy demand and fast depletion of fossil fuels have led to search for alternative energy sources and efficient energy conversion technologies. [1][2][3] Sustainable and renewable energy generation technologies such as water splitting cells, fuel cells, and metal-air batteries 4,5 are regarded as promising approaches. These electrochemical technologies generally involve the OER and ORR, which have attracted much interest in the energy conversion research. The OER occurs on the anode side of an electrochemical water splitting cell, while as a reverse reaction of the OER, the ORR occurs on the cathode side of a fuel cell and a metal-air battery. Currently, Ru/Ir oxides 4,6 and Pt oxides as well as their alloys [7][8][9] are commonly used as the most efficient OER and ORR catalysts, respectively. However, these noble metals are rather expensive, and only a small number of surface sites can serve as catalytically active sites, which make their use rather inefficient. Therefore, the design and development of novel families of low-cost electrocatalysts whose catalytic activities are comparable or even higher than those of noble metal oxides are of signicant importance.
Recently, single-atom catalysts have attracted extensive attention, as metal atoms individually dispersed on supports can be promising for the maximum metal element utilization. 10 However, the isolated metal atoms are easy to aggregate to form clusters or nanoparticles due to their high surface energies. 11 To maintain the stability of SACs, the previous experimental and theoretical studies have demonstrated that isolated metal atoms can bind at appropriate defects or hollow sites of 2D supporting materials to avoid cluster aggregation. 12,13 Moreover, due to the connement of single atoms in the appropriate supports, the electronic properties of doped single metal atoms can be tuned with different transition metal elements which can enhance their catalytic activity. [14][15][16] In the last few years, single metal electrocatalysts for the OER/ORR have been the subjects of extensive studies. Especially, metal and nitrogen co-doped carbon (M-N x -C) materials have been demonstrated in many investigations to show OER/ORR performance with comparable activities to those of Ir/Pt oxides, 17-21 thus suggesting a promising way to replace noble metal oxides by M-N x -C catalysts. For example, our previous study 22 found that the OER overpotential roughly decreases with the increase of the coordination number of N-TM. Wang et al. 23 reported that for Co and N codoped graphene (CoN x -gra, x ¼ 1-4), both OER and ORR overpotentials roughly decrease with the increase of the doping concentration of nitrogen. The root of the M-N x -C materials used as electrocatalysts can be traced back to the discovery of the capability of M-N 4 macrocycle complexes (e.g. porphyrin) toward the ORR. 24,25 To date, a series of metal-doped phthalocyanine (MPc) monolayers, whose structure is similar to that of metal-doped porphyrins with one metal atom connecting to four nitrogen atoms, have been successfully synthesized in experiments with single metal atoms orderly and strongly anchored into the pores of the Pc. [26][27][28][29][30] The synthesis procedure of MPc is exible so that metal atoms can be replaced by other TM atoms. 31 More importantly, due to the large surface area, unique atomic structures, and intrinsic properties of dispersed metal sites, MPc monolayers and their derivatives have been predicted to be potential candidates for spintronics, 32-35 gas capture, 36 and energy conversion. [37][38][39][40][41][42][43][44][45] Using rst principles calculations, Zhao et al. 46 found that the 2D Cr-Pc monolayer exhibits high catalytic activity toward CO oxidation at room temperature. Wang et al. 47 theoretically reported that FePc shows high activity for the ORR in an acidic solution due to the stable binding of Fe atoms with Pc monolayers and the coordinative unsaturated state of the doped active center, which has been conrmed by recent experimental work. 48 All these previous studies suggest that the experimentally available TM@Pc monolayers with uniformly distributed single-atom active sites can be promising candidates for SACs.
In this work, the catalytic performance of 2D transition metal doped Pc monolayers (TM@Pc, TM ¼ Mn, Fe, Co, Ni, Cu, Ru, Rh, Pd, Ir, and Pt) as potential electrocatalysts toward the OER and ORR will be systematically screened with density functional theory calculations. This systematic study will be useful to help experimentalists to understand and select these transition metal elements, as these metals have exhibited OER and ORR activities in various other 2D materials. 17,49,50 The advances in DFT have made it possible to accurately describe catalytic reactions, 51 and the computational investigation on the OER/ ORR performance of TM@Pc monolayers can shed some light on developing low-cost, effective and stable electrocatalysts, or on future improvements of the systems. Our calculations demonstrate that Rh@Pc monolayers exhibit efficient catalytic activity toward both the OER and ORR with stability both thermodynamically and kinetically.

Computational methods
The spin-polarized density functional theory calculations were carried out with the Vienna Ab initio Simulation Package (VASP) code. 52,53 The projector augmented wave (PAW) pseudopotentials were used to describe the electron-ion interactions. 54 The Perdew-Burke-Ernzerhof (PBE) 55 functional of the generalized gradient approximation (GGA) 56 was used to describe the electron-correction interactions. The van der Waals (vdW) interactions were described using the Grimme's DFT-D3 correction method. 57 A plane-wave cutoff energy of 500 eV was adopted for all the computations to describe all atoms' valence electrons.
Geometry optimizations were performed until the atomic force became less than 10 À2 eVÅ À1 . The energy was converged to be less than 10 À5 eV. A vacuum space of 20Å was used to prevent the interaction between the periodic images. The Brillouin zone was sampled with 5 Â 5 Â 1 Monkhorst-Pack k-meshes. 58 Bader charge analysis was performed to evaluate the charge transfer process. 59 To investigate the stabilities of the screened out catalysts that possess excellent performance for the OER and ORR, the TM diffusion barriers were calculated by the climbing image nudged elastic band (CINEB) method. 60,61 Ab initio molecular dynamics (AIMD) calculations were also carried out to evaluate the corresponding dynamic stability. The algorithm of the Nose thermostat was used to calculate the canonical ensemble 62 at 300 K for 10 ps with a time step of 2 fs. The implicit solvent model was used to simulate the solvent environment throughout the whole process using VASPsol under the water conditions. 63 The details of OER and ORR calculations are provided in the ESI † as in our previous work. 17 The adsorption Gibbs free energy is dened as the following eqn (1): here, G adsorbent+catalyst , G catalyst , and G adsorbent are the Gibbs free energy of the adsorbent on the catalyst, the isolated catalyst, and the isolated adsorbent, respectively.  26,28,32,33,46 In the optimized TM@Pc monolayer structure, one TM atom binds with four inwardly projecting N atoms, forming four TM-N bonds. The binding energy (E b ) is dened as

Results and discussion
where E TM@Pc and E Pc are the total energies of the TM@Pc system and the Pc substrate, respectively. m TM is the chemical potential of the TM atom calculated from its bulk crystal. Since m TM is referenced with respect to its bulk metal, negative values of E b (Fig. 1b) indicate that the TM atoms in the Pc monolayers are stable against clustering. As shown in Fig. S1, † the strong hybridization between the 2p orbital of N and the d orbital of TM atoms demonstrates the chemical bonding of N and TM, which further explains the strong interaction between TM and Pc monolayers. The calculated bond length of TM-N, TM binding energies, and Bader charge transfer are listed in Table S1.   Fig. S4. † The data can be found in Table S2. † We can conclude that with the increase of the d-electron numbers of the TM atoms in the same row of the periodic table, the adsorption Gibbs free energies of intermediates decrease. This is also consistent with the position of 3 d as exhibited in Fig. S3. † Thus, there is a negative correlation between 3 d and Gibbs free energy of intermediates, at least when the TM atoms are in the same row of the periodic table. This phenomenon was also observed in previous experimental and theoretical studies. 17,18,71,72 Accordingly, the interaction strength can be modulated to an optimal value by tuning the TM doped on the Pc sheet toward OER and ORR performance.
As proposed by Nørskov et al., 73 the adsorption Gibbs free energies of three intermediates (HO*, O* and HOO*) on the TM@Pc catalyst govern the intrinsic catalytic activity toward the OER and ORR. Three descriptors, Gibbs free energies of adsorbed HO*, O* and HOO*, are used to evaluate the catalytic activity of an OER/ORR electrocatalyst. According to the Sabatier principle, 74 both too strong and too weak interaction strength between the intermediates and catalysts lead to adverse effects on the catalytic performance. Therefore, identifying promising electrocatalysts with moderate interaction of the reaction intermediates is one of our goals. As an ideal catalyst at the U ¼ 0 V condition which can occur at its thermodynamic limit, it requires that the free energy barriers between two adjacent intermediate states for all the mentioned four electron transfer steps (equations from Sa to Sd in the ESI †) should be the same, that is 4.92 eV/4 ¼ 1.23 eV. Thus, according to the equations from (S2a) to (S2e), † we can conclude that the adsorption Gibbs free energy of HO*, O* and HOO* should be 1.23, 2.46 and 3.69 eV, respectively. 17,71 Thus, both the OER and ORR can occur at their thermodynamic limit and the corresponding overpotential h will be zero. In reality, the Gibbs free energy differences between two adjacent intermediate states are not equal. The overpotential of the OER (h OER ) is determined by the maximum difference between the two adjacent Gibbs free energies, while the overpotential of the ORR (h ORR ) is determined by the minimum difference between the two adjacent Gibbs free energies. Obtaining the relationship among these three Gibbs free energies will simplify the analysis and the search for the optimal catalysis. 75,76 Here, the scaling relationship of DG HO* vs. DG HOO* for all the considered TM@Pc catalysts is shown in Fig. 2. We found that DG HOO* can be expressed as a function of DG HO* via equation DG HOO* ¼ 0.82DG HO* + 3.14, with a high coefficient of determination (R 2 ¼ 0.992). The slope close to unity in the correlated adsorption free energies of HO* vs. HOO* reects the fact that both HO* and HOO* have a single bond between the O atom and TM, which is similar to the cases of metal and metal oxide surfaces. 77,78 Thus, the overpotential (h OER and h ORR ) as a function of four variables (DG a , DG b , DG c and DG d ) can be reduced to two independent variables (another constraint is DG a + DG b + DG c + DG d ¼ 4.92 eV, the standard Gibbs free energy of H 2 O formation from O 2 and 2H 2 ). 73,78 As shown in the following: 82DG HO* + 3.14) À DG O* and DG d ¼ 4.92 À DG HOO* ¼ 4.92 À (0.82DG HO* + 3.14). Thus, knowing only two descriptors, DG HO* and DG O* À DG HO* , is sufficient for us to describe the catalytic performance of a system toward the OER and ORR. In Fig. 3a, we plot the two-dimensional volcano to exhibit the OER activity trends through h OER as a function of two independent descriptors DG HO* and DG O* À DG HO* . The blue region of the plot shows the highest activity area with h OER reaching a minimum value of 0.21 V under the optimum condition (DG a ¼ DG b ¼ DG c ¼ 1.44 eV). Note that the minimum value of the overpotential is not zero since the relationship DG HOO* ¼ 0.82DG HO* + 3.14 excludes the ideal case of DG HO* ¼ 1.23 eV and DG HOO* ¼ 3.69 eV. Different catalysts fall on different points in the volcano plot of Fig. 3. Based on this, for the OER, Ir@Pc is the best catalyst (h OER ¼ 0.41 V) followed by Rh@Pc (h OER ¼ 0.44 V). The free energy diagrams of all the intermediate states of Rh@Pc and Ir@Pc toward the OER are shown in Fig. 4a and b at U ¼ 0 V. Notably, the calculated h OER values of Rh@Pc and Ir@Pc are comparable or even lower than those of the current best catalysts RuO 2 (h OER ¼ 0.42 V) 78 and IrO 2 (h OER ¼ 0.52 V). 77 Co@Pc also shows good activity (h OER ¼ 0.50 V) for the OER with the O* formation being the ratedetermining step. The ORR is the reverse reaction of the OER. In Fig. 3b the calculated h values for the ORR on various TM@Pc catalysts were compared. Under the optimal condition (ÀDG a ¼ ÀDG d ¼ 0.98 eV), the theoretical h ORR is found to be as low as 0.25 V. From the volcano plot in Fig. 3b, the best TM@Pc catalyst for the ORR is found to be Rh@Pc with an h ORR value of 0.44 V and the rate-determining step is the reduction of O 2 to HOO* (Fig. 4b), followed by Ir@Pc (h ORR ¼ 0.55 V). This thermodynamic limiting overpotential is even lower than that of Pt (111) (h ORR ¼ 0.48 V). 79 Fe@Pc and Co@Pc also exhibit good catalytic activities (h ORR ¼ 0.58 and 0.58 V) for the ORR with the rate-determining steps of reduction of O* to HO* and HOO* formation, respectively. Here, it should be noted that Rh@Pc can be used as the efficient bifunctional electrocatalyst for both the OER and ORR with an h OER value of 0.44 V and an h ORR value of 0.44 V, respectively. To evaluate the dynamic stability of this promising catalyst, we have performed AIMD simulations under 300 K condition for 10 ps (Fig. S5 †) for Rh@Pc. It can be seen that the energies oscillate near the equilibrium state while the structure remains unchanged, suggesting the kinetic stability of the Rh@Pc catalyst.
In another aspect, we have also calculated the reaction free energy (DG diss ) to evaluate the stability of the doped metal centers against dissolution due to the proton attack of the active region using the equation: 75 TM@Pc + nH + / nH@Pc + TM n+ . Where n refers to the oxidation state for the TM atom, nH@Pc refers to the Pc monolayer with the TM vacancy adsorbed by n number of hydrogen atoms (Fig. S6 †). The dissolution energy can be calculated as: DG diss ¼ G (nH@Pc) + G (TM n+ ) À G (TM@Pc) À nG (H + ) . Here, G (nH@Pc) and G (TM@Pc) can be calculated directly, and DG (H + ) ¼ 0.5 Â G (H 2 ) À ln 10 Â kT Â pH ¼ 0.5 Â [E (H 2 ) + ZPE (H 2 ) À TS (H 2 ) ] À ln 10 Â kT Â pH. Here, TS (H 2 ) ¼ 0.41 eV for the H 2 gas phase at 298 K, and ZPE (H 2 ) ¼ 0.26 eV. 17 As for G (TM n+ ) , we take the experimental ion formation of TM n+ , which is dened as: DG (TM n+ ) ¼ Fig. 2 The scaling relationship between the adsorption Gibbs free energies of HO* and HOO* for all the considered TM@Pc systems.  Table. S3. † Using this approach, the dissolution energy under pH ¼ 0 condition, DG diss (0), is calculated and shown in Fig. 5. We can conclude that Rh@Pc, Pd@Pc, Ir@Pc, and Pt@Pc catalysts are stable against dissolution under the pH ¼ 0 condition. For the other TM, their DG diss (0) values are less than zero, which means they will be unstable against dissolution. However, as the pH value increases, DG diss (pH) will also increase, thus, there will be a critical pH value, above which the TM@Pc catalyst will be stable. The critical value can be calculated by pH min ¼ ÀDG diss (0)/(n Â 0.0591) (Table S3 †). Thus, when the system is sufficiently alkaline, it will always be stable against dissolution, and Rh@Pc, Pd@Pc, Ir@Pc, and Pt@Pc are stable even in a very acidic environment.

Conclusion
To summarize, we have systematically screened a series of single TM atom doped Pcs as potentially efficient and stable SAC bifunctional electrocatalysts for both OER and ORR catalytic processes using a computational screening approach. Based on the computations of binding of TM atoms with Pc monolayers, we found that all the considered TM atoms exhibit a strong interaction with Pc monolayers as potentially stable SACs with high diffusion energy barriers. With the increase of the d-electron number of the doped TM atom on Pc monolayers that leads to the lower d-band center, the interaction strength between intermediates and the active doped TM atoms will decrease, which allows us to select the optimal TM@Pc catalyst toward the OER/ORR by tuning the doped TM element. According to the volcano plots of the OER and ORR, among all the studied TM@Pc catalysts, the best catalyst for the OER is Ir@Pc with an h OER of 0.41 V followed by Rh@Pc with h OER ¼ 0.44 V, and for the ORR, the best catalyst is Rh@Pc with an h ORR of 0.44 V followed by Ir@Pc (h ORR ¼ 0.55 V). It should be noted that both Rh@Pc and Ir@Pc can remain stable against dissolution under all pH conditions. This study highlights a potentially efficient new class of SACs based on the Pc monolayers toward the OER and ORR.

Conflicts of interest
There are no conicts to declare.