FULL PAPER

WWW.C-CHEM.ORG

A Balanced Procedure for the Treatment of Cluster–Ligand Interactions on Gold Phosphine Systems in Catalysis Doreen Mollenhauer*[a,b] and Nicola Gaston[b] Ligand-protected metal clusters are difficult to describe within density functional theory due to the need to treat the electronic structure of the cluster, possible charge transfer between the ligands and the cluster, and weak ligand–ligand interactions in a balanced manner. We demonstrate the use of an appropriate, stepwise benchmarking process that accounts for the nonadditivity of these different contributions to stability and catalytic activity. We consider both open- and closed-shell clusters, differently charged systems, and ligands of increasing complexity for gold phosphine systems. The use of a dispersion correction to density functional calculations was found to be

crucial for both structure optimization and the calculation of binding energies. We find that PBE-D3 performs well with a variation in energetics of 0.7–10.9 kcal/mol, PBE0-D3 better with 0.0–3.3 kcal/mol, and B2PLYP-D3 the best with 0.2–2.4 kcal/mol, when compared to the best available benchmark [CCSD(T) or SCS-MP2]. Our systematic procedure clarifies that these functionals all give accurate results for certain cases, but for the total performance over a range of interactions, they perform in C 2014 Wiley Periodicals, Inc. accordance with Jacob’s ladder. V

Introduction

that relate their physical properties to their binding energies (BE) found that the local density approximation (LDA) yields the best agreement with the experimental structure in comparison to generalized gradient approximation (GGA; OPBE), meta-GGA (TPSS), hybrid (PBE0), and long-range corrected (LCxPBE, CAM-B3LYP) DFT functionals in combination with LANL2DZ and TZVP basis sets.[7] The initial growth mechanisms of gold-phosphine clusters have been studied using BP86 or BP86-D in combination with the TZVP basis set of Ahlrichs and a [1s224f[14]] frozen core for gold and a [1s222p[6]] frozen core for the Element P.[8] The authors noticed that the dispersion parameters for gold were not well-tested and applied dispersion correction to the small systems only. Furthermore, aurophilic interaction energies have been compared to wavefunction-based Mïller–Plesset perturbation theory (MP2) benchmark values using the same basis set and pseudopotentials (PPs), resulting in a positive interaction at the DFT level and an overestimation of binding by 130% by DFT-D. Goel et al. have utilized SVWN5, TPSS, PBE, CAM-B3LYP, and B3LYP (TZVP and LANL2DZ) to describe the BE of gold clusters containing four atoms and a phosphine shell,[9] whereas Kilmartin et al. describe the removal of the ligands from

Ligand-stabilized gold nanoparticles have been found to show many interesting electronic and optical properties.[1,2] These properties show promise for applications in fields as diverse as catalysis, photonics, nanoelectronics, biolabeling, and sensing.[1–4] In recent years, significant progress has been made in the synthesis and characterization of monodisperse small gold nanoparticles or clusters (sized below 3 nm). These small systems exhibit properties which are different from their larger counterparts. For example, they have been shown to be effective catalysts for various reactions such as the oxidation and hydrogenation of alkenes, for the low temperature CO to CO2 oxidation or the water gas shift reaction.[3–5] Especially impressive is their performance in catalysis at low temperatures, which enables the possibility of both energy savings and efficiency due to high selectivity for certain reactions. However, due to the high surface area of the clusters, the influence of the protecting ligands is very important. The ligand shell strongly affects the electronic and optical properties which makes understanding essential step to enable controlled finetuning of properties. Typical ligands for small gold clusters in catalysis include thiolate or phosphine groups and halides. There are various X-ray crystal structures that have been determined, as summarized in a recent review.[6] For gold phosphine systems, in particular, a number of X-ray structures confirm that phosphine ligands bind to only one gold atom in the cluster. This is in contrast to gold thiolate systems, for example. Electronic structure methods such as density functional theory (DFT) have found significant applications in understanding the structure, reactivity, and properties of gold phosphine clusters;[6] however, a range of different DFT functionals have been applied to goldphosphine systems. Phosphine-stabilized Au4 cluster studies 986

Journal of Computational Chemistry 2014, 35, 986–997

DOI: 10.1002/jcc.23578

[a] D. Mollenhauer Callaghan Innovation, 69 Gracefield Road, Lower Hutt, 5010 Wellington, New Zealand E-mail: [email protected] [b] D. Mollenhauer, N. Gaston School of Chemical and Physical Sciences, The MacDiarmid Institute for Advanced Materials and Nanotechnology, Victoria University of Wellington, P.O. Box 600, 6140 Wellington, New Zealand Contract grant sponsor: MacDiarmid Institute for Advanced Materials and Nanotechnology (D.M. and N.G.); Ministry of Science and Innovation’s Research Infrastructure programme C 2014 Wiley Periodicals, Inc. V

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Au6(PPh3)621 with PBE method.[10] The different methodologies used to date lead to a difficulty in comparing and evaluating the results of different computational studies. Our goal is to study the influence of phosphine ligand shells on the electronic and catalytic properties of synthesized Aun clusters (n  6). Therefore, accurate values of the BE are important to describe the removal or partial removal of ligands during catalysis. Furthermore, the number of ligands necessary to stabilize certain gold clusters can be understood in general terms from this information. Usually, triphenylphosphine ligands are chosen to stabilize these gold clusters-so long-range interactions can be important for stability, both within the ligand shell and the ligand–cluster interactions. As these larger ligand-cluster systems are usually tractable only with DFT methods, benchmark studies are useful to determine and evaluate the performance of DFT, as well as the importance of dispersion. A fundamental difficulty in the benchmarking of ligand-surrounded clusters arises because these systems include closed- and open-shell cases, as well as systems of different charges. Furthermore, dispersion corrections are more important for large ligands, whereas benchmark calculations are typically limited to small model systems due to the computational limitations of accurate wavefunction-based correlation methods. To consider all these cases, we suggest a systematic methodology for the benchmarking of these systems in three steps. This is created to focus on Au-P bond lengths (BLs) and BE; however, Au-Au BLs are not treated. The first step includes the calculation of very small systems, namely Au-PH3 and Au2-PH3. These two systems represent the basic interaction in a complex (openshell) and in a small cluster (closed-shell). Here, the benchmark method coupled-cluster singles doubles [CCSD(T)] can be applied to determine accurate structural Au-P distances and BE. We compare these results to various DFT functionals, with and without dispersion corrections. The best functionals for these small systems, with respect to both structural parameters and BE, are chosen for further examination. Additionally, less computationally intensive wavefunction-based methods MP2 and spin component-scaled version (SCS-MP2) are considered in comparison to CCSD(T) for their use as reference methods for larger model systems. Second, Au-PR3 and Au2-PR3 with R5H, Me, Ph have been selected as larger test systems. As the charge on the gold clusters is important and may vary within the catalytic cycle, the same systems are examined also when positively charged. DFT(-D) results are compared to wavefunction-based methods such as SCS-MP2 and CCSD(T) when possible as well as the different functionals to each other. The importance of dispersion correction is tested in structures with more complex ligands. Third, after this test of small neutral and charged systems, DFT(-D) performance is tested for a larger gold cluster surrounded by phosphine ligands, namely Au6(PPh3)621. Thereby PPh3 is also replaced by smaller model ligands: PMe3 and PH3. Thus, the performance of DFT(-D) is analyzed for the ligand-cluster and ligand–ligand interaction in comparison to SCS-MP2 and different functionals to each other. In addition, the charge distribution using different schemes and the highest occupied

molecular orbital HOMO-LUMO gap (HOMO: highest occupied molecular orbital, LUMO: lowest occupied molecular orbital) are considered as aspects of the electronic structure important for catalysis. The goal of this three step process of benchmarking is to test the performance of DFT for gold phosphine clusters and their different components, with increasing system size, open-, and closed-shell systems and for different charged states.

Methodology The quantum chemical calculations are performed using MOLPRO, TURBOMOLE, NWCHEM, and Gaussian09 software.[11–14] The elements P, C, and H were characterized in an all electron description by the correlation consistent valence [(aug)-ccpVXZ] Dunning-type basis sets.[15,16] An energy-consistent scalar relativistic PP for simulating the [Kr] 4d[10] core and the corresponding basis sets [(aug)-cc-pVXZ-PP] was utilized for Au.[17,18] The following wavefunction-based methods have been used: CCSD method as well as with perturbative triples corrections [CCSD(T)],[19–22] second-order MP2 plus the SCSMP2,[23,24] and the Hartree-Fock (HF) method.[25] Different DFT approximations have been used, namely: LDA[26]; the GGA, including PW91,[27] PBE,[28] PBErev,[29] and BP86 functionals[30,31]; meta-GGA functionals TPSS[32] and M06-L[33]; hybrid functionals: B3LYP,[34,35] PBE0,[36] and B3BP86[31,34]; hybrid-meta functionals M06,[37] M06-HF,[38,39] M06-2x[37] as recommended for transition metal phosphine complexes[40] as well as M05[41] and M05-2x[42]; long-range corrected functionals: CAMB3LYP,[43] LC-BLYP,[44] LC-PBE0,[44] and LC-PBE[44] as well as the double-hybrid functional B2PLYP.[45] Open-shell systems were calculated with spin-restricted theory for wavefunction-based methods and spin-unrestricted theory for DFT calculations. The spin contamination deviation was always less than 0.006 except in single cases with 0.02. All closed-shell systems were considered using spin-restricted closed-shell theory instead. Several calculations were repeated (both wf-based and DFT) and result in nearly the same BE and binding distance (deviation less than 0.2 kcal/mol) when employing unrestricted theory on closed-shell systems and spin-unrestricted (wf-based) as well as spin-restricted (DFT) theory for open-shell systems. This behavior is consistent with previous reports for a similar system, namely pyridine gold complexes.[46] Thus, we assume that the different results obtained with spin-unrestricted and restricted theory are comparable with small differences. In part A, the structure of phosphine is kept fixed at the experimental one adapted from Ref. [47]. Potential energy curves varying the Au phosphine dis˚ have been calculated. To detertance between 1.8 and 3.6 A mine the Au-phosphine distances and BE, the calculated potential energy curves are fitted to polynomials up to sixth order. BEs at every single point have been corrected for the basis set superposition error (BSSE) by a counterpoise correction.[48] Basis set extrapolations have been performed for the correlation part of the energy using the formula of Ref. [48]. The BSSE for CCSD(T) and aug-cc-pVTZ(-PP)/aug-cc-pVQZ(-PP) extrapolation was in the order of 1 kcal/mol. DFT calculations Journal of Computational Chemistry 2014, 35, 986–997

987

FULL PAPER

WWW.C-CHEM.ORG

Table 1. Binding energies (BE) and equilibrium Au-P distances (BL) as obtained with different wavefunction-based approaches as well as DFT functionals for phosphine gold complex. Method CCSD(T) CCSD MP2 ROHF UHF B2PLYP PBE0 B3LYP B3P86 CAM-B3LYP LC-BLYP LC-PBE0 LC-PBE M06 M06-HF M06-2x M05 M05-2x TPSS M06-L PBErev PBE BP86 PW91 LDA

BL (A˚)

BE (kcal/mol)

2.332 2.353 2.275 nb nb 2.349 (2.348) 2.349 (2.340) 2.398 (2.391) 2.246 2.380 2.329 2.292 2.279 2.482 2.149 2.291 2.485 2.365 2.339 (2.338) 2.348 2.364 (2.353) 2.344 (2.337) 2.346 (2.335) 2.346 2.277

217.4 213.9 221.3 nb nb (216.7, (218.0, (214.3, 218.0 214.0 220.7 222.4 226.7 212.2 217.5 29.2 29.4 213.7 (220.4, 214.7 (218.2, (222.1, (221.8, 221.8 234.5

216.0 217.3 213.1

219.4 216.5 221.4 220.3

217.7) 218.7) 216.0)

221.3) 219.7) 222.9) 223.3)

have been performed using aug-cc-pVTZ(-PP) basis sets, and the dispersion correction was applied in Grimme’s D3 scheme with zero- or BJ-damping.[49,50] The BSSE was tested for several cases and found to be very small (around 0.1 kcal/mol). In part B, structure optimizations have been performed at the PBE(-D3)/cc-pVTZ(-PP) level of theory, in part C at PBE(-D3)/cc-pVDZ(-PP). Furthermore, the structural optimization in part C was performed with the multipole acceleration resolution of identity approximation for additional efficiency.[51–53] Test calculations show no difference in the structures due to this approximation. All minimum energy structures of parts B and C were verified with frequency calculations and additional single point calculations were calculated using cc-pVTZ(-PP) basis set for GGA and hybrid functionals. Using AuPPh3 as a test system indicated that the BSSE increases to 0.4 kcal/mol. For doublehybrid functionals, we used aug-cc-pVTZ(-PP). The sufficiency of these basis sets to describe the BE has been shown for a similar system, namely gold pyridine complexes.[46] Three different types of charges were calculated: Mulliken,[54] which extracts the charge information from the basis functions and shows a strong basis set dependency; the natural population analysis (NPA),[55] an improved scheme of Mulliken which relates the molecular density to the sum of atomic orbital contribution resulting in a reduced basis set dependency; and finally, Bader charges[56,57] for which the overall Journal of Computational Chemistry 2014, 35, 986–997

Results and Discussion Benchmark calculations for Au-PH3 and Au2-PH3

For CCSD, CCSD(T) and MP2 extrapolated values were given at basis aug-cc-pVTZ/QZ(-PP) with Hartree–Fock results as converged with augcc-pVQZ(-PP). For DFT calculations, the aug-cc-pVTZ(-PP) basis set was utilized. Dispersion correction was applied if possible, and the corresponding values are denoted in brackets. For the BL, the zero-damping is given, for the BE first the value derived with zero-damping than with BJ-damping. The abbreviation nb means nonbinding.

988

space is split into atomic domains and the atomic charge is derived by integrating the density within these basins. The Wiberg bond index (BI)[55] has been calculated in the NBO basis as the sum of quadratic nondiagonal elements of the density matrix between two atoms and can be interpreted as the quantum chemical analog of a bond order. Charges as well as Wiberg bond indices were in all cases calculated using PBE0/cc-pVTZ(-PP). However, test calculations with PBE give very similar results.

Benchmark calculations using wavefunction-based correlation methods were used to determine DFT functionals with and without dispersion correction for the calculation of binding distances and energies for the basic interaction. Furthermore, it was tested whether the performance for the open- and closed-shell systems is comparable. Therefore, Au-PH3 and Au2-PH3 have been chosen as benchmark systems. The BE between the ligand and the gold atom or dimer and the Au-P BL were considered and compared to benchmark values. For reliable values that enable comparison of different functionals, all methods have been used at the basis set limit. For CCSD(T) and MP2 an aug-cc-pVTZ/QZ(-PP) extrapolation was taken, while for the different DFT functionals we have used aug-ccpVTZ(-PP), which was tested to give values at the basis set limit. The results are shown in Tables 1 and 2 as well as for better comparison in Figures 1 and 2. In addition, the electronic structure has been incorporated into our analysis, by calculating the HOMO–LUMO gap (see Fig. 3 and Supporting Information, Table S1). Wavefunction-based methods. The CCSD(T) results for the phosphine gold complex and dimer (T1 never > 0.02) demonstrate that the BE is twice as strong for the phosphine dimer ˚ ) BL. The perturbative interaction, with a shorter (around 0.1 A triples contribute 20% for the complex and half as much for the dimer interaction which means a constant decrease of BE. As can be expected, no binding for the phosphine complex has been found at HF level. For the phosphine dimer interaction, the HF yields a binding contribution of 42% of the CCSD(T) BE. Therefore, the increase of BE for the system is mainly due to the HF part (see also Fig. 1) which implies a change in the binding nature. Our CCSD(T) value of the phosphine gold complex agrees well with previous calculated onecomponent (scalar) CCSD(T) with aug-cc-pVTZ basis set of 15.9 kcal/mol.[58] The MP2 method, conversely, overestimates the BE (by 20%) in both cases. However, SCS-MP2 gives improved values, rather close to the benchmark method for the phosphine gold dimer. This suggests that one should use SCS-MP2 rather than MP2 as a reference method for larger systems for which degenerate states are not occurring. Density functional theory(-D3). For the gold phosphine com-

plex, the Au-P BL can be best described by BP86-D3 and PBED3, in comparison to CCSD(T). The BL was slightly improved WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

Table 2. Binding energies (BE) and equilibrium Au-P distances (BL) as obtained with different wavefunction-based approaches as well as DFT functionals for phosphine gold dimer. Method CCSD(T) CCSD MP2 SCS-MP2 ROHF UHF B2PLYP PBE0 B3LYP B3P86 CAM-B3LYP LC-BLYP LC-PBE0 LC-PBE M06 M06-HF M06-2x M05 M05-2x TPSS M06-L PBErev PBE BP86 PW91 LDA

BL (A˚)

BE (kcal/mol)

2.255 2.267 2.203 2.229 2.397 2.394 2.273 (2.269) 2.270 (2.268) 2.296 (2.288) 2.268 2.283 2.261 2.238 2.227 2.343 2.171 2.276 2.349 2.261 2.243 (2.246) 2.245 2.281 (2.275) 2.269 (2.263) 2.272 (2.265) 2.268 2.211

237.6 234.0 246.2 239.5 216.1 216.0 (236.5, (236.5, (232.0, 236.2 233.2 241.6 243.5 247.9 228.5 240.2 227.2 224.6 233.2 (238.1, 231.0 (235.1, (239.2, (239.4, 238.6 253.1

235.4 235.4 230.1

236.4 232.8 238.2 237.1

237.4) 237.1) 233.7)

239.0) 236.7) 240.0) 240.9)

For CCSD, CCSD(T) and (SCS-)MP2 extrapolated values were given at basis aug-cc-pVTZ/QZ(-PP) with Hartree–Fock results as converged with aug-cc-pVQZ(-PP). For DFT calculations, the aug-cc-pVTZ(-PP) basis set was utilized. Dispersion correction was applied if possible, and the corresponding values are denoted in brackets. For the BL, the zerodamping is given, for the BE first the value derived with zero-damping than with BJ-damping.

by Grimme’s dispersion correction, where applied. In these cases, both zero- and BJ-damping gave very similar results. ˚ . Comparable Overall, the BL varies between 2.14 and 2.48 A results to CCSD(T) with a BE within 1 kcal/mol deviation are

Figure 1. Potential energy surfaces for the phosphine gold complex and dimer as obtained at CCSD(T)/aug-cc-pVTZ/QZ(-PP) and HF/aug-cc-pVQZ(PP) level of theory. The correlation energy is derived as difference between CCSD(T) and HF values. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FULL PAPER

obtained with B2PLYP-D3, PBE0, B3P86, and PBErev-D3 for AuPH3, with PBE0 giving the best agreement. B2PLYP, B3LYP, and PBErev yield an improvement when the dispersion correction is applied. For all other functionals, the dispersion correction leads to a stronger overestimation than the uncorrected results. In general, hybrid functionals obtain best agreement with BJ-damping, while GGA functionals do better with zerodamping. If we analyze only the values closest to the benchmark result, hybrid and double-hybrid functionals underestimate the BE, whereas standard GGA functionals overestimate it. Overall, the DFT functionals vary within 19/28 kcal/mol compared to CCSD(T). The hybrid long-range corrected functionals overestimate BE noticeably (more than 20%) and underestimate the BL compared to CCSD(T), with the exception of CAM-B3LYP where the opposite occurs. The hybridmeta functionals show an underestimation of BE (1–47%). For the phosphine gold dimer, the performance of more DFT functionals is generally better relative to CCSD(T). Nevertheless, the overall deviation of BE increases to 110/213 kcal/ mol. Again B2PLYP-D3 and PBErev-D3 result in an excellent description of the BE [within a 1 kcal/mol of CCSD(T)] and PBE0 is equally good with inclusion of the dispersion correction. Again, the dispersion correction improves the BE only for the hybrid functionals and PBErev, as was seen for the phosphine complex. Similar trends for the overestimation and underestimation of BE are found for the phosphine gold dimer as for the complex. HOMO–LUMO gap. Due to the interest in these gold phosphine systems for catalysis, a good understanding of their electronic structure is imperative. A comparison of the DFT performance for the HOMO–LUMO gap of our test systems is, therefore, given. Values have been calculated based on the CCSD(T) determined structural parameters (see Fig. 3, explicit values in Supporting Information, Table S1). HOMO–LUMO gaps derived from structures obtained with the same functional vary not more than 0.2 eV. The gap results are qualitatively very different, giving values from 2.6 to 7.5 eV for the Au-PH3 complex and 3.1–8.7 eV for the Au2-PH3. The Hybrid functionals, meta-GGA, and GGA functionals show consistent values within the functional class. Plotting the HOMO–LUMO gap obtained with different functionals of the gold phosphine complex against the gold dimer phosphine demonstrates an approximately linear dependency. This shows that the different functionals as well as Hartree–Fock calculations deliver quantitatively strongly different, yet qualitatively similar results. All functionals show similar ratios. Furthermore, open-shell restricted and unrestricted calculations result in nearly the same values for the HOMO–LUMO gap. Conclusion. As it is expected that the BE varies only slightly by calculating at the optimal BL, a method combination for bigger systems is suggested. Therefore, the best agreement for BL for both systems is given by PBE-D3 and BP86-D3. For the BE, the best performance is in order of decreasing accuracy: B2PLYP-D3, PBE0/PBE0-D3, and PBErev-D3. Therefore, a structural optimization with PBE-D3, followed by the calculation of energetic parameters with B2PLYP-D3 would be Journal of Computational Chemistry 2014, 35, 986–997

989

FULL PAPER

WWW.C-CHEM.ORG

Figure 2. Binding energy over bond length obtained at different level of theories. In cases where dispersion correction improves the results, the DFT-D3 value is given, otherwise just the DFT value. The upper graphs show all values (except LDA), the lower graphs a zoom with a BE range of 62 kcal/mol, and BL range of 61% regarding the CCSD(T) value. The values are corresponding to the Tables 1 and 2.

preferred in this case. This will be tested in the second part of this article. The HOMO–LUMO gap shows significant variation between the different functionals, but all functionals result in qualitatively similar effects when increasing the number of gold atoms.

Figure 3. The HOMO–LUMO gap is plotted as Au2-PH3 over Au-PH3 for different DFT-functionals and Hartree–Fock. The basis set applied is aug-ccpVTZ(-PP). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

990

Journal of Computational Chemistry 2014, 35, 986–997

Benchmark calculations for Au(1)-PR3 and Au2(1)-PR3 The Au-PH3 and Au2-PH3 complexes have served as the first small benchmark systems in the previous section. Experimentally, gold clusters are stabilized with larger phosphine ligands such as triphenylphosphine. We require that the DFT functional should also show a reasonable performance for this larger ligand. To analyze the behavior of DFT, the series of PH3, PMe3, and PPh3 is used in a comparative study using the complex and gold dimer interaction. We compared also the effect of changing the charge of the gold system to 11. The systems were structurally optimized using PBE both with and without dispersion correction. The influence of the dispersion correction on the different systems was analyzed in detail. The PBE-D3 optimized structures were used for single point calculations by PBE0(-D3) and B2PLYP(-D3) which gave the best results as shown previously in Introduction Section. Additionally, PBE(-D3) was considered because it was used for the structure optimization. In each case, both the zero- and BJdamping versions of the D3 correction were used and compared. If possible (for the smaller systems), we also applied CCSD(T) at the basis set limit and SCS-MP2 at or close to the basis set limit. In the latter case, we include results of a basis set dependency study in the Supporting Information. To analyze the DFT results for application to catalysis, the HOMO–

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Figure 4. Binding energies (BE) of optimized structures at PBE-D3/cc-pVTZ(-PP) level of theory: Au-PR3 and Au6PR3 (R5H, Me, Ph). CCSD(T)/aug-cc-pVTZ/ QZ, SCS-MP2/aug-cc-pVTZ/QZ, B2PLYP(-D3)/aug-cc-pVTZ and other DFT(-D3)/cc-pVTZ is applied. Dispersion correction was applied with zero- (z) or BJdamping. BEs are given in kcal/mol, charges in e, and the HOMO–LUMO gap in eV. Furthermore, the Au-P distance in A˚ as well as Wiberg bond index (BI) is presented. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

LUMO gap and different kinds of partial charges were calculated. All results are plotted in Figures 4 and 5. Tables of the data in the Figures can be found in the Supporting Information. To describe how the ligand–cluster interaction contributes to the total BE, the nature of binding was analyzed through the following measures. The Wiberg bond indices were calculated for the PH3 systems and compared to the molecular orbital diagram. The simple Au-PH3 system shows a fully occupied binding Au 5d 2z /P lone pair orbital and a singly occupied corresponding antibonding orbital. This indicates a binding order of approximately 0.5 which agrees with the Wiberg BI. Furthermore, it should be mentioned that charge transfer plays a crucial role for this interaction.[58] The Au-P BL is also relativistically shortened due to the large relativistic enhancement of the electron affinity of gold. This leads to a considerably better overlap between the valence 6s orbital of the metal, and the lone-pair orbital of the ligand than, for example, AgAPH3 or Cu-PH3.[58] Description of the nature of cluster-ligand binding.

However, for Au2-PH3 a similar bonding orbital as for the corresponding complex is fully occupied, whereas no occupied antibonding orbital exists, resulting in a bond order of around 1 which is also in agreement with the Wiberg BI. There is also an increase in the BE by a factor of 2 consistent with this interpretation. For both the related charged systems, a bond order of around 1 is evident in the MO diagram analysis and in the Wiberg bond indices, with a fully occupied binding Au 5d 2z /P lone pair orbital able to be identified. However, the charged gold phosphine systems show an even stronger increase in the BE arising from the HF part. Thus, there is an electrostatic contribution and an enhancement in charge transfer in addition to the covalent bonding contribution. The four test systems can, therefore, be interpreted as showing three different bonding situations. The influence of the dispersion correction on the structure optimization was analyzed. The influence of D3 on the Au-P distance was found to be smaller than 3% in all cases, while the Au-P-C(H) angle, crucial for overall cluster

Structural parameters.

Journal of Computational Chemistry 2014, 35, 986–997

991

FULL PAPER

WWW.C-CHEM.ORG

Figure 5. Binding energies (BE) of optimized structures at PBE-D3/cc-pVTZ(-PP) level of theory: Au2-PR3 and Au36PR3 (R5H, Me, Ph). CCSD(T)/aug-cc-pVTZ/ QZ, SCS-MP2/aug-cc-pVTZ/QZ, B2PLYP(-D3)/aug-cc-pVTZ and other DFT(-D3)/cc-pVTZ is applied. Dispersion correction was applied with zero- (z) or BJdamping. BEs are given in kcal/mol, charges in e, and the HOMO–LUMO gap in eV. For the given charges, Au(1) is the atom which binds directly to PH3. ˚ as well as Wiberg bond index (BI) is presented. [Color figure can be viewed in the online issue, Furthermore, the Au-P and Au-Au distance is given in A which is available at wileyonlinelibrary.com.]

structure, was affected by less than 0.2 . The C4–C4 distances ˚ . Thus, the disperof the PPh3 ligand changed by up to 0.3 A sion correction changes the structure only slightly, but with increasing changes along increasing complexity of PR3 ligands. This correction becomes most important for the experimentally relevant PPh3 ligands. The Au-Au distance is hardly influenced by the dispersion correction. Binding energy. Interestingly, the dispersion correction (zeroand BJ-damping) shows similar contributions to the BE for all systems, namely for systems with PH3 an increase of the BE of around 1 kcal/mol, for PMe3 of around 2–3 kcal/mol, and for PPh3 of around 4–5 kcal/mol. One exception is the dispersion correction with BJ-damping for the charged and neutral phosphine gold dimers. There the BE is increased for increasing ligand size by 3, 6, and 11 kcal/mol, respectively. The strong influence of dispersion correction on the BE demonstrates the 992

Journal of Computational Chemistry 2014, 35, 986–997

importance of the use of the dispersion correction for phosphine ligands of increasing complexity. The use of zerodamping achieves better agreement with benchmark values for small ligands and charged systems, BJ-damping for larger ligands and especially for the interaction with the neutral gold dimer. Considering the neutral gold phosphine systems in general, a sharp increase in the magnitude of the BE is observed from Au-PH3 to Au-PMe3, whereas the Au-PPh3 system is then slightly more strongly bound. For the Au-PH3 complex, PBE0D3 and B2PLYP-D3 present the best agreement with the CCSD(T) benchmark value (deviation less than 1 kcal/mol). The BE differences between the different DFT functionals remain constant for the different phosphine ligands. This picture arises also with small deviations for the phosphine gold dimer interaction. Thereby, the BE differences between the different DFT functionals decrease with increasing complexity of the ligands, WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

especially PBE and B2PLYP come closer together. SCS-MP2 results for systems with the PH3 ligand overestimate the binding in comparison to CCSD(T). The overestimation by SCS-MP2 remains for the systems with more complex ligands in comparison to the DFT values. For the charged gold phosphine systems also after a sharp increase in BE from Au1-PH3 to Au1-PMe3, a slight further increase in binding is seen for Au1-PPh3. SCS-MP2 and CCSD(T) show a very good agreement for Au1-PH3. PBE0 and B2PLYP yield without the dispersion correction a deviation of 1–2 kcal/mol from CCSD(T); however, the dispersion correction decreases the performance to 2–3 kcal/mol deviation. For the different phosphine ligands, B2PLYP-D3 and PBE0-D3 gave a consistently good agreement with SCS-MP2. PBE shows a continuous overestimation of the results by more than 10 kcal/ mol. Similar results are obtained for the charged phosphine dimer, but here PBE0-D3 and B2PLYP-D3 show already a very good agreement with CCSD(T) for PH3. Furthermore, similar ratios of the DFT values are achieved for the more complex ligands. PBE consistently overestimates energies, and increasingly so with increasing ligand size. HOMO–LUMO gap and atomic charge. The HOMO–LUMO gap decreases consistently along the series of PR3 ligands. However, the different DFT functionals produce quantitatively different values of the gap, but with systematic behavior observed for all functionals when increasing the complexity of the ligand. The ratio of the HOMO–LUMO gaps for clusters with different ligands is almost constant for all functionals for the PH3 and PMe3 ligand, with only a few exceptions for the PPh3 ligand. This demonstrates again that similar trends can be observed with all functionals for the variation of ligand complexity. In general, the PMe3 effects a stronger charge transfer for neutral systems than PPh3 as judged by the total charge difference (calculated with Mulliken, NPA, Bader) on the metal atoms. Considering the charges in detail, the outcome looks a bit more complex. For the neutral Au-PR3 complexes, Mulliken and NPA charges show a similar trend of slight decrease for PMe3 and increase for PPh3, whereas the absolute values of the Mulliken and NPA charges differ by around 0.3 e. The Bader charges are generally intermediate between them and exhibit a different trend. However, all charges are negative and indicate charge transfer from the ligand to the atom. The neutral phosphine dimers show a similar behavior, except that Mulliken and Bader charges are in better agreement. The charge for the gold atom which binds to phosphine (Au1) is decreasing going to PMe3 and increasing for PPh3 (Mulliken, Bader) as well as negative, but the NPA charge is about 0.3 e higher, positive in sign and shows a nearly constant behavior. For the other gold atom (not bound to phosphine, Au2), the results are consistently negative and show different trends. The total charge of the dimer is negative for all methods, ranging from 20.2 to 20.6 e. The positively charged complex shows a steeper decrease with the ligand complexity, with Mulliken and Bader again in good agreement and the NPA charges higher by about 0.2 e.

For the charged phosphine gold dimer systems, the charge on the gold atom bound to phosphine (Au1) decreases for the larger ligands and remains nearly constant for PMe3 and PPh3. The charge on the second gold atom (Au2) is consistently more positive. All three types of charge calculation show different charge differences for the two gold atoms. This difference remains roughly constant for NPA and Mulliken charges and increases for Bader. In all cases, the total charges are positive (0.1–0.6 e). In summary, the dispersion correction (using PBE-D3) leads to only small changes in the Au-P BL and is not very important for systems with PH3 or PMe3 in comparison to the pure functional. More significant differences can be found for the systems with PPh3, as here the p interaction between the phenyl rings requires a good description of dispersion. The BE is best described with B2PLYP-D3 (within 2 kcal/mol of benchmark values), followed by PBE0-D3 with agreement within 3 kcal/mol. PBE generally overestimates the BE, especially for charged systems. The neutral phosphine dimer is one exception, however, where PBE-D3 shows very good results in comparison to CCSD(T). Using zero instead of the BJ-damping version of the dispersion correction gave in some cases (about half ) a better agreement with benchmark values. The zerodamping performed better for small ligands and charged systems, whereas BJ damping is to be preferred for larger ligands, especially for the interaction to the gold dimer. The HOMO– LUMO gap varies considerably for different functionals, but shows systematic trends for the different considered systems. The different calculated charges (Mulliken, NPA or Bader) produce qualitatively similar trends. The absolute values of Mulliken and Bader are generally consistent, whereas NPA charges are mostly larger. Conclusions.

Benchmark calculations for a large gold phosphine system The benchmark calculations were extended to the larger Au621(PPh3)6 system. The X-ray structure from Ref. [59] served as the starting structure. The cluster structure itself is an edgesharing bitetrahedral structure, with each gold atom bound to a phosphine atom. The phenyl was replaced with methyl and hydrogen to obtain reference values. The Au621(PR3)6 system with R 5 H, Me, Ph was optimized in structure using PBE/cc-pVDZ(-PP). Furthermore, the dispersion correction D3 (zero) was utilized. There were no imaginary frequencies for Au621(PH3)6, but small ones (below 30 cm21) for the others. These were identified as translational and rotational contributions, and the structures were assigned to be minima. Frequency calculations without the dispersion correction show a decrease of the number and magnitudes (below 15 cm21) of the imaginary frequencies. The ligand shell-surrounded cluster with phenyl groups, as optimized with and without the dispersion correction is compared to the experimental structure in Table 3 and Figure 6. The average Au-Au BL of the clusters is consistently described by PBE and PBE-D3 for the PH3 ligand. For PMe3 Structural parameters.

Journal of Computational Chemistry 2014, 35, 986–997

993

FULL PAPER

WWW.C-CHEM.ORG

Table 3. Bond length (BL) of optimized structures for Au621(PR3)6 (R5H, Me, Ph) at PBE-D3/cc-pVDZ(-PP) level of theory. Au621(PH3)6 Au6[2]1(PMe3)6 Au6[2]1(PPh3)6 Exp.[59] PBE BL Average Au-P BL Average Au-Au BL min Au-Au BL max Au-Au PBE-D3 BL Average Au-P BL Average Au-Au BL min Au-Au BL max Au-Au

2.350 2.825 2.673 2.883

2.354 2.828 2.685 2.879

2.361 2.833 2.678 2.913

2.299 2.759 2.651 2.839

2.342 2.821 2.690 2.870

2.339 2.813 2.714 2.858

2.333 2.800 2.694 2.874

2.299 2.759 2.651 2.839

˚. BL are given in A

˚ . There is an and PPh3, the difference is increasing to 0.03 A increase of the average Au-Au BL with the complexity of the substituents for PBE, whereas the expected decrease related to ligand donor properties required the PBE-D3 description. For the minimal and maximal Au-Au BL in the cluster, similar trends are obtained for PBE and PBE-D3. However, the minimal BL is increased by applying the dispersion correction, whereas the maximum decreases. For the average Au-P BL, there is an increase along R 5 H, Me, and Ph with the pure PBE functional, and a decrease for PBE-D3. This shows that the dispersion correction is necessary to obtain the right substituent behavior, in agreement with the BE. Comparing these structural values with experimental parameters, PBE-D3 is best except for the minimal Au-Au distance. Using PBE-D3 resulted in an overestimation compared to ˚ , Au-Au BL s: 0.02– experimental data (Au-P BL s: 0.03–0.04 A ˚ ). Thereby the largest system with PPh3 ligands is 0.06 A ˚ . The PBE-D3 optimized structures described within 0.04 A served as starting points for the following paragraphs. However, the BL overestimation compared to experimental data using PBE functional without dispersion correction is about ˚ for Au-Au. 0.05–0.06 A˚ for Au-P BLs and 0.02–0.07 A To test the sufficiency of the basis set, a cc-pVTZ basis set was applied for R 5 H, Me. This results in a decrease of the BLs of about 0.01–0.02 A˚ in the direction of the experimental values for PH3, and a slightly weaker adjustment for PMe3. Thus, better basis sets will give an improvement in structural parameters, but packing effects in the crystal might also play a role for these small differences. In any case, the cc-pVDZ basis set produces reliable results in a realistic computational time.

It is also interesting to visually compare the overall structure to the experimental value, as shown in Figure 6D. The good agreement achieved with the dispersion correction is also demonstrated relative to the pure functional in Figures 6A–6C. The structural agreement of PBE and PBE-D3 is very good for the system with PH3, decreases for PMe3 and shows big differences for PPh3, especially in the relative positions of the phenyl rings. Binding energy. The BE of the Au621(PR3)6 systems have

been tested in three different situations. To obtain a meaningful reference value for the larger systems, the interactions were decomposed into the Au621-ligand and ligand–ligand contributions. Calculations using PBE, PBE0, B2PLYP, and the dispersion correction could then be compared to wavefunction-based SCS-MP2 values for all ligands (see Fig. 7). Furthermore, the average BE for Au-PR3 has been calculated for the total system. Here, we use SCS-MP2 only for the smallest ligands. The ligand-gold BE was determined for Au621-PR3 using the optimized structure of the full system at the PBE-D3/cc-pVDZ level of theory. The other ligands have been removed. The BE is increased compared to the average BE, as can be expected. It can be seen that applying the dispersion correction gives results closer to SCS-MP2 for all functionals. The dispersion correction increases the BE for R 5 H, to Me and Ph by approximately 3, 6, and 11 kcal/mol. Interestingly, B2PLYP-D3 and PBE-D3 show best agreement with SCS-MP2 for Au621-PH3 and Au621-PMe3 within an accuracy of 1 kcal/mol, in contrast PBE0-D3 gives the best value for Au621-PPh3, with similar accuracy. It should be mentioned that the SCS-MP2 basis set extrapolation is expected to underestimate the BE by 1–2 kcal/ mol. Overall, BJ-damping is seen to give the better agreement with benchmark values. The calculation of the ligand–ligand interaction used three ligands from the overall optimized structure, calculating the BE between one ligand and the others. As expected, there is no interaction for PH3 ligands and only a very small interaction for PMe3. However, a significant BE was found for the PPh3 ligands. All dispersion corrected DFT values agree very well with SCS-MP2 results. Neglecting the dispersion correction would effectively ignore this ligand–ligand interaction in DFT calculations. This ligand–ligand interaction between one ligand to all others is about 17 kcal/mol for PPh3. This shows that the interaction within the ligand shell contributes significantly to

Figure 6. The optimized structures of Au621(PR3)6 with R5H, Me, and Ph. Blue color corresponds to an optimization of PBE/cc-pVDZ, red color to PBE-D3/cc-, and green color to experimental structure as used from Ref. [59]. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

994

Journal of Computational Chemistry 2014, 35, 986–997

WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

FULL PAPER

Figure 7. Results of Au621(PR3)6 average bond length, the Au621(PR3) interaction as well as ligand–ligand interaction are presented on different level of theories. Structure optimizations of Au621(PR3)6 took place at PBE-D3/cc-pVDZ(-PP) level of theory. Single point calculations of BE and properties are performed using SCS-MP2/aug-cc-pVDZ/TZ (cc-pVDZ/TZ), B2PLYP/aug-cc-pVTZ (cc-pVTZ), and other DFT/cc-pVTZ is applied. BEs are given in kcal/mol.

its overall stability, and the use of dispersion correction is especially important for investigating the removal of ligands: as not only the gold-ligand bond but also the ligand–ligand interaction has to be overcome. For the full Au621-PR3 systems, the average BE is consistent across all dispersion corrected functionals, within 1.4 kcal/mol for R 5 Me and 0.2 kcal/ mol for R 5 Ph. The deviation between PBE/PBE0-D3 and B2PLYP-D3 is largest for PH3. The cc-pVDZ(-PP)/cc-pVTZ(-PP) basis set induced underestimation of the SCS-MP2 value, of roughly 2 kcal/mol, means that we achieve best agreement for Au621-PH3 with B2PLYP-D3 and PBE-D3. HOMO–LUMO gap and atomic charge. The HOMO–LUMO gap is very different in magnitude for the different functionals as seen previously, but again the trends with increasing ligand complexity are systematic (see Fig. 8). The overall charge on the gold cluster is negative, according to Mulliken and Bader,

but positive according to the NPA. Thus, similar trends are obtained for the charge transfer from the ligand shell to the cluster, as we saw previously for the small systems. In summary, it could be shown that the dispersion correction improves the structure optimization significantly, especially for the ligand phenyl rings but also for the Au-Au and Au-P distances. The dispersion correction is crucial to describe the ligand–ligand interaction in the systems with PPh3. For the Au–ligand interaction, the correction consistently improves the BE in comparison to benchmark values. All density functional results are in general agreement, and it can be concluded that B2PLYP-D3 and PBE-D3 give the best results for systems including PH3 and PMe3 ligands, and PBE0-D3 for PPh3 ones. As the system is overall a closed-shell system, a better performance can be expected for PBE-D3 than for open-shell cases. In all cases, the use of BJ-damping gives Conclusions.

Journal of Computational Chemistry 2014, 35, 986–997

995

FULL PAPER

WWW.C-CHEM.ORG

Figure 8. HOMO–LUMO gaps and partial charges of Au621(PR3)6 with R5H, Me, and Ph calculated at PBE0-D3/cc-pVTZ(-PP)//PBE-D3/cc-pVdZ(-PP). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

better performance than zero-damping for the dispersion correction of the gold–ligand interaction.

Conclusions The performance of DFT, both with and without dispersion corrections, was investigated for gold phosphine systems relevant for catalysis. Benchmark calculations were done according to a three step procedure. First, the basic interactions were investigated in small test systems, demonstrating which functionals work best for both closed- and open-shell systems, in comparison to highly correlated wavefunction-based methods. This also demonstrated that SCS-MP2 is suitable for benchmark calculations in larger systems. However, these test systems were too small to take into account the importance of dispersion correction, with the correction improving results only in certain cases. Based on these first results, three functionals of different classes were chosen for the second step, namely one GGA (PBE), one hybrid (PBE0), and one doublehybrid functional (B2PLYP). The GGA in particular was chosen as it gave best results for structural parameters. A significant detail is that the performance of PBE was dramatically improved in closed-shell systems, whereas the performance of PBE0 and B2PLYP worsen slightly. However, dispersion corrections then deteriorated the result for PBE, whereas an improvement was obtained for B2PLYP and PBE0. For the description of the HOMO–LUMO gap, different functionals produce significantly different values, but a qualitatively similar behavior was observed with increasing cluster size. Using more complex ligands for the gold complex and dimer (both neutral and charged) in our second benchmarking step, the importance of the dispersion correction could be illustrated. All functionals provide similar trends for the BE with respect to the phosphine ligands. PBE0-D3 and B2PLYPD3 agree very well in general for the BE with the variation of 996

Journal of Computational Chemistry 2014, 35, 986–997

ligands, and also agree best with benchmark values. In cases where the system was charged, the PBE functional yields a strong overestimation of the BE, with the dispersion correction even decreasing performance. The zero-damping version of D3 gave overall better agreement with benchmark results for small ligands and charged systems, but the results of BJdamping were better for larger ligands and especially for the interaction to the gold dimer. The previously observed systematic behavior of density functionals for the HOMO–LUMO gap was again observed for the increase in ligand complexity. Calculated charges based on different approaches yield mostly similar trends, but slightly different absolute values. The same functionals used in our second step were then extended to the calculation of an Au6 cluster surrounded by a phosphine ligand shell. The importance of dispersion correction was clearly demonstrated; in addition to the Au–ligand interaction, it is also crucial for the interaction within the ligand shell. It was, therefore, already necessary to include the correction for the structural optimization. The ligand interaction within the ligand shell was described very accurately by all functionals with the inclusion of dispersion correction. Concerning the gold–ligand interaction, B2PLYP-D3 and PBE-D3 give the best agreement for PH3 and PMe3 ligands, whereas PBE0-D3 is best for PPh3 ligands. In each case, BJ-damping resulted in values closer to benchmark ones. For the total (averaged) Au-ligand binding, all functionals result in similar values. The HOMO–LUMO gap and calculated charges followed the same trend as in our previous steps. For our selected DFT functionals with dispersion correction, comparison to benchmark values can be expected within: a variation of 0.7–10.9 kcal/mol (PBE-D3), 0.0–3.3 kcal/mol (PBE0D3), and 0.2–2.4 kcal/mol (B2PLYP-D3). The case of Au2PR2 is an exception, where the SCS-MP2 value is overestimated. This shows that all three dispersion corrected functionals are able to give very accurate results for certain cases, but the overall accuracy varies in accordance with Jacob’s ladder. This investigation demonstrates the need for systematic, stepwise benchmark calculations for different scales and types of interactions in these complex ligand-protected cluster systems. At this stage, we recommend PBE-D3 (zero) for the structural optimization of phosphine gold clusters and B2PLYP-D3 or PBE0-D3 for calculating the energetics. The BE is best described by B2PLYP-D3, where the necessary MP2 calculations are possible. However, PBE0-D3 must be the method of choice for calculating the BE for larger systems, in general in combination with BJ-damping for the dispersion correction.

Acknowledgments DM and NG thank the University of Canterbury’s BlueFern supercomputer facility. The authors wish to acknowledge the contribution of NeSI high-performance computing facilities to the results of this research. NZ’s national facilities are provided by the NZ eScience Infrastructure. The authors gratefully acknowledge Dr. Vladimir Golovko for fruitful discussions. Keywords: benchmarking
cluster
DFT
catalysis
gold WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

How to cite this article: D. Mollenhauer, N. Gaston. J. Comput. Chem. 2014, 35, 986–997. DOI: 10.1002/jcc.23578

]

Additional Supporting Information may be found in the online version of this article.

[1] [2] [3] [4]

W. J. Jahn, Struct. Biol. 1999, 127, 106. M.-C. Daniel, D. Astruc, Chem. Rev. 2004, 104, 293. Y. Zhu, H. Qian, R. Jin, J. Mat. Chem. 2011, 21, 6793. R. Meyer, C. Lemire, S. K. Shaikhutdinov, H.-J. Freund, Gold Bull. 2004, 37, 72. B. Hvolbak, T. V. W. Janssens, B. S. Clausen, H. Falsig, C. H. Christensen, J. K. Norskov, Nanotoday 2007, 2, 14. C. M. Aikens, Mol. Simul. 2012, 38, 607. S. A. Ivanov, I. Arachchige, C. M. Aikens, J. Phys. Chem. A 2011, 115, 8017. E. B. Guidez, A. Hadley, C. M. Aikens, J. Phys. Chem. C 2011, 115, 6305. S. Goel, K. A. Velizhanin, A. Piryatinski, S. Tretiak, S. A. Ivanov, J. Phys. Chem. Lett. 2010, 1, 927. J. Kilmartin, R. Sarip, R. Grau-Crespo, D. D. Tommaso, G. Hogarth, C. Prestipipo, G. Sankar, ACS Catal. 2012, 2, 957. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Sch€ utz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. K€ oppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklas€s, D. P. O’Neill, P. Palmieri, D. Peng, K. Pfl€ uger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, M. Wang, MOLPRO. TURBOMOLE V6.4 2012, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, N. J. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. van Dam, D. Wang, J. Nieplocha, E. Apra, T. Windus, W. de Jong, Comput. Phys. Commun. 2010, 181, 1477. T. H. J. Dunning, J. Chem. Phys. 1989, 90, 1007. R. A. Kendall, T. H. J. Dunning, R. J. Harrison, J. Chem. Phys. 1992, 96, 6796. D. Figgen, G. Rauhut, M. Dolg, H. Stoll, Chem. Phys. 2005, 311, 227. K. Peterson, C. Puzzarini, Theor. Chem. Acc. 2005, 114, 283. C. Hampel, K. Peterson, H.-J. Werner, Chem. Phys. Lett. 1992, 190, 1. P. J. Knowles, C. Hampel, H.-J. Werner, J. Chem. Phys. 1993, 99, 5219.

[5] [6] [7] [8] [9] [10] [11]

[12]

[13]

[14]

[15] [16] [17] [18] [19] [20]

[21] P. J. Knowles, C. Hampel, H.-J. Werner, J. Chem. Phys. 2000, 112, 3106. [22] M. J. O. Deegan, P. J. Knowles, Chem. Phys. Lett. 1994, 227, 321. [23] R. D. Amos, J. S. Andrews, N. C. Handy, P. J. Knowles, Chem. Phys. Lett. 1991, 185, 256. [24] S. Grimme, J. Chem. Phys. 2003, 118, 9095. [25] R. Lindh, U. Ryu, B. Liu, J. Chem. Phys. 1991, 95, 5889. [26] S. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. [27] J. P. Perdew, Y. Wang, Phys. Rev. B 1992, 45, 13244. [28] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. [29] Y. Zhang, W. Yang, Phys. Rev. Lett. 1998, 80, 890. [30] A. D. Becke, Phys. Rev. A 1988, 38, 3098. [31] J. P. Perdew, Phys. Rev. B 1986, 33, 8822. [32] J. M. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Phys. Rev. Lett. 2003, 91, 146401. [33] Y. Zhao, D. G. Truhlar, J. Chem. Phys. 2006, 125, 194101. [34] A. D. Becke, J. Chem. Phys. 1993, 98, 5648. [35] C. Lee, W. Yang, R. Parr, Phys. Rev. B 1993, 37, 785. [36] J. P. Perdew, M. Ernzerhof, K. Burke, J. Chem. Phys. 1996, 105, 9982. [37] Y. Zhao, D. G. Truhlar, Theor. Chem. Acc. 2008, 120, 215. [38] Y. Zhao, D. G. Truhlar, J. Phys. Chem. 2006, 110, 5121. [39] Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2006, 110, 13126. [40] Y. Minenkov, G. Occhipinti, V. R. Jensen, J. Phys. Chem. A 2009, 113, 11833. [41] Y. Zhao, N. E. Schultz, D. G. Truhlar, J. Chem. Phys. 2005, 123, 161103. [42] Y. Zhao, N. E. Schultz, D. G. Truhlar, J. Chem. Theory Comput. 2006, 2, 364. [43] T. Yanai, D. Tew, N. Handy, Chem. Phys. Lett. 2004, 393, 51. [44] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 2001, 115, 3540. [45] S. Grimme, J. Chem. Phys. 2006, 124, 034108. [46] D. Mollenhauer, J. Floß, H.-U. Reissig, E. Voloshina, B. Paulus, J. Comput. Chem. 2011, 32, 1839. [47] D. R. Lide, CRC Handbook of Chemistry and Physics, 87th ed.; CRC Press: Boca Raton, 2006. [48] T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic Structure Theory; Wiley: Chichester, 2000. [49] S. Grimme, J. Antony, S. Ehrlich, H. J. Krieg, Chem. Phys. 2010, 132, 154104. [50] S. Grimme, S. Ehrlich, L. Goerigk, J. Comput. Chem. 2011, 32, 1456. [51] M. Sierka, A. Hogekamp, R. Ahlrichs, J. Chem. Phys. 2003, 118, 9136. [52] K. Eichkorn, O. Treutler, H. Ohm, M. Haser, R. Ahlrichs, Chem. Phys. Lett. 1995, 240, 283. [53] K. Eichkorn, F. Weigend, O. Treutler, R. Ahlrichs, Theor. Chem. Acc. 1997, 97, 119. [54] R. S. Mulliken, J. Chem. Phys. 1955, 23, 1833. [55] E. D. Glendening, A. E. Reed, J. E. Carpenter, F. Weinhold, NBO Version 3.1. [56] R. F. W. Bader, Atoms in Molecules—A Quantum Theory; Oxford University Press, 1990. [57] W. Tang, E. Sanville, G. Henkelman, J. Phys.: Condens. Matter 2009, 21, 084204. [58] J. Granatier, M. Urban, A. Sadlej, Chem. Phys. Lett. 2010, 484, 154. [59] C. E. Briant, K. P. Hall, D. P. Mingos, J. Organomet. Chem. 1983, 254, 18.

Received: 29 November 2013 Revised: 2 February 2014 Accepted: 7 February 2014 Published online on 20 March 2014

Journal of Computational Chemistry 2014, 35, 986–997

997