Free energy of reaction by density functional theory: oxidative addition of ammonia by an iridium complex with PCP pincer ligands

Abstract

The Gibbs energy of reaction of oxidative addition of ammonia to an iridium complex in diethyl ether was calculated by seven density functional methods, in particular B3LYP, PBE, CAM-B3LYP, M05, M06, M06-L, and ωB97X. The calculated free energies, based on geometry optimization and frequency calculations in both the gas phase and solution, were compared with the experimental result, −1.3 kcal mol⁻¹, obtained by Hartwig and coworkers. The M06-L method gives the best result: −1.4 kcal/mol.

---

1. Introduction

Computational investigation of catalytic processes has advanced greatly in the last 10–15 years, and now density functional theory is widely used for interpreting mechanisms.^1,2 However many approximate density functionals in wide use are not quantitatively reliable for complex processes.³ This is sometimes ameliorated by adding molecular mechanics dispersion corrections, as in a recent study of rhodium pincer complexes as catalysts for hydroamination.⁴ In the present article (except for a brief appendix) we consider only density functionals without empirical molecular mechanics corrections, and we provide a quantitative test of modern^5–9 and popular^10,11 density functionals for the free energy of an oxidation addition reaction¹² at Ir (a sixth-period transition metal) with pincer ligands. The calculation of the free energy, ΔG⁰₂₉₈ is important because it allows comparison to the experimentally measured equilibrium constant, K, by

K = exp(−ΔG⁰₂₉₈/RT) (1)

Scheme 1

In 2005 Hartwig and coworkers published an experimental study of the oxidative addition of ammonia to an iridium complex in which 1-pentene is replaced NH₂ and H ligands;¹² see Scheme 1. They reported that the equilibrium constant K equals 9, which yields ΔG⁰₂₉₈ = −1.3 kcal/mol. This reaction is important because it activates the strong N–H bond. An interesting question is whether quantum theory has advanced to the point where such numbers can be calculated with useful accuracy for the rational design of catalysts.

Density functional theory (DFT) is the most promising quantum chemical method for such work; however, as mentioned above, many approximate density functionals are not accurate enough. During last six years new DFT methods have been presented and shown to provide higher accuracy for calculations on catalytic processes and models of catalytic systems,^7,13–18 including catalytic systems with up to 120 atoms.¹³ In this paper we investigate the performance of several DFT methods, two old ones, B3LYP¹⁰ and PBE,¹¹ and five newer functionals, CAM-B3LYP,⁹ M05,⁵ M06-L,⁶ M06,⁷ and ωB97X.⁸ Three of these density functionals (M05, M06-L, and M06) were developed in our group and are called Minnesota density functionals. The other four functionals are chosen for study because they yield useful accuracy for a large number of other problems. Because molecule 1 has 84 atoms, and molecule 2 has 73 atoms, the present application is a challenging test of the ability of modern density functionals to treat large and complex catalytic systems.

2. Methods

Since there are two reactants R and two products P, we calculated the free energy of reaction in the liquid (liq) aswhere, for a given reactant or product X, we have

G⁰₂₉₈(X,gas) = E(X,gas) + G⁰_trans(X,gas) + G_rot(X,gas) + G_vib(X,gas) + G_con(X,gas) (3)

in the gas phase and

G⁰₂₉₈(X,liq) = G⁰₂₉₈(X,gas) + ΔG⁰_S(X) (4)

in the liquid, were E is the Born–Oppenheimer energy at the equilibrium structure, G⁰_trans, G_rot, G_vib, and G_con are respectively the translational, rotational, vibrational, and conformational free energies, and ΔG⁰_S is the free energy of solvation. We assume that the conformational free energy due to the alkyl groups cancels between reactants and products, so we omit G_con. The calculations of the other terms are explained next.

The equilibrium structure and its Born–Oppenheimer energy were calculated by density functional theory with one of the seven density functionals mentioned in the introduction and with the def2-TZVP basis set from the library¹⁹ of TURBOMOLE basis sets. The translational and rotational free energies were calculated by the standard formulas, assuming rigid rotation. The gas-phase vibrational contribution was calculated for the lowest-energy conformer by the quasiharmonic approximation, which assumes the harmonic oscillator formula for the partition function but with effective frequencies.

The effective frequencies in the gas phase were calculated by a three-step process. First we calculated the harmonic frequencies by DFT. Then we scaled the frequencies by standard scaling factors^20,21 of 0.985, 1.011, 0.977, 0.982, 0.976, 0.970 and 0.976 for B3LYP, PBE, M05, M06, M06-L, ωB97X, and CAM-B3LYP, respectively, to account for anharmonicity and systematic errors in high-frequency modes. Finally, since the lowest-energy conformers of molecules 1 and 2 each have several (6 to 10) frequencies below 100 cm⁻¹, and since the harmonic approximation greatly overestimates the partition function for low-frequency modes, we raised all low frequencies to 100 cm⁻¹.

We approximate the free energy of solvation by

ΔG⁰_S(X) = ΔG_ENP + G_CDS + ΔG⁰_conc + ΔG_vib (5)

where ΔG_ENP is the electronic-nuclear-polarization free energy contribution, G_CDS is the cavity–dispersion–solvent-structure (CDS) term, ΔG⁰_conc accounts for the change in concentration on passing from the gas-phase to solution, and ΔG_vib is the change in vibrational free energy on passing from the gas-phase to solution. Because the number of moles does not change in the reaction under study, ΔG⁰_conc cancels out in ΔG⁰₂₉₈ and we need not be concerned with it. To calculate the other terms we first optimized the reactants and products in solution with the SMD²² solvation model and calculated solution-phase harmonic vibrational frequencies.

The ΔG_ENP term is based on the bulk dielectric properties of the solvent, and it includes the unfavorable change in solute internal energy that accompanies orbital relaxation in the presence of solvent (the E contribution, where E denotes electronic), the favorable change in free energy due to geometry relaxation in solution (the N contribution, where N denotes nuclear), and the net contribution from the favorable solvent polarization energy minus the cost of polarizing the solvent (the P contribution, where P denotes polarization). These terms were calculated with the SMD solvation model, whose parameters are called intrinsic Coulomb radii. We used standard SMD radii for solutes in diethyl ether, which are 1.42 Å for H, 1.85 Å for C, 1.89 Å for N, and 2.12 Å for P. There is no SMD parameter for Ir so we used the universal force field²³ (UFF) value of 1.420 Å, which is the default in the Gaussian program.

The CDS term accounts for first-solvation-shell effects. The SMD model parameters are van der Waals radii (taken from Bondi²⁴ for common elements: 1.20 Å for H, 1.70 Å for C, 1.55 Å for N, and 1.80 Å for P), an effective solvent radius (taken as 0.4 Å^22,25), and a set of geometry-dependent and solvent-dependent atomic surface tensions. No parameters were available for Ir so we took a radius of 2.0 Å and a surface tension of zero. The atomic surfaces tensions are taken from the SMD parameterization except for Ir. The lack of validated parameters for Ir is not serious for this problem because the solvent-accessible surface area²⁶ (SASA) of Ir is very small. In particular, based on the CDS and solvent radii given above, the total molecular SASAs of 1 and 2 are 521.4 and 470.8 Ų, respectively, but only 1.49 and 1.59 Ų of these SASAs are on Ir (these values are only 2.1% of SASA of a bare iridium atom). Since the SASA of Ir changes only 0.10 Ų, the contribution of the surface tension terms of Ir to the free energy of reaction is negligible.

The experiment was carried out with excess ammonia, which leads to additional hydrogen bonding possibilities²⁷ as compared to diethyl ether. The present calculations treat the solvent as diethyl ether without excess ammonia, and it is possible that this affects the extent of the agreement or disagreement of the various theoretical calculations with experiment.

We calculated ΔG_vib as ΔG_vib in solution minus ΔG_vib in the gas. We calculated ΔG_vib in solution by the same quasiharmonic approximation as used in the gas phase, but with solution-phase frequencies.

3. Results and discussion

The calculated free energies of reaction without and with solvation are given in the fourth and fifth two rows of Table 1. All seven methods predict that the solvent shifts the equilibrium toward products, with ΔG⁰₂₉₈ becoming less positive or more negative by an amount ranging from 0.6 to 1.8 kcal/mol. The M06-L density functional shows the best agreement with experiment, astonishingly differing from experiment by only 0.1 kcal/mol. We cannot expect such good agreement in general, but it is nevertheless encouraging, and it is consistent with the fact that M06-L often has the best performance for transition-metal systems. The seven methods are arranged in order of their extent of agreement with experiment, with B3LYP and CAM-B3LYP performing worst. Although B3LYP is very popular, we are no longer surprised by its poor performance because it is often very inaccurate, especially for mid-sized and large systems. One noteworthy result is that M05 has a smaller error than M06. Although the M06 functional has better average performance than M05,⁷ it is not surprising that M05 sometimes performs better for a specific system, especially one containing a transition metal, because M05 is sometimes better than M06 for systems with large near-degeneracy correlation effects.

Table 1 Free energy of reaction and its components (kcal mol⁻¹)

	M06-L	M05	M06	PBE	ωB97X	B3LYP	CAM-B3LYP
ΔE(gas)	7.0	8.1	10.6	−1.5	0.1	−6.1	−5.9
ΔH^0298(gas)	4.7	5.6	8.4	−3.6	−2.5	−8.2	−8.1
−TΔS^0298(gas)	−4.3	−3.9	−4.1	−3.8	−4.4	−3.7	−4.0
ΔG^0298(gas)	0.4	1.7	4.3	−7.5	−6.9	−11.9	−12.1
ΔG^0298(liq)	−1.4	0.4	3.0	−8.1	−8.7	−12.9	−13.2
|Error|	0.1	1.7	4.3	6.8	7.4	11.6	11.9

---

The second and third rows of Table 1 show the breakdown of the gas-phase free energy of reaction into enthalpic and entropic components, and the first row gives the energy of reaction based on Born–Oppenheimer energies ΔE. We see that the difference between the predictions of the various density functionals is primarily due to the differences in ΔE and hence in the enthalpic component rather than the entropic. In fact the difference of ΔH⁰₂₉₈ from ΔE is in the range 2.1–2.5 kcal mol⁻¹ in all cases, and TΔS⁰₂₉₈ is in the range 3.7–4.4 kcal mol⁻¹ in all cases.

4. Conclusions

The calculations were carried out for both the gas-phase reaction and the reaction in diethyl ether. Our calculations show that solvation lowers the free energy of reaction, and the resulting free energies of reaction calculated with the Minnesota density functionals M06-L, M05, and M06 agree with experiment within 0.1, 2, and 4 kcal/mol. The B3LYP, PBE, CAM-B3LYP, and ωB97X density functionals misestimate the Gibbs energy of reaction by 6.8 kcal mol⁻¹ or more. The good performance of the Minnesota functionals on this problem correlates with good performance on a number of other catalytic systems,^2,13–18 and the methods employed successfully here may be recommended for application to general problems in catalysis.

5. Appendix

A reviewer asked us to check the performance of the B3LYP-D²⁸ method, which attempts to improve the results of B3LYP calculations by adding parametrized molecular mechanical damped dispersion-like corrections, added in a post-SCF fashion. Although that is beyond the original objective of our work, we present work along the lines of this request in the present appendix. Unfortunately, B3LYP-D, like most such parametrized models, is not defined for Ir. However, recently Grimme²⁹ has optimized post-SCF dispersion-like parameters (called D3 parameters) for most of the periodic table for several density functionals, with different parameters for each functional to make up, in an average way, for that part of the attractive noncovalent interaction energy that is not properly included in the density functional under consideration. Thus the parameters make a smaller correction for functionals like M06 that already include a large part of the attractive noncovalent interaction energy in the correlation functional. D3 parameters are available for six of the seven functionals tested here. Although there is no program available that will carry out full SMD free energy calculations with these D3-corrected functionals, we can get a good estimate of whether and how much the D3 corrections improve the predictions by single-point calculations of E(X,gas) because, as discussed in section 3, at the level of accuracy involved in the present study the major variation in ΔG⁰₂₉₈ between functionals comes from the variations in ΔE. (Furthermore, we tested that the dispersion contribution to ΔE is not sensitive to differences among the various geometries optimized in the present work.)

The first row of Table 2 shows calculations of ΔE with single-point D3 corrections added, in each case using geometries and parameters^29,30 optimized without dispersion for the given functional. The difference between the first row of Table 2 and the first row of Table 1 is then used to as a correction to estimate ΔG⁰₂₉₈ for that method, and this estimate is shown in the second row. The absolute value of the error based on the second row is shown in the third row. We see that the molecular mechanics terms decrease the accuracy of the three most accurate functionals, where the terms range from 1.3 to 3.7 kcal mol⁻¹, although they improve the accuracy of the less accurate functionals, where they range from 5.1 to 7.2 kcal/mol. The large role of attractive noncovalent interactions inorganometallic chemistry has been emphasized previously,^13a,17 and it is encouraging that the better performing functionals model these reasonably well⁷ without including post-SCF corrections. The improvement in the less accurate functionals by adding such corrections is also encouraging.

Table 2 Free energy of reaction and its components (kcal mol⁻¹), including post-SCF damped dispersion terms

	M06-L	M05	M06	PBE	B3LYP	CAM-B3LYP
a Estimated as explained in the Appendix. b N.A. denotes parameters not available.
Original D3 method
ΔE(gas)	8.3	11.8	13.1	3.6	1.1	−0.4
ΔG^0298(liq)^a	0.0	4.1	5.5	−3.0	−5.8	−7.6
|error|	1.3	5.4	6.8	1.7	4.5	6.3
D3BJ method
ΔE(gas)	N.A.^b	N.A.	N.A.	4.1	1.9	0.3
ΔG^0298(liq)^a	N.A.	N.A.	N.A.	−2.5	−4.9	−7.0
|error|	N.A.	N.A.	N.A.	1.2	3.6	5.7

---

Very recently an alternative form of the damping, labeled D3BJ, has been proposed, and parameters are available for three of the functionals studied here.^30,31 Results obtained with the new damping are shown in the final three rows of Table 2, which shows that the D3BJ parametrization improves the accuracy of the correction for the three cases where it is available.

Acknowledgements

This work was supported in part by the Air Force Office of Scientific Research (AFOSR).

References

1. Transition State Modeling for Catalysis, ACS Symposium Series 721, ed. D. G. Truhlar and K. Morokuma, American Chemical Society, Washington, D.C., 1999.
2. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757.
3. J. N. Harvey, Annu. Rep. Prog. Chem., Sect. C, 2006, 102, 203.
4. A. Uwe, M. Hölscher and W. Leitner, Chem.–Eur. J., 2010, 16, 9203.
5. Y. Zhao, N. E. Schultz and D. G. Truhlar, J. Chem. Phys., 2005, 123, 161103.
6. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.
7. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215.
8. J. D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
9. T. Yanai, D. Tew and N. Handy, Chem. Phys. Lett., 2004, 393, 51.
10. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623.
11. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
12. J. Zhao, A. S. Goldman and J. F. Hartwig, Science, 2005, 307, 1080.
13. (a) Y. Zhao and D. G. Truhlar, Org. Lett., 2007, 9, 1967; (b) Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2009, 5, 324; (c) H.-C. Yang, Y.-C. Huang, Y.-K. Lan, T.-Y. Luh, Y. Zhao and D. G. Truhlar, Organometallics, 2011, 30, 4196.
14. (a) E. A. Amin and D. G. Truhlar, J. Chem. Theory Comput., 2008, 4, 75; (b) A. Sorkin, E. A. Amin and D. G. Truhlar, J. Chem. Theory Comput., 2009, 5, 1254.
15. Y. Zhao and D. G. Truhlar, J. Phys. Chem. C, 2008, 112, 6860.
16. K. Yang, J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2010, 132, 164117.
17. B. B. Averkiev, Y. Zhao and D. G. Truhlar, J. Mol. Catal. A: Chem., 2010, 324, 80.
18. (a) Y. Zhao and D. G. Truhlar, Rev. Mineral. Geochem., 2010, 71, 19; (b) Y. Zhao and D. G. Truhlar, Chem. Phys. Lett., 2011, 502, 1.
19. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
20. I. M. Alecu, J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2010, 6, 2872.
21. J. Zheng, I. M. Alecu, B. J. Lynch, Y. Zhao and D. G. Truhlar, Database of Frequency Scale Factors for Electronic Structure Methods. University of Minnesota, Minneapolis, 2010, http://comp.chem.umn.edu/freqscale/index.html (accessed June 16, 2011).
22. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378.
23. A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard and W. M. Skiff, J. Am. Chem. Soc., 1992, 114, 10024.
24. A. Bondi, J. Phys. Chem., 1964, 68, 441.
25. J. D. Thompson, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. A, 2004, 108, 6532.
26. D. A. Liotard, G. D. Hawkins, G. C. Lynch, C. J. Cramer and D. G. Truhlar, J. Comput. Chem., 1995, 16, 422.
27. L. Brammer, Dalton Trans., 2003, 3145–3157.
28. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
30. http://toc.uni-muenster.de/DFTD3/index.html (accessed Aug. 22, 2011).
31. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456.

---