Benchmark Study of Density Functional Theory for Neutral Gold

Article pubs.acs.org/JPCA

Benchmark Study of Density Functional Theory for Neutral Gold Clusters, Aun (n = 2−8) Heehyun Baek, Jiwon Moon, and Joonghan Kim* Department of Chemistry, The Catholic University of Korea, Bucheon 14662, Republic of Korea S Supporting Information *

ABSTRACT: Neutral gold clusters, Aun (n = 2−8), were optimized using coupled cluster singles and doubles with perturbative triples (CCSD(T)) with a triple-ζ-level basis set to develop reliable reference values for their structural and energy parameters in order to assess the performance of density functionals. The performance of 44 density functional theory (DFT) methods for calculating molecular structures and relative energies is assessed with respect to CCSD(T). In addition, their performance when calculating vertical ionization potentials (vIPs) of Aun (n = 2−8) is also assessed by comparison with experimental values. The revTPSS functional shows good performance for calculating both the structural and energy properties of Aun (n = 2−8), whereas B3P86 shows a remarkable performance in calculating the vIPs. The quadruple-ζ-level valence basis set is necessary for obtaining accurate energy values in CCSD(T) calculations.

1. INTRODUCTION

Recently, thiolated gold clusters have attracted much attention because of their fascinating catalytic and optical properties.16−20 The only practical theoretical tools for investigating such materials are DFT methods. However, systematic benchmark DFT studies of the molecular properties of neutral gold clusters are still lacking. Therefore, the abovementioned high-level ab initio calculations using CCSD(T) with reasonably large basis sets are required to provide valuable references for assessing DFT exchange-correlation functionals to find the optimal DFT functionals. Additionally, although the vertical ionization potentials (vIPs) of neutral gold clusters have been well characterized experimentally,21 systematic benchmark studies of DFT methods for calculating vIPs of neutral gold clusters are scarce. In this work, geometry optimizations were performed for Aun (n = 2−8) using CCSD(T) with the triple-ζ level basis set. Moreover, the relative energies of Aun (n = 2−8) were also calculated using CCSD(T) with the quadruple-ζ-level basis set. Reliable geometrical and thermochemical references for Aun (n = 2−8) have thereby been developed. To the best of our knowledge, this is the first time the structures of Aun (n = 2−8) have been fully optimized using CCSD(T) with the triple-ζlevel basis set. On the basis of the developed references, assessments of DFT functionals for the prediction of the molecular structures and relative energies of Aun (n = 2−8) have been performed in order to find the optimal DFT

Gold clusters have attracted considerable attention for their selective catalytic reactions.1 Surface roughening may play an important role in determining the catalytic reactivity of gold clusters because the nonplanarity of gold clusters localizes the electron density and promotes chemical reactivity.2 As such, the planarity or nonplanarity of gold clusters has been extensively investigated.3−11 The molecular structure of the ground state of Au8 was a long-standing issue from the first coupled-cluster singles and doubles with perturbative triples (CCSD(T)) single-point calculations, which gave results that contradicted those from density functional theory (DFT)12,13 calculations; the planar structure was predicted as the ground state of Au8 by DFT.9 However, this controversy has recently been resolved because experimental evidence has shown that the planar structure is a minimum.14 The reason this issue was controversial for such a long time was a lack of geometry optimization using high-level ab initio methods such as CCSD(T) with reasonably large basis sets; most CCSD(T) calculations were performed using singlepoint energy calculations on molecular structures that had been optimized by Møller−Plesset second-order perturbation theory (MP2) or DFT.5,8,10,15 Geometry optimizations using CCSD(T) were performed with a double-ζ basis set.5 In addition, most CCSD(T) calculations excluded the 5s and 5p electrons from the valence shells,5,8,9 which reduced the reliability of the calculations. Therefore, geometry optimizations of neutral gold clusters using CCSD(T) with larger basis sets that include the 5s and 5p electrons in the valence shells are necessary. © XXXX American Chemical Society

Received: November 25, 2016 Revised: March 6, 2017 Published: March 14, 2017 A

DOI: 10.1021/acs.jpca.6b11868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. List of the Exchange-Correlation Density Functionals Used in This Work GGA functional BLYP BP86 BPW91 mPWPW91 PBE SOGGA11 N12

meta-GGA refs 34, 34, 55, 31, 54 66 63

52 53 56 59

global-hybrid GGA

global-hybrid meta-GGA

range-separated hybrid

functional

refs

functional

refs

functional

refs

functional

refs

τ-HCTH TPSS M06-L revTPSS M11-L MN12-L

38 69 77 57 60 64

B3LYP B3P86 B3PW91 mPW1PW91 PBEh1PBE B97-1 B97-2 B98 PBE0 BHandHLYP X3LYP SOGGA11-X APF APFD

35, 52 35, 53 35 30 41 42 73 67 32 36 74 61 33 33

B1B95 τ-HCTHhyb TPSSh BMK M05 M06

37 38 68 39 76 78

OHSE1PBE HSE03 HSE06 CAM-B3LYP LC-ωPBE ωB97X ωB97XD HISSb M11 N12-SX MN12-SX

43, 45−47, 49−51 43, 45−47, 49−51 43, 45−47, 49−51 75 70−72 40 40 44 62 65 65

calculations with those basis sets. The calculated results are summarized in Figure S1, and the discussion is in the Supporting Information. In addition, T1 diagnostics of CCSD(T)/dhf-TZVPP of Aun (n = 2−8) were calculated to examine the reliability of CCSD(T), and the results are summarized in Table S1 in the Supporting Information. As can be seen from Table S1 in the Supporting Information, the values of the T1 diagnostic of almost all structures of Aun (n = 2−8) are less than 0.02, indicating that the reference wave functions (RHF and ROHF) of CCSD(T) are reasonably good and the results of CCSD(T) are reliable.29 Geometry optimizations and subsequent harmonic vibrational frequency calculations were performed using 44 DFT methods.30−78 All DFT methods used in this work are summarized in Table 1. In all DFT calculations, dhf-TZVPP was used for the geometry optimizations and subsequent harmonic vibrational frequency calculations. For single-point energy calculations, dhf-QZVPP was used, as in the CCSD(T) calculations. In addition, the vIPs of neutral Aun (n = 2−8) were calculated using single-point calculations with dhf-QZVPP on the structures optimized by dhf-TZVPP. The zero-point energies (ZPEs) of neutral Aun (n = 2−8) calculated by dhfTZVPP were used. An UltraFine integration grid option with 99 radial shells and 590 angular points per shell was used for all DFT calculations in Gaussian 09 to ensure accuracy. All DFT calculations were performed using the Gaussian 09 program.28 The spin−orbit (SO) effect on the vIPs was examined using the two-component spin−orbit DFT (SODFT) method.79 In the SODFT calculations, only the B3P86 functional was used on the structure optimized by B3P86/dhf-TZVPP. The [2p1d] functions were augmented into the valence basis functions of dhf-QZVPP in the SODFT calculations.24 All SODFT calculations were performed using the NWChem6.6 program.80

functional. In addition, the performance of DFT for calculating vIPs was also assessed by comparison with experimental values.21 This work will provide a reference for the further development of new DFT functionals and guidelines for selecting DFT functionals for calculating the molecular properties of large Au-containing clusters.

2. COMPUTATIONAL DETAILS CCSD(T)22 based on the restricted Hartree−Fock (RHF) reference wave function (singlet) and unrestricted CCSD(T) (UCCSD(T))23 based on the restricted open-shell HF (ROHF) reference wave function (doublet) were used for even-numbered and odd-numbered neutral gold clusters, respectively. Hereafter, we denote both methods as CCSD(T). Geometry optimizations of Aun (n = 2−8) were performed using the CCSD(T) method. The dhf-TZVPP24 relativistic effective core potential (RECP) was used for Au atoms in order to consider the scalar relativistic effect. It is noted that the pseudopotential of dhf-TZVPP is the same as that of the augcc-pVTZ-pp (ECP60MDF) RECP;25 60 electrons of Au were treated as a core, small-core RECP. The basis set of the valence electrons is a segmented contract triple-ζ level, [6s5p3d2f1g]. The 5s, 5p, 5d, and 6s electrons of Au were correlated in all CCSD(T) calculations; in other words, the core−valence correlation effect was considered in this work. It is noted that the core sizes of the CCSD(T) calculations in this work are smaller than those in previous theoretical investigations,5,8,9 indicating that the CCSD(T) calculations in this work are quite reliable. Subsequent single-point energy calculations using CCSD(T) with quadruple-ζ basis sets (dhf-QZVPP, [7s6p4d4f2g])24 were performed on the optimized structures by CCSD(T)/dhf-TZVPP (except for structures 7.a and 7.c, see the main text) to obtain more accurate energetics. No vibrational frequency calculations were carried out for any of the CCSD(T) calculations because of computational cost. All CCSD(T)22 calculations with dhf-TZVPP were performed using the Molpro2012 program,26 and all CCSD(T)27 singlepoint calculations with dhf-QZVPP were performed using the Gaussian 09 program.28 The reference wave function of UCCSD(T)/dhf-QZVPP for odd-numbered gold clusters when using Gaussian 09 is an unrestricted HF (UHF) wave function. To examine the performance of dhf-T(Q)ZVPP compared to aug-cc-pVT(Q)Z-pp for the structural and energy parameters of Au3 and Au4, we performed the CCSD(T)

3. RESULTS AND DISCUSSION 3.1. Results of Aun (n = 2−8) using CCSD(T). Geometry optimizations of Aun (n = 2−8) (except for structures 7.a and 7.c) were performed using CCSD(T)/dhf-TZVPP, and the molecular structures of Aun (n = 2−8) are shown in Figure 1. (xyz coordinates of all structures are listed in the Supporting Information (SI).) In this section, the results calculated in this work using CCSD(T) are compared to the results of high-level ab initio calculations and experimental findings (for Au2 only).81 High-level ab initio calculations for odd-numbered B


Article

The Journal of Physical Chemistry A

Figure 1. Optimized molecular structures of Aun (n = 2−8) using CCSD(T)/dhf-TZVPP, with the exception of structures 7.a and 7.c, which were optimized using revTPSS/dhf-TZVPP. (See the main text.) Values in parentheses are the relative electronic energies (in kcal/mol) using CCSD(T)/ dhf-TZVPP and CCSD(T)/dhf-QZVPP//CCSD(T)/dhf-TZVPP [CCSD(T)/dhf-QZVPP//revTPSS/dhf-TZVPP for 7.c] in italics.

gold clusters, in particular, Au5 and Au7, are scarce because of the computational cost of the use of unrestricted formalism. Thus, the discussion will focus on the even-numbered neutral gold clusters. The Au−Au bond length of Au2 (2.490 Å) is reasonably close to the experimental value (2.472 Å).81 In addition, the optimized bond length of Au2 in this work is in better agreement with the experimental value than that

obtained using the other triple-ζ basis set (CCSD(T)/ccpVTZ-PP, 2.497 Å).15 The ground state of Au3 is a distorted triangular structure (3.a) due to the Jahn−Teller distortion. The geometry optimization of the ground state (2B2) of Au3 was performed using CCSD(T) with the large-core RECP in conjunction with the [8s4p5d3f] valence basis set in a previous theoretical study.82 The optimized bond length of structure 3.a C


Article


Figure 2. MAEs and maximum deviations of the calculated (a) Au−Au bond lengths (47 bonds total) and (b) ∠Au−Au−Au bond angles (30 angles and 1 dihedral angle total) with respect to the values of CCSD(T)/dhf-TZVPP.

(2.607 Å) is larger than our optimized value (2.581 Å), which can be ascribed to the use of large-core RECP with a smaller valence basis set. Despite the overestimation of the bond length, the bond angle (65.4°) is in good agreement with ours (65.2°). The relative energy between structures 3.a and 3.b was estimated using CCSD(T) with small-core RECP with the [5s5p4d3f2g] valence basis set including SO effects.83 The estimated relative energy is about 0.17 eV (3.9 kcal/mol), which is in good agreement with our calculated value (2.1 kcal/ mol). The relative energies of three isomers of Au4 (4.a, 0.0; 4.b, 2.9; and 4.c, 17.3 kcal/mol) have been investigated at the CCSD(T)/aug-cc-pVTZ-PP//SO-TPSSh/aug-cc-pVTZ-PP level,84 and the findings are in excellent agreement with our calculated results (4.a, 0.0; 4.b, 2.8; and 4.c, 16.5 kcal/mol). Olson et al.9 investigated isomers of Au6 using a single-point CCSD(T) calculation via SBKJC RECP augmented with one

set of f functions on the structures optimized by MP2 with the same basis set. Because of the poor MP2 structure and low-level basis set, the relative energies of Au6 (6.a, 0.0 and 6.b, 15.0 kcal/mol) are significantly underestimated when compared to our calculated values (6.a, 0.0 and 6.b, 22.2 kcal/mol). However, recent CCSD(T) calculations with an improved basis set gave relative energies of Au6 that are quite close to ours (6.a, 0.0 and 6.b, 22.5 kcal/mol).4 The planar or nonplanar structure of the ground state of Au8 is a long-standing issue in the field of computational chemistry, but recent experimental evidence has indicated that the planar structure (8.a) is the ground state.14 It is noted that this work contains the first full geometry optimizations of isomers of Au8 using CCSD(T) with a triple-ζ basis set. The CCSD(T)/dhf-TZVPP calculations give a planar structure, 8.a, as the ground state, which is in line with experimental observations. As can be seen from Figure 1, D


Article


Figure 3. MAEs and maximum deviations of the relative energies (in kcal/mol) with respect to the values of CCSD(T)/dhf-QZVPP//CCSD(T)/ dhf-TZVPP [CCSD(T)/dhf-QZVPP//revTPSS/dhf-TZVPP for the 7.c structure, see the main text].

CCSD(T). MN12-L (0.009/0.032 Å) also shows comparable performance to revTPSS. It is noted that these two DFT functionals are local functionals. Both DFT functionals (revTPSS and MN12-L) show better performance for predicting the bond lengths of Aun (n = 2−8) than their counterpart hybrid functionals (TPSSh and MN12-SX). As shown in Figure 2a, TPSS (meta-GGA, 0.022/0.036 Å) and TPSSh (hybrid-meta, 0.020/0.033 Å) exhibit reasonable performance for optimizing bond lengths. These results indicate that the TPSS family (TPSS, revTPSS, and TPSSh) generally work well for predicting the molecular structures of Aun (n = 2−8). As can be seen from Figure 2a, M05 (0.140/ 0.211 Å) has the worst performance for predicting the bond lengths of Aun (n = 2−8), significantly overestimating the Au− Au bond lengths. M06 (0.104/0.144 Å) follows M05, having the second worst performance. These results indicate that typical Minnesota functionals (M05 and M06) are not adequate functionals for predicting the molecular structures of Aun (n = 2−8). Recent hybrid Minnesota functionals such as SOGGA11X, M11, N12-SX, and MN12-SX improve the performance, but they are still worse than revTPSS and MN12-L for calculating the molecular structures of Aun (n = 2−8). BLYP (0.089/0.201 Å) also gives poor results, just following the M05 functionals. B3LYP (0.072/0.118 Å), a hybrid version of BLYP, also gives poor performance, although their results are better than those of BLYP. One should be careful to use B3LYP, the most popular DFT functional, for optimizing the bond lengths of Aun (n = 2−8). According to the MAEs, most DFT methods work well for predicting the bond angles of Aun (n = 2−8), with MAE values being less than 1.0°. The hybrid revised form of the pure PBE functional, PBEh1PBE (0.5/2.1°) shows the best performance for calculating the bond angles of Aun (n = 2−8). In addition, almost all hybrid versions of PBE functionals (PBE0, OHSE1PBE, HSE03, and HSE06) provide reasonable performance. The TPSS family of functionals (TPSSh (0.5/3.5°), TPSS (0.6/6.7°), and revTPSS (0.8/6.8°)) also exhibits good performance. In contrast, M11 (6.3/72.7°) has the worst performance for calculating the bond angles. BLYP and τ-

using a quadruple-ζ basis set leads to the relative energies being closer together and inverts the energy ordering between structures 8.b and 8.c compared to that found using CCSD(T)/ dhf-TZVPP. In addition, the order of the relative energies of 8.b (4.60 kcal/mol) and 8.c (4.39 kcal/mol) calculated by CCSD(T)/cc-pVTZ-PP//CCSD(T)/cc-pVDZ-PP in recent theoretical work5 differs from that in this work (8.b, 4.5 and 8.c, 4.8 kcal/mol). These results indicate that a quadruple-ζ basis set is necessary to obtain accurate energetics for neutral gold clusters. The relative energies of Au8 as calculated by CCSD(T)/dhf-QZVPP//CCSD(T)/dhf-TZVPP (8.a, 0.0; 8.b, 4.5; and 8.d, 6.4 kcal/mol) are quite close to those found using CCSD(T)/cc-pVQZ-PP (8.a, 0.0; 8.b, 4.3; and 8.d, 6.3 kcal/ mol) and CCSD(T) with a complete basis set (CBS) limit (8.a, 0.0; 8.b, 3.8; and 8.d, 6.3 kcal/mol).15 These results indicate that the level of theory used in this work is reliable. 3.2. Performance of DFT for Determining Molecular Structures of Aun (n = 2−8). We considered 47 bond lengths and 30 bond angles when assessing the performance of the DFT functionals for predicting the geometrical parameters of Aun (n = 2−8) in this study. One dihedral angle has also been considered to clarify the dimensionality issue in structure 6.b. The geometrical parameters of 7.a and 7.c have been excluded from the assessments because the optimized molecular structures of 7.a and 7.c using CCSD(T) were not available. The mean absolute errors (MAEs) of the Au−Au bond lengths and ∠Au−Au−Au bond angles with respect to the CCSD(T)/ dhf-TZVPP reference values are shown in Figure 2a,b, respectively. (All MAEs and maximum deviations of geometrical parameters are summarized in Tables S2 and S3 in the SI.) In Figure 2a,b, the two vertical axes on the left and right have different scales, showing MAE and the maximum deviation, respectively. As can be seen in Figure 2a, revTPSS (meta-generalized gradient approximation (GGA)) exhibits the best performance for predicting the bond length (MAE: 0.008/maximum deviation: 0.042 Å, hereafter we use MAE/maximum deviation notation), indicating that the bond lengths of Aun (n = 2−8) as optimized by revTPSS are quite close to those obtained by E


Article


Figure 4. MAEs and maximum deviations of the calculated vIPs (in eV) with respect to experimental values (taken from ref 21).

3.3. Performance of DFT for Calculating the Relative Energies of Aun (n = 2−8). Figure 3 shows the MAEs and maximum deviations of the relative energies of Aun (n = 2−8) with respect to those of CCSD(T)/dhf-QZVPP//CCSD(T)/ dhf-TZVPP, with the exception of structures 7.a and 7.c. All MAEs and maximum deviations of relative energies are summarized in Table S4 in the SI. For structure 7.c, the single-point energy calculation using CCSD(T)/dhf-QZVPP was performed on the structure optimized by revTPSS/dhfTZVPP because the revTPSS functional gives molecular structures of Aun (n = 2−8) that are closest to those of CCSD(T) according to the results in the above section. However, for structure 7.a, which has Cs symmetry, the CCSD(T)/dhf-QZVPP single-point energy calculation could not performed because the computational cost is too high to be feasible using our computer facilities. Thus, only the relative energies of structures 7.b and 7.c are included for the assessment of the energy parameters of Au7. As shown in Figure 3, revTPSS, the best functional for calculating the bond lengths, also shows the best performance for calculating the relative energies (1.6/4.2 kcal/mol) of Aun (n = 2−8). B1B95 (2.4/4.6 kcal/mol) and TPSSh (2.6/5.2 kcal/mol) follow it. TPSS also exhibits good performance (3.0/ 6.0 kcal/mol); therefore, the TPSS family (TPSS, revTPSS, and TPSSh) perform well in calculating the relative energies as well as for predicting the molecular structures of Aun (n = 2−8). Although MN12-L gives very good results for predicting the molecular structures of Aun (n = 2−8) that are comparable to those obtained using revTPSS, it shows only modest results when calculating the relative energies. M11, the worst functional for calculating bond angles, also shows the worst performance for energy calculations (11.4/26.8 kcal/mol). The M11-L local functional also shows a large MAE value. It is noted that M11 and M11-L exhibit poor performance for calculating the energetics of Ru−Ru bonds contained in organometallic complexes in a recent study.85 Thus, the M11 family of functionals is not recommended for calculating energetics related to the bonds between transition metals. τHCTH (10.1/19.6 kcal/mol), BLYP (9.9/19.1 kcal/mol), and

HCTH also give large MAEs. As shown in Figure 2b, six DFT functionals in particular (BHandHLYP, CAM-B3LYP, LCωPBE, ωB97X, SOGGA11, and M11) give distinctly larger MAEs and maximum deviations when compared to other functionals, attributed to their incorrect prediction that the structure of 6.b would be planar. Thus, one should be careful when using these six DFT functionals to predict the turning points from 2D to 3D structures for neutral gold clusters. In addition, as can be seen from Figure 2a, these six DFT functionals also give large MAE values when predicting the Au−Au bond lengths. Therefore, they are strongly not recommended for calculating the molecular structures of neutral gold clusters. The importance of considering the empirical dispersion effect when predicting the molecular structures of neutral gold clusters can be identified by comparing the MAE values of ωB97X and APF with those of ωB97XD and APFD. As can be seen from Figure 2a, a consideration of the dispersion interaction slightly improves the prediction of the Au−Au bond length. The same results can be observed in the bond angles of ωB97X and ωB97XD. However, contrasting results are observed in the MAEs of the bond angles of APF and APFD, where considering the dispersion interaction leads to further deviation from the reference values. These results indicate that the dispersion interaction has a minor effect on the molecular structures of neutral gold clusters. In summary, the revTPSS functional works quite well for predicting the molecular structures of Aun (n = 2−8), so it is recommended for predicting reasonable molecular structures for neutral gold clusters. Therefore, the revTPSS functional could be useful for predicting the molecular structures of large neutral gold clusters. The good performance of revTPSS for predicting the preferences of neutral gold clusters for 2D or 3D molecular structures has been reported in a recent study.6 In that study, the accuracy of revTPSS was established by comparing the results of revTPSS with those of the random phase approximation (RPA). The calculated results in this work are in line with the previous study.6 F


Article


performance (0.47/0.84 eV) than almost all of the other DFT functionals. When comparing B3P86 and BP86, mixing exact exchange seems to improve the performance for calculating the vIPs. However, as can be seen from Figure 4, the MAE of another GGA and its counterpart hybrid functionals show that mixing exact exchange cannot guarantee better performance. For example, the MAE of PBE0 (hybrid) is larger than that of PBE (GGA). These results indicate that the B exchange and P86 correlation functional are themselves optimal for calculating the vIPs. M11-L (1.34/1.80 eV) shows the worst performance for calculating vIPs. Subsequently, M05 (1.04/ 1.26 eV), M06-L (0.99/1.17 eV), BHandHLYP (0.97/1.31 eV), and LC-ωPBE (0.95/1.33 eV) give large MAEs and maximum deviations. Because of the large spin−orbit coupling (SOC) of Au, SOC may decrease the MAE of vIPs, leading to better agreement with the experimental values. To examine this possibility, SODFT calculations with B3P86 (SO-B3P86) were performed to calculate the vIPs of Aun (n = 2−4). The calculated results for the vIPs of Aun (n = 2−4) using B3P86 and SO-B3P86, are summarized in Table S7 in the SI, in which the SO effect on the vIPs is negligible and the difference between B3P86 and SOB3P86 is less than 0.08 eV. These results can be attributed to the cancelation of SOC on both Aun and Aun+ (n = 2−4). According to these results, the source of the poor performance when calculating the vIPs of Aun (n = 2−8) is not the SO effect but the DFT functional itself. Thus, an improvement in DFT functionals will be necessary to predict the vIPs of Aun correctly, with these being closely related to their catalytic properties.

SOGGA11 (9.6/34.3 kcal/mol), which exhibit poor performance for calculating bond angles, give large MAEs and maximum deviations for energy calculations as well. B3LYP also shows poor performance (8.8/19.2 kcal/mol). As mentioned in the above section, B3LYP cannot provide accurate structural parameters; therefore, it is not recommended for calculating the molecular properties of Aun (n = 2− 8). Functionals that contain LYP correlation functionals (BLYP, B3LYP, BHandHLYP, X3LYP, and CAM-B3LYP) generally have large MAE values. As shown in Figure 3, hybrid functionals generally show slightly better performance than their counterpart GGA functionals, with the exception of recent Minnesota functionals such as SOGGA11-X, N12-SX, and MN12-SX. As can be seen by comparing the results of ωB97X and APF with those of ωB97XD and APFD, the inclusion of empirical dispersion results in a slight increase in MAE. Therefore, the dispersion interaction does not have a major effect on the relative energies in neutral gold clusters, as in the molecular structures. In addition, B1B95, M06-L, APFD, M11L, and MN12-L did not give structure 8.a (D4h) as the ground state of Au8. Although B1B95 works well for calculating the relative energies of Aun (n = 2−8), it cannot be used to predict planar−nonplanar turning points. The revTPSS functional provides sufficiently good performance for calculating both the structural and energy properties of Aun (n = 2−8). These results are positive because revTPSS is a local functional, and thus the computational cost of using revTPSS is cheaper than that of using hybrid functionals; therefore, the application of revTPSS to large neutral gold clusters is feasible. In addition, revTPSS exhibited good performance for predicting the 2D−3D crossover points of cationic (Au8+) and anionic (Au12−) gold clusters.6 It could be expected that revTPSS may also work well for predicting the properties of charged gold clusters. Thus, the use of the revTPSS functional with a reasonable-sized basis set is strongly recommended for calculating the properties of large neutral and charged gold clusters. 3.4. Performance of DFT for Calculating vIPs of Aun (n = 2−8). Figure 4 shows the MAEs and maximum deviations of the calculated vIPs compared to experimental values.21 All calculated vIPs, experimental values, MAEs, and maximum deviations of Aun (n = 2−8) are summarized in Table S5 in the SI. We also calculated vIPs of Aun (n = 2−5) using CCSD(T)/ dhf-QZVPP//CCSD(T)/dhf-TZVPP, and the results are shown in Table S6 in the SI. Because of the computational cost of CCSD(T), vIPs of up to Au5 are available in the CCSD(T) calculations. As shown in Table S6 in the SI, the results of CCSD(T) are modest. This result prompts us to find an alternative DFT functional that exhibits a better performance for calculating vIPs than CCSD(T). The first thing to note is that almost none of the DFTs work well for calculating the vIPs of Aun (n = 2−8), with the MAEs of almost all of the DFTs reaching ∼0.8 eV. Although revTPSS gives satisfactory results for the molecular structures and relative energies of Aun (n = 2−8), it shows poor performance for calculating the vIPs. Even the recently developed Minnesota functionals (SOGGA11, M11-L, N12, MN12-L, SOGGA11X, M11, N12SX, and MN12-SX) exhibit poor performance. However, as shown in Figure 4, the B3P86 functional shows remarkably good performance in the calculations of vIPs (0.18/0.38 eV). HSE03 (0.47/0.72 eV) and BP86 (0.47/0.84 eV) follow it, but a large gap exists between the MAE of B3P86 and theirs. It is noted that a rather old GGA functional, BP86, shows better

4. CONCLUSIONS In this work, reliable geometrical and energy references have been developed for Aun (n = 2−8) from optimized structures using CCSD(T) with the triple-ζ-level basis set and relative energies using CCSD(T) with the quadruple-ζ-level basis set. In addition, the performance of various DFT functionals when calculating molecular structures and relative energies was assessed with respect to CCSD(T). The performances when calculating vIPs of Aun (n = 2−8) was also assessed by comparing the results with experimental values. The local functional revTPSS provides sufficiently good performance for calculating both the structural and energy properties of Aun (n = 2−8). Therefore, revTPSS with a reasonably large basis set is strongly recommended for calculating the structural and energy properties of large gold-containing clusters. B3P86 demonstrates superb performance for calculating the vIPs. In addition, using a quadruple-ζ basis set leads to the relative energies being closer and the energy order being inverted at Au8, and as such, quadruple-ζ basis sets are necessary to obtain accurate energetics of neutral gold clusters. On the other hand, it has been demonstrated that SOC has a minor effect when calculating the vIPs of neutral gold clusters.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b11868. Performance of dhf-T(Q)ZVPP compared to aug-ccpVT(Q)Z-pp for the structural and energy parameters of Au3 and Au4, values of the T1 diagnostic of CCSD(T)/ dhf-TZVPP of Aun, MAEs and maximum deviations of G


Article


■

(14) Gruene, P.; Butschke, B.; Lyon, J. T.; Rayner, D. M.; Fielicke, A. Far-IR spectra of small neutral gold clusters in the gas phase. Z. Phys. Chem. 2014, 228, 337−350. (15) Diefenbach, M.; Kim, K. S. Spatial Structure of Au8: importance of basis set completeness and geometry relaxation. J. Phys. Chem. B 2006, 110, 21639−21642. (16) Abbas, M. A.; Kim, T. Y.; Lee, S. U.; Kang, Y. S.; Bang, J. H. Exploring Interfacial Events in Gold-Nanocluster-Sensitized Solar Cells: Insights into the Effects of the Cluster Size and Electrolyte on Solar Cell Performance. J. Am. Chem. Soc. 2016, 138, 390−401. (17) Chen, Y.; Zeng, C.; Kauffman, D. R.; Jin, R. Tuning the Magic Size of Atomically Precise Gold Nanoclusters via Isomeric Methylbenzenethiols. Nano Lett. 2015, 15, 3603−3609. (18) Jiang, D. E. The expanding universe of thiolated gold nanoclusters and beyond. Nanoscale 2013, 5, 7149−7160. (19) Knoppe, S.; Boudon, J.; Dolamic, I.; Dass, A.; Bürgi, T. Size Exclusion Chromatography for Semipreparative Scale Separation of Au38(SR)24 and Au40(SR)24 and Larger Clusters. Anal. Chem. 2011, 83, 5056−5061. (20) Malola, S.; Lehtovaara, L.; Knoppe, S.; Hu, K.-J.; Palmer, R. E.; Bürgi, T.; Häkkinen, H. Au40(SR)24 Cluster as a Chiral Dimer of 8Electron Superatoms: Structure and Optical Properties. J. Am. Chem. Soc. 2012, 134, 19560−19563. (21) Jackschath, C.; Rabin, I.; Schulze, W. Electronic Structures and Related Properties. Electron Impact Ionization Potentials of Gold and Silver Clusters Men, n≤ 22. Ber. Bunsenges. Phys. Chem. 1992, 96, 1200−1204. (22) Hampel, C.; Peterson, K. A.; Werner, H. J. A comparison of the efficiency and accuracy of the quadratic configuration interaction (QCISD), coupled cluster (CCSD), and Brueckner coupled cluster (BCCD) methods. Chem. Phys. Lett. 1992, 190, 1−12. (23) Watts, J. D.; Gauss, J.; Bartlett, R. J. Coupled-cluster methods with noniterative triple excitations for restricted open-shell Hartree− Fock and other general single determinant reference functions. Energies and analytical gradients. J. Chem. Phys. 1993, 98, 8718−8733. (24) Weigend, F.; Baldes, A. Segmented contracted basis sets for oneand two-component Dirac−Fock effective core potentials. J. Chem. Phys. 2010, 133, 174102. (25) Figgen, D.; Rauhut, G.; Dolg, M.; Stoll, H. Energy-consistent pseudopotentials for group 11 and 12 atoms: adjustment to multiconfiguration Dirac−Hartree−Fock data. Chem. Phys. 2005, 311, 227− 244. (26) Werner, H.; Knowles, P.; Knizia, G.; Manby, F.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G. MOLPRO, version 2012.1, A Package of Ab Initio Programs; 2012. http://www.molpro.net (accessed June 18, 2014). (27) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479−483. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (29) Lee, T. J.; Taylor, P. R. A diagnostic for determining the quality of single-reference electron correlation methods. Int. J. Quantum Chem. 1989, 36, 199−207. (30) Adamo, C.; Barone, V. Toward reliable adiabatic connection models free from adjustable parameters. Chem. Phys. Lett. 1997, 274, 242−250. (31) Adamo, C.; Barone, V. Exchange functionals with improved long-range behavior and adiabatic connection methods without adjustable parameters: The mPW and mPW1PW models. J. Chem. Phys. 1998, 108, 664−675. (32) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (33) Austin, A.; Petersson, G. A.; Frisch, M. J.; Dobek, F. J.; Scalmani, G.; Throssell, K. A density functional with spherical atom dispersion terms. J. Chem. Theory Comput. 2012, 8, 4989−5007.

Au−Au bond lengths, MAEs and maximum deviations of ∠Au−Au−Au bond angles, MAEs and maximum deviations of the relative energies, MAEs and maximum deviations of vIPs, vIPs of Aun calculated by CCSD(T)/ dhf-QZVPP//CCSD(T)/dhf-TZVPP, vIPs calculated by B3P86 and SO-B3P86, and xyz coordinates of Aun optimized by CCSD(T)/dhf-TZVPP (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +82-2-2164-4338. Fax: +82-2-2164-4764. ORCID

Joonghan Kim: 0000-0002-7783-0200 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2014R1A1A1007188). This work was also supported by the National Institute of Supercomputing and Network/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2016-C1-0002).

■

REFERENCES

(1) Hayashi, T.; Tanaka, K.; Haruta, M. Selective Vapor-Phase Epoxidation of Propylene over Au/TiO2Catalysts in the Presence of Oxygen and Hydrogen. J. Catal. 1998, 178, 566−575. (2) Mills, G.; Gordon, M. S.; Metiu, H. Oxygen adsorption on Au clusters and a rough Au (111) surface: The role of surface flatness, electron confinement, excess electrons, and band gap. J. Chem. Phys. 2003, 118, 4198−4205. (3) Götz, D. A.; Schäfer, R.; Schwerdtfeger, P. The performance of density functional and wavefunction-based methods for 2D and 3D structures of Au10. J. Comput. Chem. 2013, 34, 1975−1981. (4) Han, Y. K. Structure of Au8: Planar or nonplanar? J. Chem. Phys. 2006, 124, 024316. (5) Hansen, J. A.; Piecuch, P.; Levine, B. G. Communication: Determining the lowest-energy isomer of Au8: 2D, or not 2D. J. Chem. Phys. 2013, 139, 091101. (6) Johansson, M. P.; Warnke, I.; Le, A.; Furche, F. At what size do neutral gold clusters turn three-dimensional? J. Phys. Chem. C 2014, 118, 29370−29377. (7) Li, X. B.; Wang, H. Y.; Yang, X. D.; Zhu, Z. H.; Tang, Y. J. Size dependence of the structures and energetic and electronic properties of gold clusters. J. Chem. Phys. 2007, 126, 084505. (8) Olson, R. M.; Gordon, M. S. Isomers of Au8. J. Chem. Phys. 2007, 126, 214310. (9) Olson, R. M.; Varganov, S.; Gordon, M. S.; Metiu, H.; Chretien, S.; Piecuch, P.; Kowalski, K.; Kucharski, S. A.; Musial, M. Where does the planar-to-nonplanar turnover occur in small gold clusters? J. Am. Chem. Soc. 2005, 127, 1049−1052. (10) Serapian, S. A.; Bearpark, M. J.; Bresme, F. The shape of Au8: gold leaf or gold nugget? Nanoscale 2013, 5, 6445−6457. (11) Xiao, L.; Tollberg, B.; Hu, X.; Wang, L. Structural study of gold clusters. J. Chem. Phys. 2006, 124, 114309. (12) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864. (13) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133. H


Article


approximation for exchange and correlation. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 4978−4978. (57) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Workhorse semilocal density functional for condensed matter physics and quantum chemistry. Phys. Rev. Lett. 2009, 103, 026403. (58) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 45, 13244−13249. (59) Perdew, J. P.; Ziesche, P.; Eschrig, H. Electronic Structure of Solids’ 91; Akademie Verlag: Berlin, 1991; Vol. 11. (60) Peverati, R.; Truhlar, D. G. M11-L: A local density functional that provides improved accuracy for electronic structure calculations in chemistry and physics. J. Phys. Chem. Lett. 2012, 3, 117−124. (61) Peverati, R.; Truhlar, D. G. Communication: A global hybrid generalized gradient approximation to the exchange-correlation functional that satisfies the second-order density-gradient constraint and has broad applicability in chemistry. J. Chem. Phys. 2011, 135, 191102. (62) Peverati, R.; Truhlar, D. G. Improving the accuracy of hybrid meta-GGA density functionals by range separation. J. Phys. Chem. Lett. 2011, 2, 2810−2817. (63) Peverati, R.; Truhlar, D. G. Exchange−correlation functional with good accuracy for both structural and energetic properties while depending only on the density and its gradient. J. Chem. Theory Comput. 2012, 8, 2310−2319. (64) Peverati, R.; Truhlar, D. G. An improved and broadly accurate local approximation to the exchange−correlation density functional: The MN12-L functional for electronic structure calculations in chemistry and physics. Phys. Chem. Chem. Phys. 2012, 14, 13171− 13174. (65) Peverati, R.; Truhlar, D. G. Screened-exchange density functionals with broad accuracy for chemistry and solid-state physics. Phys. Chem. Chem. Phys. 2012, 14, 16187−16191. (66) Peverati, R.; Zhao, Y.; Truhlar, D. G. Generalized gradient approximation that recovers the second-order density-gradient expansion with optimized across-the-board performance. J. Phys. Chem. Lett. 2011, 2, 1991−1997. (67) Schmider, H. L.; Becke, A. D. Optimized density functionals from the extended G2 test set. J. Chem. Phys. 1998, 108, 9624−9631. (68) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes. J. Chem. Phys. 2003, 119, 12129−12137. (69) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta−generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 146401. (70) Vydrov, O. A.; Heyd, J.; Krukau, A. V.; Scuseria, G. E. Importance of short-range versus long-range Hartree-Fock exchange for the performance of hybrid density functionals. J. Chem. Phys. 2006, 125, 074106. (71) Vydrov, O. A.; Scuseria, G. E. Assessment of a long-range corrected hybrid functional. J. Chem. Phys. 2006, 125, 234109. (72) Vydrov, O. A.; Scuseria, G. E.; Perdew, J. P. Tests of functionals for systems with fractional electron number. J. Chem. Phys. 2007, 126, 154109. (73) Wilson, P. J.; Bradley, T. J.; Tozer, D. J. Hybrid exchangecorrelation functional determined from thermochemical data and ab initio potentials. J. Chem. Phys. 2001, 115, 9233−9242. (74) Xu, X.; Goddard, W. A. The X3LYP extended density functional for accurate descriptions of nonbond interactions, spin states, and thermochemical properties. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 2673−2677. (75) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange− correlation functional using the Coulomb-attenuating method (CAMB3LYP). Chem. Phys. Lett. 2004, 393, 51−57. (76) Zhao, Y.; Schultz, N. E.; Truhlar, D. Exchange-correlation functional with broad accuracy for metallic and nonmetallic

(34) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (35) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (36) Becke, A. D. A new mixing of Hartree−Fock and local densityfunctional theories. J. Chem. Phys. 1993, 98, 1372−1377. (37) Becke, A. D. Density-functional thermochemistry. IV. A new dynamical correlation functional and implications for exact-exchange mixing. J. Chem. Phys. 1996, 104, 1040−1046. (38) Boese, A. D.; Handy, N. C. New exchange-correlation density functionals: The role of the kinetic-energy density. J. Chem. Phys. 2002, 116, 9559−9569. (39) Boese, A. D.; Martin, J. M. Development of density functionals for thermochemical kinetics. J. Chem. Phys. 2004, 121, 3405−3416. (40) Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom−atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (41) Ernzerhof, M.; Perdew, J. P. Generalized gradient approximation to the angle-and system-averaged exchange hole. J. Chem. Phys. 1998, 109, 3313−3320. (42) Hamprecht, F. A.; Cohen, A. J.; Tozer, D. J.; Handy, N. C. Development and assessment of new exchange-correlation functionals. J. Chem. Phys. 1998, 109, 6264−6271. (43) Henderson, T. M.; Izmaylov, A. F.; Scalmani, G.; Scuseria, G. E. Can short-range hybrids describe long-range-dependent properties? J. Chem. Phys. 2009, 131, 044108. (44) Henderson, T. M.; Izmaylov, A. F.; Scuseria, G. E.; Savin, A. Assessment of a middle-range hybrid functional. J. Chem. Theory Comput. 2008, 4, 1254−1262. (45) Heyd, J.; Peralta, J. E.; Scuseria, G. E.; Martin, R. L. Energy band gaps and lattice parameters evaluated with the Heyd-ScuseriaErnzerhof screened hybrid functional. J. Chem. Phys. 2005, 123, 174101. (46) Heyd, J.; Scuseria, G. E. Efficient hybrid density functional calculations in solids: assessment of the Heyd−Scuseria−Ernzerhof screened Coulomb hybrid functional. J. Chem. Phys. 2004, 121, 1187− 1192. (47) Heyd, J.; Scuseria, G. E. Assessment and validation of a screened Coulomb hybrid density functional. J. Chem. Phys. 2004, 120, 7274− 7280. (48) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207− 8215. (49) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys.118, 8207 (2003)]. J. Chem. Phys. 2006, 124, 219906. (50) Izmaylov, A. F.; Scuseria, G. E.; Frisch, M. J. Efficient evaluation of short-range Hartree-Fock exchange in large molecules and periodic systems. J. Chem. Phys. 2006, 125, 104103. (51) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 2006, 125, 224106. (52) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (53) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (54) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (55) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 6671−6687. (56) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Erratum: Atoms, molecules, solids, and surfaces: Applications of the generalized gradient I


Article

The Journal of Physical Chemistry A compounds, kinetics, and noncovalent interactions. J. Chem. Phys. 2005, 123, 161103. (77) Zhao, Y.; Truhlar, D. G. A new local density functional for maingroup thermochemistry, transition metal bonding, thermochemical kinetics, and noncovalent interactions. J. Chem. Phys. 2006, 125, 194101. (78) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (79) Nichols, P.; Govind, N.; Bylaska, E. J.; De Jong, W. A. Gaussian Basis Set and Planewave Relativistic Spin− Orbit Methods in NWChem. J. Chem. Theory Comput. 2009, 5, 491−499. (80) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. NWChem: A comprehensive and scalable opensource solution for large scale molecular simulations. Comput. Phys. Commun. 2010, 181, 1477−1489. (81) Huber, K. G. Herzberg Constants of Diatomic Molecules; Van Nostrand Reinhold: New York, 1979. (82) Wesendrup, R.; Hunt, T.; Schwerdtfeger, P. Relativistic coupled cluster calculations for neutral and singly charged Au3 clusters. J. Chem. Phys. 2000, 112, 9356−9362. (83) Rusakov, A. A.; Rykova, E.; Scuseria, G. E.; Zaitsevskii, A. Importance of spin-orbit effects on the isomerism profile of Au3: An ab initio study. J. Chem. Phys. 2007, 127, 164322. (84) Shi, Y. K.; Li, Z. H.; Fan, K. N. Validation of Density Functional Methods for the Calculation of Small Gold Clusters. J. Phys. Chem. A 2010, 114, 10297−10308. (85) Paranthaman, S.; Moon, J.; Kim, J.; Kim, D. E.; Kim, T. K. Performance of Density Functional Theory and Relativistic Effective Core Potential for Ru-Based Organometallic Complexes. J. Phys. Chem. A 2016, 120, 2128−2134.

J


Benchmark Study of Density Functional Theory for Neutral Gold

Recommend Documents