Adsorption Site Screening on a PGM-Free Electrocatalyst: Insights from Grand Canonical Density Functional Theory

: Platinum-group-free Fe-N-C catalysts have been shown to be a promising class of electrocatalysts for the oxygen reduction reaction (ORR) in acidic media, with the Fe atom being regarded as the site most responsible for ORR activity. Despite previous density functional theory (DFT) modeling of this class of catalysts, the effect of the self-consistent applied potential and relevance of competing adsorbates to the ORR are still not well understood. In this work, we used grand canonical DFT to calculate the adsorption thermodynamics of H, O, OH, O 2 , H 2 O, NO, SO 4 , HSO 4 , and ClO 4 ligands to symmetrically unique Fe, C, and N sites for an FeN 4 moiety hosted within a graphene monolayer (FeN 4 @G). We find that the use of the applied potential within the grand canonical ensemble (GCE) can significantly alter the adsorption energies of these ligands in comparison to their canonical adsorption energies for all adsorption sites. Under a GCE applied potential, O, OH, H 2 O, NO, SO 4 , and HSO 4 can all have a lower adsorption energy to Fe than O 2 , suggesting that other ligands may be persistently bound to Fe atoms represented by this active site model during the ORR. Lateral spectator OH ligands can adsorb near the Fe site, and these ligands can modulate the grand canonical adsorption energies by approximately 0.1 − 0.4 eV. In addition, the use of a more oxidative applied potential affects the adsorption energies of ligands qualitatively differently, with the adsorption energies of H, O, O 2 , and NO being destabilized, while the adsorption energies of the other ligands are stabilized. A term-by-term comparison of the grand canonical and canonical adsorption energies shows that despite both the clean catalyst and adsorbed ligand systems oxidizing to similar (significant) degrees under an applied potential, this effect only partially cancels when calculating grand canonical adsorption energies. Additionally, we show that the choice of the molecular referencing scheme is important, and multiple schemes using both charged and neutral species with different implicit solvation models are compared.


INTRODUCTION
−10 Fundamental insights into the ORR mechanism on these catalysts are hampered, however, by the difficulty of modeling or observing reaction site mechanisms and surface structure under reaction conditions. 3his problem is compounded by the wide variety of experimental parameters that can potentially impact the mechanism, such as active site structure, pH, applied potential, and electrolyte concentration.
One promising PGM-free catalyst currently being studied involves planar Fe x N y sites substituted into graphene monolayers (Fe x N y @G) for two C atoms. 1 In this catalyst, it has been proposed that the Fe atoms bound within the Fe x N y moieties are the catalytically active sites for the ORR and that the reactants and intermediates in the ORR bind to Fe during reaction.Broadly, the Fe x N y sites can be either substituted into an extended graphene monolayer (a "bulk"-hosted site) or near/at the edge of the monolayer (an edge-hosted site). 1,10dditional defects in the graphene such as C defects or N substitutions can further modify the exact structure of an active site.7][8][9]11 The variety of possible Fe x N y @G sites, competing solvated species, and layered catalyst structure under reaction conditions (commonly: water solvent, U = 0.00 to 1.23 V vs SHE; 298 K; 1 atm; 0.5 M H 2 SO 4 [pH = 0]) greatly complicates both in situ monitoring of the catalyst during reaction and computational modeling. Foexample, other species, such as H + , OH − , H 2 O, or S-containing ligands (from H 2 SO 4 or Nafion), could all compete for adsorption onto Fe or other locations on the catalyst prior to O 2 adsorption.These species could also all interact sterically with a second adsorbing molecule or modify the electronic structure of the electrocatalyst via co-adsorption on the catalyst.Although this class of catalysts has been the subject of many recent experimental and computational studies, the inclusion of self-consistent applied bias in density functional theory (DFT) calculations has only had limited investigation, with the primary focus being the thermodynamics of the ORR intermediates.4,5,8,9,12 In this work, we modeled the FeN 4 @G electrocatalyst under solvation and experimentally relevant self-consistent applied potentials within the grand-canonical density functional theory (GCDFT) formalism, which has been shown to accurately capture the effects of potential in a variety of electrochemical systems.13−15 We investigated how H, O, OH, O 2 , H 2 O, NO, SO 4 , HSO 4 , and ClO 4 ligands adsorb to different Fe, N, and C sites on the FeN 4 @G catalyst to determine the thermodynamically favored species at each chosen applied bias.These ligands are commonly present either from the solution phase with sulfuric/perchloric acid (H, H 2 O, SO 4 , HSO 4 , and ClO 4 ), ORR mechanism (O 2 , O, and OH), or in experiments used to probe for the presence of Fe (NO).16 We explain the origins of shifts in adsorption energy with potential, quantify the role of the molecular referencing scheme for charged and neutral species in different implicit solvation models, and predict the impact of lateral spectator ligands on modifying adsorption to Fe.

METHODS
All joint DFT calculations were performed with the JDFTx code using a 50 Ryd planewave cutoff, a 3 × 3 × 1 k-point mesh, spin-polarization, the GBRV ultrasoft pseudopotentials, and Perdew−Burke−Ernzerhof (PBE) generalized gradient approximation (GGA) functional with Grimme's D3 van der Waals correction. 17−20 Solvation effects were included via the CANDLE implicit solvation model with the concentration of the electrolyte set to 0.5 M. 21 All potentials were referenced to the calibrated CANDLE point of zero charge. 21Coulombtruncation along the axis perpendicular to the graphene sheet was used. 22Calculations including potential were carried out within the grand canonical formalism by setting the applied potential to either 0 V vs SHE or 1 V vs SHE. 13 Because the model is a 2-dimensional material, the lattice of the FeN 4 @G catalyst in the planar lattice directions was relaxed under each applied potential condition.All adsorption calculations under each potential condition used the corresponding relaxed lattice.We note that the maximum changes in lattice vector lengths and angles among the three potential conditions were 0.02 Å and 0.1°, indicating that the lattice vectors of this model are not sensitive to potential choice in this range.The symmetrically unique adsorption sites on the FeN 4 @G catalyst were determined with Pymatgen. 23,24−27 The molecular reference states for all adsorbates were relaxed as molecules in a box with a 1 × 1 × 1 k-point mesh and 3-dimensional Coulomb truncation and otherwise the same settings as described above.The molecular reference states used for each adsorbate and ensemble are listed in Table 1.We emphasize that the requirement to maintain constant charge when calculating adsorption energies means that canonical reference schemes do not completely capture both the surface charging behavior and presence of charged reference species seamlessly described by the GCE.For our GCE calculations, H The energies for all charged molecules were calculated by specifying a charge state for these molecules and applying the linear free energy approximation, −μN, to calculate their energies at a specific applied potential because they do not interact with the applied potential via the catalyst surface when in bulk solvent.The CE molecular reference states that use the energy of 1 2 H 2 as the energy of a proton/electron pair are shifted by a potential-dependent (vs SHE) post-processing correction −eUN e , where N e is the number of proton/electron pairs subtracted.Because ORR experiments are commonly run at pH = 0 for this catalyst, the energy of H 3 O + and OH − was shifted by −kT ln (10)pH (and thus 0 for these calculations) and −kT ln (10)pOH (−0.83 eV at 298 K), respectively, which accounts for the concentration dependence of the entropy for these molecules. 1,28−33 We find that this energy difference is similar among functionals and that the PBE functional likely benefits from a cancellation of errors in predicting water equilibrium (Table S2).The oxidation of a graphene monolayer under the GCE applied potential was found to be similar using both the PBE and TPSS functionals.
The CANDLE model has been shown to perform well for a wide variety of systems, and implicit solvation models in general are a computationally expedient electrochemical modeling approach because they bypass the requirement for configurational sampling of explicit water molecules. 21,34urther development of these approaches is still needed however to accurately describe solvation, electric field, and hydrogen bonding effects across all types of electrochemical systems.Different implicit solvation models may describe the physics of these systems differently.To quantify this, we carried out similar adsorption energy calculations for Fe adsorption using the non-empirical spherically averaged liquid susceptibility ansatz (SaLSA) implicit solvation model but The Journal of Physical Chemistry C otherwise the same computational settings. 34The results for this study are shown in Table S3.The effect of the solvation model on all molecular reference states used in this work is shown in Table S4.

RESULTS AND DISCUSSION
3.1.FeN 4 @G Active Site Model.The planar bulk-hosted FeN 4 @G active site model was chosen to provide an initial screening due to its relative simplicity as compared to multi-Fe sites or edge-hosted sites and because it has been observed experimentally.All symmetrically unique adsorption sites for the FeN 4 @G electrocatalyst were determined, and the subset shown in Figure 1 was chosen for our screening because it includes Fe, N, and C atoms in several unique local bonding environments.The adsorption sites are denoted here as being either atop sites (located above 1 atom in FeN 4 @G), bridge sites (located between 2 atoms in FeN 4 @G), or hollow sites (located in the middle of a Fe/N/C ring in FeN 4 @G).Each adsorbate was initialized at each adsorption site shown in Figure 1 except for the H atom, which was only initialized at atop sites due to it only forming single bonds.OH was initialized in an O down configuration, H 2 O in a flat configuration, O 2 in a bidentate configuration, NO in an N down configuration, while HSO 4 , SO 4 , and ClO 4 were centered over each site.O 2 relaxed to both a monodentate and bidentate configuration on Fe during the screening however, thus accounting for both major binding modes.The effect of the applied potential was quantified by comparing canonical ensemble (CE) calculations (calculations that did not depend on applied potential during energy minimization but instead only via a post-processing correction for each electron in the chosen molecular reference) with implicit solvation to grand canonical ensemble (GCE) calculations with implicit solvation at an applied potential of either 0 V vs SHE or 1 V vs SHE (see Section 2 for details).0 and 1 V vs SHE were chosen as the studied potentials in order to span most of the potential range commonly used in ORR experiments.The charge-asymmetric nonlocally determined local-electric (CANDLE) implicit solvation method was used due to its accuracy in describing a wide variety of neutral molecules, cations, and anions. 21igure S1 shows the structure of the minimum energy relaxed adsorbate ligands on Fe under an applied potential of 0 V vs SHE.The minimum energy structures for the other potential conditions were qualitatively similar with the primary difference being changes in adsorbate bond lengths of typically 0.2 Å or less.
Without a potential within the CE, this clean FeN 4 @G model must retain a net charge state of 0 (Table 2).To realize this, the 54 C atoms all slightly oxidize to a total charge state of 0.60, the 4 N atoms reduce to a total charge state of −1.24, and the Fe atom oxidizes to a charge state of 0.64.The Fe atom has a net spin moment of approximately 2 and is the only atom with a significant net moment.Within the GCE, however, the net charge state of this clean FeN 4 @G model is 0.85 and 2.22 under 0 V vs SHE and 1 V vs SHE, respectively (Table 2).This result immediately demonstrates the impact of surface charging on this catalyst and implies that the point of zero charge (PZC) for this catalyst occurs at negative potentials vs SHE.Additionally, Table 2 shows that the change in the total system charge state can be primarily attributed to each C atom slightly oxidizing under 1 V vs SHE, resulting in the total C charge state increasing from 1.38 to 2.52.The slight oxidation of all C atoms under 0 and 1 V vs SHE is consistent with the behavior of a same-sized graphene monolayer with the FeN 4  Values in parentheses are the average charge state or spin moment for each type of atom.The 0 and 1 V vs SHE columns were calculated within the grand canonical ensemble.

The Journal of Physical Chemistry C
site replaced by C atoms, which increased its net charge state from 0.44 under 0 V vs SHE to 1.31 under 1 V vs SHE.The spin moments for all atoms are relatively insensitive to potential, although we note that the Fe spin moment does slightly increase as the potential increases.
The unequal change in C, N, and Fe charge states described above is reflected in the projected density of states (PDOS) for this catalyst under each potential condition (Figure 2).Without a potential, the highest energy valence states have Fe, N, and C character.Under GCE 0 V vs SHE and 1 V vs SHE, the catalyst is oxidized relative to CE calculation, with the highest energy valence states as calculated within the CE shifting above the Fermi energy (E F ). Notably, the states do not rigidly shift above the E F such that all atoms oxidize in proportion to their respective number of states that cross E F .Instead, the self-consistent nature of these GCDFT calculations allows the density of states to renormalize such that states with primarily C character become unoccupied.This behavior is most visible by comparing the top and bottom panels in Figure 2. Specifically, a state with Fe character remains immediately below E F in both panels, while the C states initially below E F without a potential completely shift above E F in the bottom panel.Similar non-rigid shifts in the Fe d orbitals can be observed, which are described further in the Supporting Information and Figure S2.
3.2.Adsorption in the Grand Canonical and Canonical Ensembles.Within the GCE at an electrochemical potential μ, the total energy of each system, Φ[μ], can be written as eq 1: −37 Consequently, grand canonical adsorption energies, Φ ads , are calculated with eq 2: where Φ A* [μ] is the grand canonical total energy of the surface with the adsorbed species of interest, Φ Surf [μ] is the grand canonical total energy of the surface (e.g., the FeN 4 @G electrocatalyst), and Φ MolcRef [μ] is the grand canonical total energy of the adsorbate's molecular reference (see Section 2 and Table 1 for the referencing scheme).For CE calculations that did not depend on the applied potential during energy minimization, all adsorption energies, E ads , are calculated using an equation in the same form as eq 2, but instead using Helmholtz total energies rather than grand canonical total energies.A post-processing correction for the standard hydrogen electrode potential is then applied if a proton/ electron pair is involved in the adsorbate's molecular reference via a H 2 molecule. 28See Table S5 for details.

The Journal of Physical Chemistry C
The qualitative results of the screening described in the previous section are shown via adsorption maps in Figure S3, and the quantitative results are summarized in Figure 3. Figure S3 shows which species adsorb to FeN 4 @G or desorb after relaxation.Apart from a GCE potential of 1 V vs SHE allowing SO 4 adsorption to C, the qualitative adsorption/desorption behavior of each ligand does not change with the potential condition.Species that adsorbed to FeN 4 @G without a potential also adsorbed to the same site(s) under GCE 0 and 1 V vs SHE, but more extreme potentials could change this behavior.Most importantly, all species initialized on Fe remain adsorbed on Fe.The adsorption of a ligand displaced Fe toward the adsorbate such that it was out-of-plane with the rest of the N-C sheet by 0.1−0.5 Å (Figure S4).Overall, the adsorption energy for NO adsorption to Fe is the lowest as it is over 1 eV more favorable than O 2 for all potential conditions.No species remained adsorbed atop N except for the H atom under a potential.Instead, the adsorbates migrated in order to directly bind to Fe.
In contrast to N adsorption, there was wide qualitative variation in whether species could remain adsorbed to C sites (Figures S3 and S5).For instance, NO and H 2 O selectively adsorbed to only Fe, while H, O, and OH bound to nearly every C atop site included in this screening.The bonding model of selective NO adsorption to Fe, for example, can be briefly summarized by only the Fe atom d orbitals having the correct symmetry to allow electron back-bonding (see Figure S6 and the nearby text for details).Because of this and the low predicted NO adsorption energies, experiments using NO to probe for the presence of FeN 4 sites should expect selective uptake to Fe if one assumes a similar structure to the model studied here.Although some species such as O and OH can remain adsorbed to C, C adsorption is at least 1 eV less favorable than Fe adsorption for the same adsorbate.Consequently, favorable OH adsorption to C only occurs above GCE 1 V vs SHE, while favorable O adsorption to C only occurs below 0 V vs SHE.This suggests that graphene oxidation would be driven by either O or OH adsorption, but not both, at potentials outside a window of approximately 0−1 V vs SHE.The O and OH C adsorption energies are affected by the local C bonding environment (Figure 3), which is discussed further for OH in the Supporting Information (Figure S5 and the nearby text).The OH adsorption energy is the most strongly affected by the choice of molecular reference.
For adsorption to Fe, we find that the adsorption energy varies significantly with the adsorbing species and potential condition/ensemble choice.The molecular references used here for O, O 2 , H 2 O, and NO (Table 1) do not include a proton/electron pair, meaning that their CE adsorption energies are not potential-dependent.For the rest of the adsorbates, their CE molecular references do include the energy of a proton/electron pair via half the energy of a H 2 molecule, leading to changes in their adsorption energies when the potential is not 0 V vs SHE due to the post-processing correction.These molecules' GCE adsorption energies change with potential due to surface charging and are also affected by a change in molecular references that include charged cations and anions.We note that the choice of molecular reference species for charged adsorbates is not unique, and other molecular referencing schemes exist, with two such examples shown in Figure S7. 38The accuracy of our implicit solvation model in describing charged molecular reference species and the specific impact of surface charging on adsorption energy is quantified further in the next section and Section 2.
The adsorption energy for H, O, O 2 , and NO is all less favorable under GCE 1 V vs SHE than under GCE 0 V vs SHE, while the adsorption energy for OH, H 2 O, SO 4 , HSO 4 , and ClO 4 is all more favorable.The adsorption energy for H was positive for both applied potentials.These trends impact which ligands can compete with O 2 for available Fe sites.Monodentate O 2 adsorption to Fe is 0.15 to 0.28 eV, which is more favorable than bidentate adsorption; however, depending on the potential, monodentate O 2 adsorption can be less favorable than OH, NO, SO 4 , and O adsorption.The H 2 O and HSO 4 adsorption energies under GCE 0 V vs SHE (−0.56 and −0.30 eV, respectively) are less favorable than O 2 adsorption (−0.99 eV) but are more favorable than O 2 (−0.62 eV) under GCE 1 V vs SHE (−0.79 and −0.78 eV, respectively).The stronger H 2 O adsorption to Fe as compared to O 2 under GCE 1 V vs SHE is particularly notable because it implies that H 2 O formed from ORR OH hydrogenation may be persistent and thus block O 2 from subsequent Fe adsorption and reduction.Furthermore, SO 4 binds even more strongly to Fe than HSO 4 at GCE 1 V vs SHE (−1.10 eV).It may not occupy as many sites, however, because the HSO 4 pK a is 1.96, meaning that HSO 4 will be approximately 92 times more abundant in a pH = 0 solution than SO 4 .ClO 4 selectively adsorbs to Fe with a similar structure to SO 4 but has an adsorption energy of 0.02 eV higher than O 2 under a GCE 1 V vs SHE potential, suggesting that it will not poison Fe sites to the same extent as SO 4 and HSO 4 .
Overall, the similar adsorption energies of OH, SO 4 , O, H 2 O, SO 4 , and HSO 4 to the O 2 adsorption energy imply that this active site may not be highly active or selective toward the ORR, in contrast with the high experimentally observed activities for this type of electrocatalyst.However, this discrepancy could be explained by multiple factors that all warrant future detailed investigation.Axial spectator ligands, explicit water molecules/solvation model choice, other active site models, and DFT functional choice could all significantly alter the predicted adsorption energies of these ligands.Explicit water molecules, for example, have been shown to interact with adsorbed intermediates, bind as an axial spectator to the Fe, and affect the predicted ORR mechanism in recent canonical ensemble calculations. 39The results shown here provide insight for this catalyst model and set of computational choices.
One structure modification that emerges from the inclusion of grand canonical applied potential with this set of molecular references is the presence of one or two lateral spectator OH groups adsorbed on a C atom close to Fe. Spectator OH adsorption on C is favorable at slightly higher GCE applied potentials (−0.17 eV under 1.23 V vs SHE), and we hypothesize that defects would allow OH adsorption to C at lower potentials.Here, we calculate how the adsorption energies of the above adsorbates change on two new active sites: the first with one lateral spectator OH placed on the lowest energy C adsorption site (Figure S8A) and another with a second spectator OH placed on the 180°rotationally and symmetrically equivalent C atom on the same side of the graphene layer (Figure S8C).To estimate whether the shifts in adsorption energy originate from a direct steric effect or steric effects and other effects such as charge transfer, the OH groups were placed either on the same side of the graphene layer as the adsorbate or on the opposite side of the graphene layer.

The Journal of Physical Chemistry C
One lateral spectator OH slightly increases the spin moment of Fe under GCE 0 V vs SHE, while two lateral spectator OH slightly increase the spin moment under GCE 0 and 1 V vs SHE.
Summaries of the changes in adsorption energies for O, OH, monodentate O 2 , SO 4 , and HSO 4 as a function of potential are shown in Figure 4.The CE shift will be the same for either the CE 0 or 1 V vs SHE correction.The effects of spectator OH on other ligands are shown in Figure S8B,D.We find that 2 spectator OH ligands typically weaken adsorption to Fe by 0.1 to 0.6 eV, while one spectator OH ligand has a smaller effect.The effect of the GCE potential on adsorption energy shift is the opposite as the trend in clean site adsorption energy with potential for all adsorbates except O 2 .For example, SO 4 adsorption on clean FeN 4 @G becomes more favorable with increased GCE potential, while OH spectators destabilize SO 4 more with increased GCE potential.For systems with one OH spectator, the side of spectator adsorption typically has a small effect.The NO molecule adsorption energy changes most and is likely related to the Fe−N−O angle changing from nearly 180°with no spectators to ∼150°with one OH spectator.For systems with two OH spectators, spectators on the same side as the adsorbate are most destabilizing, with SO 4 , HSO 4 , and ClO 4 being among the most affected due to their larger size.Interestingly, with two same-side spectators present, O 2 is destabilized by a similar amount (∼0.2 eV) as SO 4 , HSO 4 , and ClO 4 , while the adsorption energy of H 2 O is hardly affected, implying that future work focusing on H 2 O desorption may be needed to determine this type of active site's ORR activity.

Origin of Adsorption Energy
Changes due to an Applied Potential.The shifts in adsorption energies due to the inclusion of the GCE applied potential as compared to applying a potential within the CE, ΔΦ ads [μ], are nonuniform.For example, the adsorption energy of OH on Fe is stabilized by 1.51 eV when a GCE 0 V vs SHE potential is applied as compared to the CE 0 V vs SHE case, while the adsorption energy of O 2 on Fe is destabilized by 0.17 eV.The GCE adsorption energies shift relative to the CE adsorption energies because of a combination of system charging and concurrent structural relaxation.To quantify the energetic change due to surface charging, we calculated ΔΦ X [μ] for each term in the adsorption energy (X = surface (Surf), adsorbate references (MolcRef), and surface + adsorbate (A*)), comparing the GCE total energy of each system, Φ X [μ], to the total energy of the analogous CE system where surface charging was forbidden (see Table S5 for an example calculation).Because grand canonical total energies calculated using eq 1 include a μN term that accounts for the energy of all system electrons at a specified potential, the canonical DFT total energy of each system without any charging, A DFT, CE , must be shifted by a similar term in order to compare them.The number of electrons used for this shift will be an integer by definition.As a result, ΔΦ[μ] was calculated using eq 3: where N CE is the number of electrons in the calculation for A DFT, CE .We note that the screened ligands adsorbed in the same qualitative manner to the Fe site in both the CE and GCE and thus facilitated more relevant energy comparisons.
The results of calculating eq 3 for each term in the CE and GCE adsorption energies are summarized in Figures 5 and S9.The calculation of ΔΦ ads [μ] in this work involves the grand canonical total energy of the clean FeN 4 @G surface.Figure 5 shows that surface charging stabilizes clean FeN 4 @G for both 0 and 1 V vs SHE and that the stabilization is larger for the more positive 1 V vs SHE case further away from the FeN 4 @G PZC (blue data).Recall that FeN 4 @G is more oxidized under   3was used to calculate these energy differences to compare the GCE 0 and 1 V vs SHE calculations to the CE 0 and 1 V vs SHE calculations.Negative values indicate that surface charging (for the clean site and adsorption system energies) or the change in molecular references lowered the energy for that term in the adsorption energy equation.Note that although surface charging may lower the energy of all three components of the adsorption energy, the total ΔΦ ads [μ] may still be positive because the reference states could be stabilized more in total than the adsorbed system.The effect of surface charging on the molecular references, the clean FeN 4 @G surface, and the adsorbed systems are shown with the magenta, blue, and orange lines, respectively.The total effect of all three of these contributions is shown with the black lines using an equation in the same form as eq 2 for the adsorption energy.The lines drawn between the 0 and 1 V vs SHE data are only guides to the eye and do not represent the approximately quadratic relationship between ΔΦ[μ] and potential for surface systems.The 0 and 1 V x axis groups denote the difference between either the 0 or 1 V vs SHE grand canonical ensemble data and the canonical ensemble data.

The Journal of Physical Chemistry C
a more positive applied potential (Table 2).In general, ΔΦ[μ] is at least quadratically related to both the system's charge state due to surface charging, (N CE − N), and the difference between the system's PZC and the applied potential (U − U PZC ). 14,40,41Indeed, we observe an approximately quadratic relationship between ΔΦ[μ] and N CE − N across all of the systems with molecules adsorbed to Fe (Figure S10).As a result, energy minimization within the GCE stabilizes the FeN 4 @G energy relative to the CE to a greater extent for the more oxidative potential.
Next, Figure 5 shows that the inclusion of the applied potential for the adsorbed systems (gold data) tends to stabilize them by a similar amount as for the clean FeN 4 @G surface, which is expected given the low coverage of each ligand on FeN 4 @G and lack of large structural changes to FeN 4 @G.As a result, similar stabilizations of the clean and adsorbed systems will cancel when calculating an adsorption energy.However, this is not always the case as the adsorbed states of O, SO 4 , HSO 4 , and ClO 4 are stabilized to a much different extent than the clean electrocatalyst, particularly under 1 V vs SHE.The FeN 4 @G + HSO 4 -adsorbed system, for example, is stabilized by only 1.05 eV, as compared to the 1.79 eV stabilization of the clean surface, resulting in a net adsorption energy destabilization of 0.74 eV due to differences in surface charging.We note that the PZCs of some adsorbed systems are coincidentally close to the applied potentials.For example, the PZC for HSO 4 adsorbed on Fe (FeN 4 @G + HSO 4 ) is only 0.12 V vs SHE (Figure S11), meaning that energy minimization within the GCE results in very little surface charging, and thus energy change, relative to the CE (Figure 5, left yellow triangle for HSO 4 is only −0.02 eV).
Finally, because the molecular reference states for O, O 2 , H 2 O, and NO are the same (neutral) species (Table 1) for both the CE and GCE calculations, the μN and μN CE terms in eq 3 exactly cancel, meaning that surface charging completely explains differences in the CE and GCE adsorption energies.Adsorption energy calculations for these molecules are unambiguous because they each have a neutral, molecular reference state that can be present in the solution phase.In contrast, charged molecular references such as H + can be referenced to either the energy of half a H 2 molecule, as is commonly done within the computational hydrogen electrode (CHE) approach within the CE or by subtracting the total energy of H 2 O from H 3 O + within the GCE. 28In this work, we used charged molecular reference species for H, OH, SO 4 , HSO 4 , and ClO 4 (Table 1) for the grand canonical adsorption energies because these molecules could all exist in a bulk electrolyte in a real electrochemical environment.This choice is complemented by the CANDLE implicit solvation model used in these calculations.Its detailed treatment of the asymmetry in the solvation of cations and anions manifests as significant improvement in the description of solvated anions over other common solvation methods. 21This can be seen by its accurate prediction of the pK w of H 2 O and pK a values of H 2 SO 4 , HSO 4 − , and HClO 4 (see Section 2 and Table S1 for further details).We used only neutral molecules for the CE adsorption energies, with the energy of a proton/electron pair defined by the energy of the H 2 molecule, as is standard practice within the CHE.These two sets of molecular reference states have differing numbers of electrons associated with the molecular reference energies of H, OH, SO 4 , HSO 4 , and ClO 4 .However, the post-processing potential correction commonly applied to the CE molecular reference states for these molecules still captures the same potential dependence as for their charged molecular reference states, meaning that the result of eq 3 for the molecular reference states (magenta data) is not potential-dependent.The magenta data is thus not 0 when the total energies of an adsorbate's CE and GCE molecular references in Table 1 are not equal.For OH, this discrepancy is large, while for HSO 4 , the discrepancy is very small.
When the change in molecular reference energy is combined with the change in energy due to surface charging (black data), some adsorbed states, such as FeN 4 @G + NO and FeN 4 @G + HSO 4 , do not have a large net change in adsorption energy, while others, such as FeN 4 @G + H and FeN 4 @G + OH, do.For systems like FeN 4 @G + NO, the fortuitous cancellation of possibly large surface charging effects sometimes allows CE adsorption energies to be close in energy to their grand canonical counterparts.In general, however, this may not necessarily be the case, and the grand canonical formalism is expected to be needed in order to accurately model diverse electrochemical systems.

CONCLUSIONS
In this work, we used GCDFT to predict the thermodynamics of ligand adsorption to an FeN 4 @G electrocatalyst.This approach accounts for the surface charging that occurs in a real electrochemical cell by allowing the number of electrons in each of the screened systems equilibrate self-consistently at a specific applied potential.We found that the use of the GCE significantly alters the predicted electronic structures, charge states, and adsorption energies of the screened systems as compared to adsorption energies calculated within the canonical ensemble.The adsorption energies of several ligands under an applied potential were lower than O 2 , suggesting that other ligands may be persistently bound to the Fe during ORR.Additionally, these adsorption energies were modulated by the presence of spectator OH groups near the Fe.Comparison of GCE total energies with CE total energies reveals that changes in adsorption energy occur because surface charging alters the energetics of the reference states (clean active site and molecular reference) and adsorbed systems unequally.Notably, we find that the choice of molecular reference states is important and compared the effect of using charged or uncharged species.The accurate prediction of experimental properties using charged references with the CANDLE solvation model suggests that this approach may be a viable alternative to the use of neutrally charged species.However, as described above, the impact of the referencing scheme and advanced solvation models on electrocatalyst activity should be further studied using more advanced hybrid or "beyond-DFT" functionals that better predict physically realistic energy gaps while incorporating long-range bonding interactions. 42Future work incorporating GCE reaction barriers and microkinetic models of the ORR elementary steps, both with and without spectator ligands, is necessary to fully predict the activity of this electrocatalyst model.

The Journal of Physical Chemistry C
Analysis of the referencing scheme chosen, adsorbate bonding, and adsorbed system structures and electronic properties (PDF) Lattice and atomic coordinate data in the POSCAR format for the lowest energy structures (adsorbates bound to Fe) for the three potential conditions (ZIP) ■ AUTHOR INFORMATION

Figure 1 .
Figure 1.Adsorption sites screened in this work for the FeN 4 @G electrocatalyst.The gray box denotes the boundaries of the modeled supercell for the FeN 4 @G active site model.Gray, blue, and gold circles represent C, N, and Fe atoms, respectively.Yellow triangles denote atop sites, purple squares denote bridge sites, and green hexagons denote hollow sites.

Figure 2 .
Figure 2. PDOS of the FeN 4 @G electrocatalyst using the canonical or grand canonical ensembles.PDOS of each element type with no potential (canonical ensemble, top panel), 0 V vs SHE (grand canonical ensemble, middle panel), or 1 V vs SHE (grand canonical ensemble, bottom panel).The gray data represents C states, the blue data represents N states, and the gold data represents Fe states.

Figure 3 .
Figure 3. E ads or Φ ads for all adsorbed ligands on the FeN 4 @G electrocatalyst for all potential conditions.Only symmetrically unique relaxed structures are shown.Points without a horizontal split indicate monodentate adsorption, while points with a horizontal split indicate bidentate adsorption.O 2 and SO 4 bidentate adsorption occurs through two different O atoms to different catalyst sites but was always less favorable than monodentate adsorption.Gray points denote ligand adsorption on C, blue points denote ligand adsorption on N, and gold points denote ligand adsorption on Fe.The adsorption energies are grouped along the x axis by the potential condition they were run at.The adsorption energies are grouped along the x axis by if they were run within the canonical ensemble with a post-processing correction for potential (CE 0 V and CE 1 V), grand canonical ensemble 0 V vs SHE (GCE 0 V), or grand canonical ensemble 1 V vs SHE (GCE 1 V).

Figure 4 .
Figure 4. Effect of 2 lateral spectator OH ligands on adsorption to Fe in the FeN4 @G electrocatalyst under all potential conditions.Open circles represent 2 spectator OH spectators on the opposite side of the graphene and filled circles denote 2 spectator OH on the same side of the graphene.The adsorption energies are grouped along the x axis by if they were run within the canonical ensemble (CE; blue), grand canonical ensemble 0 V vs SHE (0 V; gold), or grand canonical ensemble 1 V vs SHE (1 V; red).These adsorption energy shifts for the canonical ensemble are independent of the post-processing correction for potential.The results for additional ligands and 1 OH spectator are shown in FigureS8.

Figure 5 .
Figure 5. Decomposition of contributions to ΔΦ ads [μ] by the adsorbed system, clean surface, and molecular references.Equation3was used to calculate these energy differences to compare the GCE 0 and 1 V vs SHE calculations to the CE 0 and 1 V vs SHE calculations.Negative values indicate that surface charging (for the clean site and adsorption system energies) or the change in molecular references lowered the energy for that term in the adsorption energy equation.Note that although surface charging may lower the energy of all three components of the adsorption energy, the total ΔΦ ads [μ] may still be positive because the reference states could be stabilized more in total than the adsorbed system.The effect of surface charging on the molecular references, the clean FeN 4 @G surface, and the adsorbed systems are shown with the magenta, blue, and orange lines, respectively.The total effect of all three of these contributions is shown with the black lines using an equation in the same form as eq 2 for the adsorption energy.The lines drawn between the 0 and 1 V vs SHE data are only guides to the eye and do not represent the approximately quadratic relationship between ΔΦ[μ] and potential for surface systems.The 0 and 1 V x axis groups denote the difference between either the 0 or 1 V vs SHE grand canonical ensemble data and the canonical ensemble data.
3 O + , OH − , SO 4 OH, and H 2 O could use an alternative referencing scheme based on the overall reduction reaction of O 2 .The energy of a proton for the GCE adsorption energies was found by subtracting the total energy of H 2 O from H 3 O + and does not depend on the GCE applied potential.This is because H 2 O and H 3 O + have the same number of electrons and thus identical potential-dependent contributions to their total energies.
−were all used as molecular references because these molecules could all exist in a bulk electrolyte in a real electrochemical environment, and their associated pK w and pK a values were predicted accurately.Although not the focus of this work, calculations comparing ORR thermodynamics between O 2 , O,

Table 1 .
Molecular Reference States for Each Ligand in the Screening a

Table 2 .
Summary Charge States and Spin Moments for the Atom Types in the FeN 4 @G Electrocatalyst under Different Potential Conditions a