Comparing B3LYP and B97 Dispersion-corrected Functionals for Studying Adsorption and Vibrational Spectra in Nitrogen Reduction.

Electrochemical ammonia synthesis is being actively studied as a low temperature, low pressure alternative to the Haber-Bosch process. This work studied pure iridium as the catalyst for ammonia synthesis, following promising experimental results of Pt-Ir alloys. The characteristics studied include bond energies, bond lengths, spin densities, and free and adsorbed vibrational frequencies for the molecules N2 , N, NH, NH2 , and NH3 . Overall, these descriptive characteristics explore the use of dispersion-corrected density functional theory methods that can model N2 adsorption - the key reactant for electrochemical ammonia synthesis via transition metal catalysis. Specifically, three methods were tested: hybrid B3LYP, a dispersion-corrected form B3LYP-D3, and semi-empirical B97-D3. The latter semi-empirical method was explored to increase the accuracy obtained in vibrational analysis as well as reduce computational time. Two lattice surfaces, (111) and (100), were compared. The adsorption energies are stronger on (100) and follow the trend EB3LYP >EB3LYP-D3 >EB97-D3 on both surfaces.

Introduction

In 1954, about 2,900,000 metric tons of ammonia were produced in the United States, 70 % of which was utilized by the agricultural sector. More than 60 years later, in 2018, U.S. ammonia production reached an estimated 12,500,000 metric tons, with the world total coming in at 140,000,000 tons. This essential industrial process currently uses over 1 % of the world's power. Ammonia is used in the synthesis of a variety of products, most notably, fertilizers. The Haber‐Bosch process, the process by which ammonia is produced, is a large source of carbon emissions primarily because the synthesis reaction is performed at very high temperatures and pressures (350–550 °C and 150–350 atm). Researchers are looking for alternative mechanistic pathways to bypass the energy requirements of Haber‐Bosch. Such processes include biological nitrogen fixation, ultraviolet promotion of N2, and electrochemical nitrogen reduction.[ , , , , ] For the industrial Haber‐Bosch process, iron and ruthenium are commonly studied[ , , ] because transition metals readily adsorb and release oxygen. Catalysts studied for the electrochemical process include ruthenium, iron oxides, and nitride compounds.[ , , , ] Osmium and platinum have also been identified as good catalysts.[ , ] Computational techniques have made it easy to screen a large number of metal catalysts and their anchoring supports for desired properties. Recent supports successfully used to study electrochemical nitrogen reduction include graphene and phosphorene.[ , ]

Experimentally, Allagui et al. found that a bimetallic PtIr nanoparticle catalyst decomposed ammonia at a 33 % higher rate than platinum alone. Le Vot et al. found that iridium lowers the overpotential of ammonia oxidation. Boggs and Botte came to the same conclusion in their study of PtIr anodes deposited on carbon substrates in alkaline solution. Estejab and Botte used small cluster calculations to show Ir3 is more active than Pt3 in nitrogen oxidation. These studies showed success with Pt‐Ir for nitrogen oxidation, while Sheets and Botte found success with Pt‐Ir for nitrogen reduction as well. These studies serve as a basis for exploring mechanistic behavior on larger clusters of pure iridium for nitrogen reduction. By better understanding the thermodynamics of the surface, new steps can be taken in the lab to increase ammonia yield.

When using clusters to represent catalysts, the computational scheme is very important. The popular B3LYP Density Functional Theory (DFT) method does not describe dispersion between molecules well and thus cannot appropriately describe N2 adsorption on a catalyst surface. Hopmann et al. computationally studied bond formation and ligand exchange reactions mediated by iridium in solution. They found that the free energies of B3LYP systems with dispersion provide moderate overall accuracy and significant improvement over B3LYP alone. Therefore, this study tests B3LYP against its dispersion‐corrected form, as well as against a less computationally expensive dispersion‐corrected functional, B97‐D3.

This paper has two primary objectives. The first is to calculate adsorption energies, binding sites, and vibrational spectra for the intermediates involved in the ammonia synthesis reaction on pure iridium. This data can be used to shed light on the reaction mechanism of ammonia synthesis on the (111) and (100) surfaces. The activation energies, which correlate to the speeds of the possible reactions, can also be calculated, but they are beyond the scope of this study. Calculations on pure iridium will serve as a direct comparison to previous calculations performed on pure platinum clusters. Key differences will be analyzed with respect to adsorption energetics and bonding sites. These calculations can also serve as a benchmark for future bimetallic calculations. Calculations on iridium are done in the absence of an electrochemical environment, similar to prior calculations on platinum. The aim of the study is to directly compare the two studies and gather insights on the relative differences between bonding characteristics on Pt and Ir. However, vibrational frequencies and adsorption energies will shift in an electrochemical environment, so the effect of voltage and solvent molecules on adsorption is left to future studies. Furthermore, the insight gained from these heterogeneous catalytic studies are a necessary foundation for future electrochemical studies.

Significant energy savings would result from an electrochemical process that takes place at ambient temperatures and pressures.[ , ] Thus, all calculations in this study are run at standard temperature and pressure (STP). The second objective is to compare the accuracy of three Density Functional Theory methods – B3LYP, B3LYP‐D3 and B97‐D3 – in predicting bonding sites, bonding lengths, bonding energies and vibrational frequencies on the catalyst surface. These characteristics are well‐modeled by cluster calculations.[ , , ]

Computational Methods

All calculations were implemented in Gaussian 09 using default convergence limits. The computational schemes compared in this study are B3LYP, B3LYP‐D3, and B97‐D3. B3LYP is the Becke‐style 3‐Parameter functional.[ , ] The dispersion corrected version, B3LYP‐D3, uses Grimme's dispersion and the Beck‐Johnson damping parameter. B97‐D3 is Grimme's functional with Becke‐Johnson damped dispersion as well.[ , , ] While B3LYP is a hybrid functional, often used to achieve bond lengths and vibrational frequencies close to experimental values, B97 is a semi‐empirical pure functional which delivers speedier results on larger clusters of atoms.

For consistency between methods and with previous calculations on platinum, we will use the LANL2DZ basis set for the iridium atoms and 6–311++g** basis sets for nitrogen and hydrogen.

The (111) and (100) surfaces of the face‐centered cubic iridium catalyst are both modeled in this study. As a default in Gaussian 09, the most abundant isotope (Ir‐193) is used for calculations. The clusters built for these calculations are kept at frozen coordinates during optimization, with an atomic Ir‐Ir spacing of 0.3839 nm. All adsorbates are unfrozen. Frequency calculations are conducted at the minimum energy geometries obtained from optimization. Bending modes with an IR intensity less than 20×10−40 esu2 cm2 or Raman activity less than 20 A4/AMU are considered difficult to observe against the cluster vibrations and are therefore not all recorded in this paper. For general comparisons, this paper classifies vibrations as weak (activity≤40), moderate (activity=40–100), strong (activity=100–200), and very strong (activity=200+).

Cluster Geometry

A cluster size of 15 Ir atoms was chosen by analyzing the cohesive energy and spin density of Irn clusters (n=10, 15, 20, 25, 30) shown in Figure 1.

Figure 1

Irn clusters modelled, where n=10,15, 20, 25, and 30. The left column depicts the (111) surfaces and the right column depicts the (100) surfaces. Atoms are numbered so they can be matched with their respective spin densities in Figure

Each Ir atom has nine valence electrons spread amongst its 5d and 6 s orbitals. An initial estimate for the number of unpaired electrons in our modeled catalyst was determined using Kua and Goddard's Interstitial Electron Model. Their model estimates the unpaired electrons for the (111) surface and hypothesizes that the surface has s1d8 configuration. This methodology was used to describe the (100) surface as well, as an initial guess. The number of unpaired electrons in the ground state and the energies of each cluster are calculated and presented in Table 1.

Table 1

Cohesive energy of iridium clusters with n atoms, where n=10, 15, 20, 25, and 30. The ground state spin multiplicity was predicted using the Interstitial Electron model and verified by testing two multiplicities above and below the ground state. The lowest energy state is presented here for each cluster size.

Ground state spin multiplicity of (111) Irn clusters.				Ground state spin multiplicity of (100) Irn clusters.
Method	Irn	Unpaired e‐	Ecohesive (kJ/mol)	Method	Irn	Unpaired e‐	Ecohesive (kJ/mol)
B3LYP	10	18	−292	B3LYP	10	18	−286
	15	27	−319		15	27	−296
	20	22	−338		20	14	−322
	25	F	F		25	F	F
B3LYP‐D3	10	20	−318	B3LYP−D3	10	18	−312
	15	27	−352		15	27	−328
	20	22	−377		20	14	−362
	25	F	F		25	F	F
B97‐D3	10	16	−385	B97‐D3	10	16	−378
	15	19	−391		15	25	−393
	20	20	−447		20	14	−436
	25	19	−467		25	27	−449
	30	F	F		30	22	−467

F ‐ Calculation failed to converge.

Cluster energies were calculated for two multiplicities above and below the ground state multiplicity to evaluate convergence to a minimum. Clusters larger than 20 atoms did not converge using B3LYP or B3LYP‐D3. Only the (100) surface of B97‐D3 converged for clusters larger than 25. As shown in Table 1, the larger the cluster, the closer the cohesive energy comes to approaching the experimental bulk cohesive energy of iridium, 670 kJ/mol.[ , ] (111) surfaces are more tightly packed and thus approach the bulk cohesive energy quicker than (100) surfaces. A cluster size of 15 atoms is ultimately chosen over a size of 20 atoms because the spin densities across the 15‐atom cluster are closer to the number of unpaired electrons per atom in the ground state d‐orbital i. e. two (see Figure 2).

Figure 2

Spin density of Irn clusters (n=10, 15, 20, 25). Columns 1, 2, and 3 show the cluster spin densities using B3LYP, B3LYP‐D3, and B97‐D3 respectively. The 30‐atom cluster only converged for the (100) surface at the B97‐D3 level of theory, so it was not included in this comparison.

The 15‐atom cluster is chosen for both the (111) and (100) surfaces; they are pictured in Figure 1c and 1d respectively. They are both two‐layer clusters with enough surface area that all principal adsorption sites on the surface are distanced from the catalyst's edge, reducing “edge effects.”

Simulation Procedure

First, the molecules N2, N, NH, NH2, and NH3 are optimized in the gas‐phase using each of the three DFT functionals. The optimized molecules are then placed on the (100) and (111) catalyst clusters at one of the preferential sites in Figure 3.

Figure 3

Preferential sites for adsorption on the (111) (above) and (100) (below) Ir15 surfaces. B=Bridge, T=Top, HCP=Hexagonal Close‐Packed, FCC=Face‐Centered Cubic, and H=Hollow.

Each molecule is optimized separately at every possible preferential site in order to determine the lowest energy position on the catalyst. For N2 optimizations, that includes testing end‐on and side‐on N−N configurations at each site. After optimization, frequency calculations are run on the lowest energy adsorption sites. This data is then compared between functionals and to other theoretical and experimental benchmarks.

Results and Discussion

Binding energy is calculated from equation 1:

where binding energies are based on electronic energies for calculated preferential adsorption sites, and on zero‐point and thermal‐corrected energies elsewhere.

Preferential Adsorption Sites

The molecules successfully adsorb to the surface of each catalyst at the sites show in Figure 4. About one hundred structures were successfully optimized, with nine failures that were due to an inability to converge or bond to the catalyst surface. When two different starting geometries converged to the same final site, the lower energy optimization was used in calculations for Figure 4.

Figure 4

Binding energy (based on electronic energies) of molecules adsorbed on the (111)[a–d] and (100)[e–h] surfaces using the three different functionals. B=Bridge, T=Top, HC=Hexagonal Close‐Packed, FCC=Face‐Centered Cubic, and H=Hollow.

Ir15 (111) Cluster

N and NH adsorb at top, hexagonal close‐packed, and face‐centered cubic positions using all methods – B3LYP, B3LYP‐D3, and B97‐D3. The HCP position is most favored, closely followed by FCC. This is very similar to the theoretical findings of Krekelberg et al., whose study used periodic DFT methods PW91 and RPBE. They found that HCP and FCC were equally favorable. NH2 adsorbs at the bridge position for all methods and the top position for B3LYP‐D3 and B97‐D3, with the bridge position being slightly more favorable. NH3 also adsorbs at top and bridge positions for all methods, and the top position is most favorable. This agrees with Krekelberg's theoretical findings, but is in conflict with an experiment that found adsorption of NH3 at FCC/HCP. Although for most of the molecules in this study, more than one adsorption site is found possible, none of the methods could find a minimum at the HCP or FCC position for NH3.

Ir15 (100) Cluster

N adsorbs at the top and bridge positions for all methods and for B97‐D3, the hollow position as well; however, the bridge position is the most favorable across all methods. For NH, B3LYP‐D3 does not converge to a bridge site, but B3LYP and B97‐D3 do. B97‐D3 shows an additional, less favorable adsorption at the top position. Both B3LYP and B3LYP‐D3 converged to hollow sites located only at the edges of the surface. Due to edge effects on this hollow site, it is likely that bridge is the most favorable site for NH. This is in agreement with adsorption of monatomic O on Ir(100), which also has two unpaired electrons and adsorbs at the bridge position. NH2 and NH3 have only one favorable site – bridge and top respectively.

Overall, all methods yield the same preferential adsorption sites on both the (100) and (111) surfaces, with the exception of NH on (100). B3LYP always gives the smallest binding energies, while its dispersion‐corrected form gives larger binding energies and B97‐D3 gives the largest binding energies.

N2

N2 comparisons are limited to B3LYP‐D3 and B97‐D3. B3LYP simulations do not have the dispersion calculations necessary to adsorb N−N to the catalyst, as DFT cannot accurately describe the types of dispersion forces that would exist between N2 and the metal surfaces. Therefore, no N2 simulations converged with B3LYP. The binding energies calculated for the other methods are shown in Figure 5 with respect to the angle of adsorption on the catalyst surface (N−N bond with Ir surface).

Figure 5

Binding energy (based on electronic energies) of N2 adsorbed on the (111) and (100) surfaces at various angles. Θ is the angle the N−N bond makes with the catalyst surface.

B97‐D3 suggests that the end‐on configuration is favorable for both (111) and (100) surfaces. This position hovers at the top site, in agreement with Krekelberg. B3LYP‐D3 suggests that there are multiple orientations of N2 within 0.05 eV that successfully adsorb, with the lowest energy configurations being mostly parallel to the surface (side‐on) for both (111) and (100).

Frequency Calculations and Bonding Characteristics

The most favored adsorption site for each of the molecules in Figure 4 – also illustrated in Figure 6 – is used for a frequency calculation.

Adsorption of N2

Diatomic nitrogen adsorbed on the two surfaces is pictured in Figure 6e and 6j. This end‐on, top position corresponds to the lowest energy optimization achieved with B97‐D3. The N−N bond makes an angle of 89° with the (111) and (100) surfaces.

Figure 6

Top view of preferred adsorption sites for adsorbed N, NH, NH2, and NH3 on the (111)[a‐d] and (100)[f–i] surface. A side view is given for N2 [e and j]. Images shown are from the B97‐D3 optimizations, but site preferences are consistent across all methods. However, the orientation of N2 is different using B3LYP‐D3.

The B3LYP−D3 N−N bonds make angles of 30° and 7° with the (111) and (100) surfaces respectively. The spin density per atom before and after adsorption of N2 can be found in Supplementary Data Figure S1. With B3LYP‐D3, the spin density profile is consistent before and after adsorption, with a slight reduction in spin on atom 5 for the (111) surface because the N−N bond is angled 30° towards that atom, providing context for slight charge transfer between the N and Ir atom. With B97‐D3, the (111) surface sees a similar effect. Notably, there is a greater variance in B97‐D3 ground state spin density before and after bonding than with B3LYP‐D3. On the (100) surface, the spin density drop is particularly steep (N2 oriented 89°). This suggests bonding of the inner N atom to Ir atom number 5. This bond results from an interaction between the Ir 5d‐states and π‐bonds polarized in the direction of Ir, thus weakening the N−N bond and laying ground for N−Ir electron pairing. Adsorbate energy is still relatively small due to σ‐bond repulsion along the axis of the bond. Covalent bond formation is also evidenced by the bond lengths in Table 2, where N2 molecules that adsorbed preferentially at the top position, perpendicular to the catalyst (B97‐D3), adsorbed greater than 1 nm closer than those adsorbed at an angled or horizontal position (B3LYP‐D3). Bond energies are all far from the experimentally determined 9.8 eV, but similar to theoretical predictions of −0.54 and −0.21 eV.

Table 2

Bond energies and bond lengths are given for the most favorable orientation of N2 on the (111) and (100) surfaces.

Lattice	Method	Bond Energy (eV)	Bond Length (Å)
111	B3LYP	F	F
	B3LYP‐D3	−0.12	3.09
	B97‐D3	−0.18	1.94
100	B3LYP	F	F
	B3LYP‐D3	−0.24	3.28
	B97‐D3	−0.32	1.93

Vibrational analysis of N2 as a free molecule is broken down by method in Table 3.[ , , ] B3LYP(‐D3) predicts the stretching mode of N2 to within 10 cm−1, when vibrational scaling factors are taken into account. The scaling factor for B3LYP(‐D3) is 0.96 and for B97D3 is 0.98. These should be used to compare the results in this paper's vibration tables to any experimental observations. After adsorption, B3LYP‐D3 predicts a ∼10 cm−1 decrease in the frequency of the stretching mode, while B97‐D3 predicts an ∼170 cm−1 decrease in the frequency. Experiment gives the adsorbed frequency to be 2185/2210/2223 cm−1, which is closest to the B97‐D3 prediction of 2,140 cm−1 (scaled).[ , ] From this we can conclude B97‐D3 correctly predicts the lowest energy position on the surface is at the top position, with the N2 molecule perpendicular to the surface.

Table 3

Vibrational modes and their corresponding frequencies (in cm‐1) for free and adsorbed N2 on both the (111) and (100) surface. Values are compared to existing experimental and theoretical observation.

Vibrational mode	Gas Phase (cm−1)[†]				Adsorbed Phase (cm−1)[†]
	B3LYP	B3LYP‐D3	B97‐D3	Other Work		B3LYP	B3LYP‐D3	B97‐D3	Other Work
stretch (N‐N)	2445	2445	2365	2359 [43]		F	2434	2184	2185[44][ac*], 2210[46][ad], 2223[40][b]
stretch⊥						F	–	373	403[40][b]
(Ir‐N)
stretch (N‐N)						F	2438	2152
stretch⊥ (Ir‐N)						F	‐	404

[a] experimental; [b] theoretical; [c] IR; [d] EELS; [*] index plane not specified; [†] tabulated frequencies from this work are unscaled

2359

–

‐

Adsorption of N

Monatomic nitrogen is pictured in Figure 6a and 6 f at its preferential positions – HCP (111) and bridge (100). Average spin density before and after adsorption of N on the cluster surface can be found in Supplementary Data Figure S2. All methods show a consistent spin distribution on non‐bonding atoms before and after binding on the (111) surface. Atoms 2, 5, and 6 are the clear sites of bonding for the three unpaired electrons of nitrogen because the spin density drops by approximately one at each location. On the (100) surface, bonding is evident on only atoms 1 and 5. The third unpaired electron on N is not clearly bound to any particular Ir atom.

Bond length as described in Table 4 is the shortest distance between the N atom and the catalyst surface. Bond length for B3LYP and its dispersion‐corrected form are 1.17/1.16 Å and 1.26 Å on the (111) and (100) surfaces respectively; B97‐D3 length is slightly longer – 1.22/1.29 Å. These are close to that computed by Krekelberg – 1.16 Å and 1.11 Å. Moving from B3LYP to B3LYP‐D3 to B97‐D3, bond energy gets progressively stronger.

Table 4

Bond energies and bond lengths are given for N at its lowest energy position – HCP for the (111) surface and bridge for the (100) surface.

Lattice	Method	Bond Energy (eV)	Bond Length (Å)
111	B3LYP	−3.5	1.17
	B3LYP‐D3	−3.8	1.16
	B97‐D3	−4.0	1.22
100	B3LYP	−3.7	1.26
	B3LYP‐D3	−4.1	1.26
	B97‐D3	−4.2	1.29

Table 5 gives the IR and Raman frequencies for various vibrational modes. The asymmetric stretch on the (100) surface is very strongly Raman active, while the symmetric stretch is moderately Raman and IR active. All vibrations on the (111) surface are weak.

Table 5

Vibrational modes and their corresponding frequencies (in cm‐1) for adsorbed N at its lowest energy position – HCP for the (111) surface and bridge for the (100) surface. The displacement of N for each stretching mode is labeled as either perpendicular or parallel to the cluster surface. Values are compared to existing experimental and theoretical observation.

Vibrational mode		Adsorbed Phase (cm−1) [†]
		B3LYP	B3LYP‐D3	B97‐D3	Other work
stretch⊥ (Ir‐N)		594	595	586	546[40][a]
stretch∥ (Ir‐N)		494	498	441
stretch∥ (Ir‐N)		456	462	423
stretch⊥ (Ir‐N)		632	632	614
stretch∥ (Ir‐N)		346	361	435

[a] theoretical – RPBE; [†] tabulated frequencies from this work are unscaled

Adsorption of NH

Figure 6b and 6 g shows NH adsorbed at its most preferential site – HCP (111) and bridge (100). Supplementary Data Figure S3 shows the average spin density before and after adsorption of NH on the cluster surface. On the (111) surface, the nitrogen atom is adsorbed between atoms 2, 5, and 6, just as lone N is. The spin density per atom takes on an almost identical shape with B3LYP, B3LYP‐D3, and B97‐D3 (at a lower average spin density for the latter). For the (100) surface, B3LYP‐D3 converges to only an erroneous edge‐effect hollow position. For the other two methods, the two unpaired e− of NH form bonds with bridging Ir atoms 5 and 8. The spin structure of B97‐D3 is noticeably more similar to B3LYP after adsorption than before adsorption.

Bond length in Table 6 is consistent between the three methods for (111) – 1.23/1.24 Å. On the (100) surface, B97‐D3 is 0.07 Å shorter. As with each of the other molecules, bond energy follows the trend EB3LYP>EB3LYP‐D3>EB97‐D3.

Table 6

Bond energies and bond lengths are given for NH at its lowest energy position – HCP for the (111) surface and bridge for the (100) surface.

Lattice	Method	Bond Energy (eV)	Bond Length (Å)
111	B3LYP	−2.7	1.24
	B3LYP‐D3	−3.1	1.23
	B97‐D3	−3.3	1.23
100	B3LYP	−2.9	1.42
	B3LYP‐D3	F	F
	B97‐D3	−3.6	1.35

Table 7 shows NH frequencies for gas phase and Ir‐adsorbed states and how these frequencies compare to benchmark experimental work.[ , ] The N−H stretch perpendicular to the catalyst surface is strongly Raman active on both surfaces. On the (100) surface, the wagging mode is strongly Raman active as well, according B97‐D3 calculations. This same mode under the B3LYP method is only moderately Raman active with a stronger IR activity. All other modes are weak.

Table 7

Vibrational modes and their corresponding frequencies (in cm‐1) for free and adsorbed NH at its lowest energy position – HCP for the (111) surface and bridge for the (100) surface. The displacement of N in Ir‐N the stretches is labeled as either perpendicular or parallel to the cluster surface. Values are compared to existing experimental and theoretical observation.

Vibrational mode

Gas Phase (cm−1) [†]

Adsorbed Phase (cm−1) [†]

[a] experimental; [b] IR; F‐failed to converge; [†] tabulated frequencies from this work are unscaled

stretch (N‐H)

3133[47][ab], 3283[48][ab]

3518

3529

3455

rock (N‐H)

798

797

794

883

881

847

617

621

610

466

473

454

stretch (N‐H)

3452

3442

rock (N‐H)

880

877

635

stretch∥ (I‐N)

399

355

370

Adsorption of NH2

Figure 6c and 6 h shows NH2 adsorbed at its preferential site – bridge. Supplementary Data Figure S4 shows the average spin density before and after adsorption of NH2 to the surface. On the (111) surface, the adsorption takes place on Ir atoms 5 and 6 as evidenced by the drop in spin density on those atoms. On the (100) surface, bonding takes place on atoms 1 and 5 (B3LYP and B97‐D3), and 5 and 8 (B3LYP‐D3). Again, the spin distribution per atom for each method looks nearly identical after bonding.

Bond lengths in Table 8 are similar between methods, varying by at most 0.04 Å. EB3LYP>EB3LYP‐D3>EB97‐D3 for both surfaces.

Table 8

Bond energies and bond lengths are given for NH2 at its lowest energy position – bridge for both the (111) surface and (100) surface.

Lattice	Method	Bond Energy (eV)	Bond Length (Å)
111	B3LYP	−1.6	1.64
	B3LYP‐D3	−2.0	1.63
	B97‐D3	−2.1	1.63
100	B3LYP	−2.4	1.61
	B3LY‐D3	−2.8	1.59
	B97‐D3	−2.9	1.57

Table 9 lists the vibrational modes of gas‐phase and adsorbed NH2.[ , , ] B3LYP and its dispersion‐corrected form match experiment for the symmetric mode slightly better than B97‐D3, while B97‐D3 matches the asymmetric mode best. On both the (111) and (100) surface, the scissoring mode is IR active, while the symmetric stretch is strong and the antisymmetric stretch is weak in Raman activity. The stretching mode between N and the catalyst surface itself is very weak on both surfaces.

Table 9

Vibrational modes and their corresponding frequencies (in cm‐1) for free and adsorbed NH2 at its lowest energy position – bridge for the (111) and (100) surface. Values are compared to existing experimental and theoretical observation. Asymmetric and symmetric stretches are abbreviated as a. stretch and s. stretch.

	Vibrational mode	Gas Phase (cm−1) [†]				Adsorbed Phase (cm−1) [†]
		B3LYP	B3LYP‐D3	B97‐D3	Other Work	B3LYP	B3LYP‐D3	B97‐D3
	a. stretch (N‐H)	3439	3439	3361	3301[49,50][ab]	3558	3562	3507
	s. stretch (N‐H)	3254	3254	3168	3219,[49,50][ac,ab] 3220[47][ab]	3453	3456	3392
	scissoring (N‐H)	1507	1507	1491	1497,[49,50][ac,ab] 1499[47][ab]	1527	1526	1484
	stretch⊥ (Ir‐N)					493	504	497
	a. stretch (N‐H)					3595	3605	3513
	s. stretch (N‐H)					3497	3503	3404
	scissoring (N‐H)					1532	1525	1496
	stretch⊥ (Ir‐N)					495	506	510

[a] experimental; [b] IR; [c] Difference Frequency Laser Spectroscopy; [†] tabulated frequencies from this work are unscaled

Adsorption of NH3

Figure 6d and 6i show NH3 adsorbed at its most preferential site – top. Supplementary Data Figure S5 gives the spin density of the catalyst before and after bonding. On the (111) and (100) surfaces, ammonia adsorbs on Ir atom 5 (a top position atom). The spin density drops by ∼0.6‐1 unpaired electron at that site (the same magnitude as the drop in density for NH2), indicating NH3 may be covalently bonding with the top Ir atom. However, the adsorption energy is still much lower than for NH2.

Table 10 gives bond lengths and energies of ammonia's adsorption. Bond lengths are similar on the (111) surface. On the (100) surface, the bond length with B3LYP is 0.12 Å longer than with B3LYP−D3, whose bond length is 0.12 Å longer than with B97‐D3. The trend in binding energy is EB3LYP>EB3LYP‐D3>EB97‐D3.

Table 10

Bond energies and bond lengths are given for NH3 at its lowest energy position – top for both the (111) and (100) surface.

Lattice	Method	Bond Energy (eV)	Bond Length (Å)
111	B3LYP	−0.3	2.17
	B3LYP‐D3	−0.8	2.14
	B97‐D3	−0.9	2.14
100	B3LYP	−0.3	2.37
	B3LYP‐D3	−0.7	2.25
	B97‐D3	−0.9	2.13

Table 11 lists the principal vibrational modes of NH3 in gas‐phase and Ir‐adsorbed phase.[ , , ] After scaling, all gas‐phase vibrations for all three methods come within ∼20 cm−1of their true value, with the exception of the B97D3 estimation of the wag. After adsorption the stretch vibrations are red‐shifted ∼30 cm−1, and the wag is blue‐shifted ∼150‐200 cm−1. Other theoretical work comes closest to B97‐D3 predictions for all except the weak linear stretch between NH3 and the surface, for which B3LYP comes closest. The wag is a strong IR mode. The symmetric modes are strongly Raman active (very strong on the B3LYP (100) cluster) and all other modes are weakly Raman active.

Table 11

Vibrational modes and their corresponding frequencies (in cm‐1) for free and adsorbed NH3 at its lowest energy position – top for the (111) and (100) surface. Values are compared to existing experimental and theoretical observation. Asymmetric and symmetric stretches are abbreviated as a. stretch and s. stretch.

	Vibrational mode	Gas Phase (cm−1) [†]				Adsorbed Phase (cm−1) [†]
		B3LYP	B3LYP‐D3	B97‐D3	Other Work	B3LYP	B3LYP‐D3	B97‐D3	Other Work
	a. stretch (H2‐N‐H)	3605	3605	3538	3414[51][ac], 3444[52][ac]	3572, 3562	3575, 3562	3500, 3512	3427[40][b]
	s. stretch (N‐H)	3479	3479	3406	3337[51,52][ac]	3441	3441	3370	3333[40][b]
	bend (N‐H)	1670	1669	1651	1627[52][ac], 1628[51][ac]	1635	1633	1604	1615[53][ac*]
	wag (N‐H)	1007	1007	1016	950[51,52][ac]	1200	1201	1175	960[40][b], 1250[53][ac*]
	stretch⊥ (Ir‐N)					352	395	398	264[40][b]
	a. stretch (H‐N‐H)					3558	3585	–
	a. stretch (H2‐N‐H)					3583	3575	3495, 3489
	s. stretch (N‐H)					3448	3447	3362
	wag (N‐H)					1104	1149	1177
	rocking (N‐H)					527	600	699
	stretch⊥ (Ir‐N)					242	242	427

[a] experimental; [b] theoretical; [c] IR; [*] index plane not specified; [†] tabulated frequencies from this work are unscaled

–

Conclusions

This study had two objectives. The first was to calculate local adsorption characteristics – binding sites, bond distances, adsorption energies, and vibrational spectra for select intermediates involved in the ammonia synthesis reaction on pure iridium. These intermediates are N2, N, NH, NH2, and NH3. Secondly, three DFT methods – B3LYP, B3LYP‐D3, and B97‐D3 – were compared to determine the effect dispersion has on the popular B3LYP method and contrast it with a faster semi‐empirical GGA, B97‐D3.

B3LYP vs. B3LYP‐D3 (effect of dispersion)

The dispersive function's primary effect is to account for long‐range forces in molecular adsorption on catalyst surfaces by the inclusion of van der Waals forces, as evidenced in Figure 4. Using Grimme's dispersion with B3LYP amounts to an increase in binding energy of about 0.5 eV for several of the intermediates. This not only brings the energy closer to experimental benchmarks but is necessary for the primary reactant in ammonia synthesis, N2, to be adsorbed on the modeled catalyst. Without it, DFT optimizations do not converge and thus cannot be included in a mechanistic study involving N2 as a reactant.

Aside from energetic differences and very small decreases in bond length, the dispersion‐corrected B3LYP has a very similar electronic structure as B3LYP. Both methods predict the same ground state spin multiplicity and very similar spin density profiles on the (111) and (100) surfaces before and after adsorption of intermediates. This is consistent with Grimme's claim that the dispersive term does not rely on nor affect electronic structure. Select vibrational modes of adsorbed molecules are affected by dispersion, but most are not. Unadsorbed molecular vibrational modes are not affected at all.

B3LYP‐D3 vs. B97‐D3 (effect of hybrid and semi‐local GGAs)

B97‐D3 serves to increase the binding strength of the system and of adsorbates more than B3LYP‐D3 because GGA bonded interactions tend to be stronger than their hybrid counterparts that include exact exchange. This paper shows a decrease in adsorption energy from B3LYP‐D3 to B97‐D3 of ∼0.2 eV. Bond lengths with B97‐D3 are slightly smaller for N2, NH, NH2, and NH3, and larger for N. Gas‐phase vibrational frequencies with B97‐D3 are closer to experimental benchmarks as well as theoretical benchmarks of non‐empirical GGAs. B3LYP‐D3 overestimates the frequency of most modes, sometimes by as much as 100–200 cm−1, as in the case of the N2 stretch. Adsorbed to the catalyst, B97‐D3 is the only method that correctly predicts the lowest energy bonding pattern and orientation of N2.[ , ] Despite both functionals including dispersion corrections, both continue to under bind dispersion‐bound N2.

A notable difference in the ground state multiplicity and spin density profile of the catalyst surface exists between B97‐D3 and the hybrid methods (figures in Supplementary Data). Despite this fact, after adsorption of an intermediate to the catalyst surface, the resulting spin density profiles follow an almost identical trend as the hybrid methods. This may imply B97‐D3 has more difficulty predicting the correct ground state structure before adsorption since the surface is a charge‐transfer complex with unpaired electrons, two characteristics often modeled better by methods with exact exchange (although dispersion helps). Since the density profiles are a match after adsorption, all three methods result in the same site preferences across all adsorbates (with the exception of N2 and NH(100) using B3LYP‐D3). This fact, combined with more accurate thermodynamic data and reduced computation time, make B97‐D3 a viable option for continued study of this system and related metal cluster systems. B97‐D3 includes minimal empiricism, sitting in the top 10 least‐parameterized semi‐empirical methods, according to one study of seventy functionals.

(111) vs. (100) (effect of lattice plane)

Adsorbates on the (100) surface are more strongly bonded than on the (111) surface, except for NH3, where the bonding energies are comparable. This is due to the difference in spin densities between the adsorption sites. For NH, adsorption on (100) appears to form only two bonds, whereas it forms three on (111). However, on (100), NH is bonding with two Ir atoms with high spin densities (atoms 5 and 8 in Figure 2 ‐ B97D3), whereas on (111) it bonds with relatively lower energy Ir atoms (atoms 5,6, & 2 in Figure 2 – B97D3). The result is that the bond formation for NH reduces the spin density of the bonded atoms further on (100) than on (111), more readily stabilizing the (100) surface. For NH, the reduction is 3.2 unpaired electrons on (100), and only 2.3 on (111)(see supplementary data). The same analysis follows for N adsorption. For NH2, not only does the spin density of the bonded atoms stabilize further on (100), it also adsorbs closer to the surface, therefore having a more pronounced effect on the stronger (100) adsorption energy. The hollow site on (100) is not a preferred binding site for any of the intermediates. The spacing between Ir atoms is too great in this region for a minimum energy to be found in the overlap between the orbitals from N and those from the four Ir nearest neighbors in the hollow position. Despite their increased binding energy, Ir‐N and Ir‐NH bond lengths are greater on the (100) surface. Frequency calculations result in a ∼5–150 cm−1 difference between vibrations of the same mode on different lattice planes.

Platinum vs. Iridium (effect of catalyst)

Previous calculations using B3LYP on a 15‐atom (111) cluster of platinum are directly comparable to the B3LYP results in this study.

In the ground state of the platinum cluster, there are 12 unpaired electrons. In this study, there are 27, in addition to there being a higher cohesive energy of Ir. Preferential catalyst positions are a close match, the only difference being this study shows a small preference for adsorption of N and NH at the HCP site over the FCC site.

Binding strength increases with the same trend NH3

Future Work

In future work, it would be beneficial to move forward with a single DFT method and perform calculations for N2H, N2H2, N2H3, and N2H4, to evaluate the dissociative and associative pathways for the electrochemical ammonia synthesis process, applying a voltage and taking coverage effects into account. Experimental follow‐up could be performed with vibrational analysis.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.