Hybrid Density Functionals for Clusters of Late Transition Metals

Article pubs.acs.org/JCTC

Hybrid Density Functionals for Clusters of Late Transition Metals: Assessing Energetic and Structural Properties Thomas M. Soini,† Alexander Genest,‡ Astrid Nikodem,† and Notker Rösch*,†,‡ †

Department Chemie and Catalysis Research Center, Technische Universität München, 85747 Garching, Germany Institute of High Performance Computing, Agency for Science, Technology and Research, 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632

‡

S Supporting Information *

ABSTRACT: We present the first application of hybrid density functional theory (DFT) methods to larger transition-metal clusters. To assess such functionals for this class of systems, we compare the performance of three modern hybrid DFT methods (PBE0, TPSSh, M06) and their semilocal counterparts (PBE, TPSS, M06L) regarding average bond distances and binding energies per atom for a series of octahedral model clusters Mn (M = Ni, Pd, Pt; n = 13, 38, 55, 79, 116). With application to large particles in mind, we extrapolated the results to their respective bulk limits and compared them to experimental values. In some cases, average nearest-neighbor distances are notably overestimated by the PBE0 and M06 hybrid functionals. Results on energies allow a grouping of the tested functionals into sets of similar behavior for the three metals studied. Among the methods examined, the TPSSh hybrid density functional shows the best overall performance. parametrization24,25 that also accounts for data of small transition-metal compounds. This approach is exemplified by meta-GGA (mGGA) functionals whose additional functional dependency on the kinetic energy density promises even higher accuracy.26,27 Nevertheless, the aforementioned static correlation error (SCE) remains present also in these new hybrid DFT methods. Of course, one must consider the clear advantages of hybrid DFT methods over semilocal functionals, particularly regarding reduced self-interaction artifacts.19,24,28,29 This leads to the question whether the impact of the SCE is small enough in more recent hybrid DFT methods to admit a sufficiently accurate treatment of the electronic structure of transitionmetal particles. This situation motivated the present work on model clusters of late transition metals. In the spirit of earlier studies,30,31 we employed model clusters32,33 of various sizes to compare the performance of several (hybrid) DFT methods for structural, energetic, and electronic properties. Cluster scaling techniques permit one to extrapolate such results to their respective bulk limits and to compare these latter data to experiment.34−39 We selected three popular semilocal density functionals and compared their results to their respective hybrid DFT counterparts. The PBE40 GGA functional does not rely on any fitted parameters, and neither does its hybrid version, PBE0, whose (exact) exchange mixing factor of 0.25 is based on theoretical considerations.23 The TPSS functional26 is a

1. INTRODUCTION Density functional theory1 (DFT) within the Kohn−Sham (KS) framework2 has seen many attempts and suggestions for improvement. Hybrid DFT, the admixture of exact exchange (EXX) to a semilocal exchange functional, is a very striking example of such an advance.3,4 Initially, this scheme was empirically motivated;3 theoretical justifications followed.5−9 The B3LYP (hybrid) functional10 soon emerged as the most widely applied DFT method.11 With increasing popularity, it also became apparent that B3LYP offers limited accuracy for systems with bonds between transition-metal centers,12−15 inferior to the computationally more efficient semilocal functionals that rely on the generalized gradient approximation (GGA). These problems were initially rationalized16 by reference to the EXX term that partially substitutes the semilocal exchange part. While the latter implicitly accounts for static correlation,16,17 the EXX term does not, because of its roots in a purely single-determinant theory. The implicit treatment of static correlationhence, the effect of almostequal contributions of several determinants to an electronic stateis therefore partially reduced in hybrid DFT methods, which may deteriorate their accuracy for metals, which are dominated by this effect.17−20 However, a recent analysis21 revealed the EXX term to be only a smaller source of error for transition metals.13 Some missing exact constraints in the LYP correlation term and the exclusive parametrization of B3LYP on main group data4,22 were found to contribute notably more.18 Newer hybrid DFT functionals try to avoid these latter deficiencies either by a less empirical construction23 or a © 2014 American Chemical Society

Received: August 3, 2014 Published: September 18, 2014 4408

dx.doi.org/10.1021/ct500703q | J. Chem. Theory Comput. 2014, 10, 4408−4416

Journal of Chemical Theory and Computation

Article

nonempirical meta-GGA, whereas its hybrid version TPSSh41 was obtained by optimizing its EXX prefactor, 0.10, on a set of experimental data. Finally, the functionals M06L42 and M0624,25 are examples for mGGA and hybrid mGGA functionals with a relatively large number of parameters that were assigned by comparison with several reference datasets. M06 has an EXX prefactor of 0.27, and M06L was obtained as its local reparameterization. Thus, in contrast to the other pairs of functionals examined in this work, M06L and M06 differ in the parameter values of their semilocal parts. While hybrid DFT methods with even larger EXX prefactors exist,24 we will denote PBE0 and M06 in the following as functionals with a high fraction of exact exchange, to distinguish them from TPSSh.

2.2. Cluster Models and Scaling Procedure. Suitable experimental references are hardly available for metal particles. Rather, comparison with calculated results often is restricted to the bulk material. To this end, we will exploit scaling laws based on cluster nuclearity.34−36,60 To achieve good linear scaling behavior, we selected cluster models Mn (M = Ni, Pd, Pt) of close structural similarity with cut-outs of the corresponding fcc bulk metal, namely clusters of truncated cuboctahedral shape. As we are aiming for extrapolation to the bulk, this choice is preferable although more stable isomers likely exist.37,61,62 Similarly, anisotropies and electronic situations with no relevance to the bulk are to be avoided in the present context. Therefore, we selected as models only clusters of Oh symmetry, of closed geometric shells and with a limited number of lowercoordinated atoms. We did not include octahedral clusters such as M19, because of their low-coordinated corner atoms. Guided by these principles, one arrives at truncated octahedra, bounded by (111) and (001) facets, but without four-coordinated corner atoms. We considered cuboctahedral cluster models Mn with up to 3 shells of atoms, n = 13, 38, 55, 79, and 116 (Figure 1). M13, M55, and M79 are atom-centered clusters, and M38 and M116 have an octahedron M6 as central moiety.

2. THEORETICAL METHODS 2.1. Computational Details. All results were obtained by calculations with the linear combination of Gaussian-type orbital fitting-functions density functional (LCGTO-FF-DF) method.43 The computations were carried out with the program package PARAGAUSS, version 4.0,44 which recently was extended by options for exact exchange, hence, hybrid DFT functionals. Although PARAGAUSS includes its own electron repulsion integral routines, in the present study, an external, highly optimized implementation45 was employed in combination with our recently developed dynamic load balancing library to achieve an efficient parallelization.46 We employed orbital basis sets of the type def2-TZVP,47 which, for the elements Pd and Pt, include the corresponding Stuttgart−Dresden small-core pseudo-potentials.48 For Ni, the def2-TZVP basis set describes all electrons. For the evaluation of the Hartree potential43,49 the electronic charge density was approximated with an auxiliary basis set.50,51 That standard set was extended by r2-type functions,43 but their effect on total energies was found to be negligible in preliminary tests on Pt13 and Pt38, suggesting near completeness of the auxiliary basis set. The semilocal exchange-correlation terms were evaluated by numerical integration. The corresponding grid was constructed as a superposition of spherical, atom-centered Lebedev-type grids.52−54 The angular parts of these grids were locally exact for spherical harmonics with angular momentum up to L = 29; the radial part for Ni, Pd, and Pt consisted of 305, 315, and 305 radial shells, respectively. When assembling the exact exchange contribution to the Hamiltonian, a density-matrix-weighted Schwarz screening procedure55 was applied with a tolerance of 10−10 a.u. The self-consistent field (SCF) iterations were converged using the direct inversion of iterative subspace method,56 until the relative change of the elements of the density matrix was below 5 × 10−6. To facilitate the convergence, a fractional occupation number technique43 was invoked using a Fermitype level broadening. During geometry optimization, the energy range of the broadening function was reduced from 0.5 eV to finally 0.05 eV, to help selecting the energetically lowest electronic configuration. From a comparison with values extrapolated to integer occupations in an earlier study57 (section 3.2), we find our results to be converged within ∼5 kJ/mol or less, with respect to the broadening parameter. Hence, we omitted an extrapolation to integer occupations. All cluster structures were allowed to relax under Oh spatial symmetry constraints until the maximum gradient acting on each atomic center was below 10−5 a.u. For this task, the molecular-dynamics-based optimizer FIRE58 was applied, available in the utility suite ParaTools.59

Figure 1. Sketches of the cluster models Mn (M = Pt, Pd, and Ni; n = 13, 38, 55, 79, and 116).

In analogy to the geometric features of the cluster models, one may also strive for unifying features of the electronic structures of the clusters in a series to be used for extrapolation. However, constraining the magnetic moments per atom of the cluster models to the corresponding bulk values (i.e., zero for Pt and Pd, 0.616 for Ni) yields variations of up to 5 pm for the extrapolated nearest-neighbor average distances and 30 kJ/mol for the cohesive energies. This is to be expected because at least the smaller model clusters still behave as molecular entities. Furthermore, the differences were found larger for the results of hybrid DFT calculations, because of their increased band gaps that enhance the aforementioned molecular behavior. Therefore, only results of unrestricted Kohn−Sham calculations will be discussed in the following. Results from spin-restricted 4409



Article

Table 1. Extrapolated Cluster Values of Average Nearest-Neighbor Distances (dav) Obtained for Bulk Metals Ni, Pd, and Pt from Calculations with the Six Functionals Examined (Also Shown Are the Corresponding Values of the Fitted Slopes (kd) and the Coefficients of Determination R2, as Well as Experimental Values)

a

PBE

PBE0

TPSS

TPSSh

M06L

M06

exp

Ni

dav (pm) kd R2

252.8 −13.6 0.946

258.7 −23.4 0.946

251.1 −11.8 0.944

252.0 −12.1 0.950

253.3 −12.4 0.959

258.9 −23.1 0.948

249a

Pd

dav (pm) kd R2

279.8 −23.9 0.984

280.8 −30.0 0.993

276.7 −22.0 0.984

276.9 −23.0 0.971

279.3 −20.4 0.969

288.1 −31.9 0.935

275b

Pt

dav (pm) kd R2

283.9 −33.9 0.987

281.6 −31.6 0.926

281.9 −32.5 0.988

281.0 −31.8 0.984

284.6 −31.4 0.989

290.9 −37.4 0.996

277c

Data taken from ref 67. bData taken from ref 68. cData taken from ref 69.

calculations are presented and briefly discussed in the Supporting Information (SI) (see Section S1). The calculations of atomic species, required for calculating cohesive energies, were carried out in C2v symmetry to admit a localized occupation of all valence shells. The resulting atomic electron configurations were found to be 3d94s1 for Ni, 4d10 for Pd, and 5d96s1 for Pt. The latter two results match the corresponding experimental findings.63 In the case of Ni, the experimental configuration 3d84s2 is known to be preferred bceause of a spin−orbit interaction64 that is not considered here. The extrapolation to the bulk limit exploits the finding that many properties X(n) of metal clusters Mn were found to scale quite nicely with n−1/3 in linear fashion:31,34−36,60,65 X(n) = kXn−1/3 + X∞

(1)

With a least-squares fit of the calculated data X(n) to this relationship, one is able to identify the slope kX as well as the predicted bulk limit X∞ for n → ∞.31,34−36,60,65 This relationship is expected to hold for sufficiently large cluster sizes n, from a minimum nscal onward, for quantities that scale with the surface-to-volume ratio of the particles, or the inverse of the effective particle radius or n−1/3.34−36,60 The specific value of nscal is dependent, of course, on the properties and the systems examined.33,65,66

Figure 2. Average nearest-neighbor distances dav (pm) calculated for the model clusters Ptn together with the fitted trend line linear in n−1/3.

study,31 which agree extremely well, within 0.1 pm. For most exchange-correlation approximations, the coefficients of determination (R2) increase along the series Ni < Pd < Pt; PBE0 represents an exception, because of its higher R2 value of 0.993 for Pd. Otherwise, R2 values typically are ∼0.95 for Ni, ∼0.98 for Pd, and generally even higher for Pt. From these results, one concludes that extrapolation of the average bond lengths to the corresponding bulk limits is generally well-justified. The bulk limits dav can generally be ordered as follows (see Table 1): d(TPSS) ≈ d(TPSSh) < d(PBE) ≈ d(M06L) < d(PBE0) ≈ d(M06). Only for Pt, PBE0 yields a comparatively low limit for dav, because of minor outliers in this data series. Comparing the semilocal DFT methods with their hybrid DFT counterparts, one finds that the functionals PBE0 and M06 (i.e., the functionals with a high fraction of single-determinant exchange) mostly overestimate bonds lengths. The extrapolated results dav(PBE0) are by 1.0 pm (Pd) and 5.9 pm (Ni) larger than the corresponding PBE values. Compared to M06L, the M06 method yields extrapolated bulk distances that are higher−by 5.6 pm for Ni, 6.0 pm for Pd, and 5.2 pm for Pt. Only the aforementioned PBE0 result for Pt is 2.1 pm smaller than the corresponding PBE result. This significant impact of the exact exchange term on structural properties can be

3. RESULTS AND DISCUSSION 3.1. Average Bond Length. Average nearest-neighbor distances dav of moderately large metal clusters are known to increase smoothly with cluster size.31,35,36,38 This behavior reflects the trend to decreasing bond distances when the average coordination number decreases because larger clusters comprise a smaller fraction of surface atoms with lower coordination numbers, hence featuring stronger interactions with their neighboring atoms.35,36 Apart from reflecting such physical variations, these changes in geometry also provide first insight into the quality of the extrapolation procedure (Table 1). Results for individual clusters are collected in Table S3 of the SI, while Figure 2 shows the corresponding linear function (eq 1) on the example of Pt. A graphical comparison of the data for all three metals and the corresponding fitting functions are provided as Supporting Information (Figure S1). As expected, the dav values increase with the nuclearity of the cluster models. We compare the PBE estimate for Pd also to the corresponding results of a recent all-electron cluster scaling 4410



Article

Table 2. Extrapolated Cluster Values of the Cohesive Energy Ecoh per Atom Obtained for Bulk Metals Ni, Pd, and Pt from Calculations with the Six Functionals Examined (Also Shown Are the Corresponding Values of the Fitted Slopes (kE) and the Coefficients of Determination (R2), as Well as Experimental Values)

a

PBE

PBE0

TPSS

TPSSh

M06L

M06

exp

Ni

Ecoh (kJ/mol) kE (kJ/mol) R2

444.6 −405.5 0.993

303.6 −304.3 0.994

459.4 −452.6 0.994

393.5 −409.6 0.996

462.2 −491.4 0.998

347.8 −387.8 0.981

428a

Pd


395.9 −381.3 0.996

337.1 −431.7 1.000

425.0 −427.5 0.997

387.8 −410.2 0.998

440.3 −478.7 0.998

315.4 −370.7 0.999

376b

Pt


555.9 −512.1 0.995

494.5 −530.3 0.991

580.1 −566.0 0.995

546.8 −547.5 0.996

579.9 −583.6 0.997

487.1 −470.0 0.995

563b

Data taken from ref 74 bData taken from refs 72 and 73.

Therefore, except for Pt, we are unable to confirm a good performance of PBE0, as reported for the bulk metals in recent plane-wave studies.21,70 Although the corresponding dav values are mostly smaller than the present results, the overall agreement between the two datasets still seems reasonable. 3.2. Cohesive Energy. The cohesive energy Ecoh per atom of a cluster Mn is determined from the total energies Etot(Mn) of the cluster and Etot(M1) of the corresponding atomic species:

rationalized by the fact that hybrid DFT functionals prefer other electronic states than semilocal functionals, as will be discussed in section 3.3. These states subsequently lead to different structures. Compared to all other methods, the slopes (kd) of the linear trend lines of PBE0 and M06 for Ni and Pd are notably more negative (see Table 1). The Pdn structures obtained with PBE0 are less expanded than the corresponding PBE and M06L geometries, so that the overestimation of the bulk limit results entirely from the steeper slope. The situation is different for TPSSh where the EXX contribution is smaller. This functional yields extrapolated dav values that are very similar to the corresponding TPSS results, with differences of

Hybrid Density Functionals for Clusters of Late Transition Metals

Recommend Documents