Benchmarking Density Functional Based Tight ... - ACS Publications

Article pubs.acs.org/JPCA

Benchmarking Density Functional Based Tight-Binding for Silver and Gold Materials: From Small Clusters to Bulk Luiz F. L. Oliveira,†,‡ Nathalie Tarrat,§ Jérôme Cuny,† Joseph Morillo,§,∥ Didier Lemoine,‡ Fernand Spiegelman,† and Mathias Rapacioli*,† Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, UMR5626 and ‡Laboratoire de Collisions Agrégats et Réactivité LCAR/IRSAMC UMR5589, Université de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France § CEMES CNRS UPR 8011, 29 rue Jeanne Marvig, BP 94347, 31055 Toulouse Cedex 4, France ∥ Université de Toulouse (UPS), 118 Route de Narbonne, F-31062 Toulouse Cedex 9, France Downloaded via EASTERN KENTUCKY UNIV on January 16, 2019 at 23:01:18 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: We benchmark existing and improved self-consistent-charge density functional based tight-binding (SCC-DFTB) parameters for silver and gold clusters as well as for bulk materials. In the former case, our benchmarks focus on both the structural and energetic properties of small-size AgN and AuN clusters (N from 2 to 13), medium-size clusters with N = 20 and 55, and finally larger nanoparticles with N = 147, 309, and 561. For bulk materials, structural, energetics and elastic properties are discussed. We show that SCC-DFTB is quite satisfactory in reproducing essential differences between silver and gold aggregates, in particular their 2D−3D structural transitions, and their dependency upon cluster charge. SCC-DFTB is also in agreement with DFT and experiments in the medium-size regime regarding the energetic ordering of the different low-energy isomers and allows for an overall satisfactory treatment of bulk properties. A consistent convergence between the cohesive energies of the largest investigated nanoparticles and the bulk’s is obtained. On the basis of our results for nanoparticles of increasing size, a twoparameter analytical extrapolation of the cohesive energy is proposed. This formula takes into account the reduction of the cohesive energy for undercoordinated surface sites and converges properly to the bulk cohesive energy. Values for the surface sites cohesive energies are also proposed.

■

use of accurate functionals such as revTPSS,47 M06,49,50 SOGGA,51 or RPA47 methods seems to have lifted the uncertainties about the structural and electronic features of even small-size clusters. For instance, the 2D−3D transition in gold clusters has been in debate for some time and only recent theoretical studies provided results in agreement with the experimental data.47 This transition was shown to depend on the cluster’s charge and to occur earlier for gold cations (N = 8) than for neutrals (N = 11) or anions (N = 12). Another example is Ag20, for which recent publications still propose different structure orderings.27,32 Force fields, generally based on manybody potentials,52−58 allow for global structural investigations of metal clusters with N as large as a few hundred particles59−66 and local minimization investigations of systems up to several thousands (tens of thousands) particles.67−69 However, they do not allow us to describe electronic properties such as ionization, electron attachment, or Jahn−Teller effects. They have, nevertheless, been fruitful in providing size-dependent structural transitions,67,68,70−72 and also for selection of candidate structures to be further optimized at a higher level of theory.27,73 Despite fast progress in the derivation of accurate functionals and efficient algorithms to be implemented on high performance

INTRODUCTION Silver and gold materials, both noble metals of group 11, have been largely studied due to their large range of applications in materials, surface science, and nanoscience. For instance, silver and gold nanoparticles exhibit some unique chemical and physical properties such as photoabsorption,1−5 fluorescence,6−8 and catalytic activity,9−15 which can be applied, among others, for biological labeling,16−19 and for the development of light emitting sources in optoelectronics.20−23 Fine tuning of these systems properties, within bulk, surfaces, and aggregates, are closely related to their geometrical arrangements and electronic structure details that are hardly accessible by experimental measurements only. Consequently, one has to resort to quantum chemical simulations to provide a better understanding of experiments such as vibrational or photoelectron spectroscopy to help in ascertaining the structure of silver and gold clusters. This has been achieved essentially either via wave function based methods for the characterization of small clusters24−27 or via density functional theory (DFT) for bulks, surfaces, and middlesize clusters. In particular, the structure of neutral, cationic, and anionic clusters with N up to ∼20−25 has been the subject of numerous DFT investigations for both silver22,27−32 and gold.33−48 Nevertheless, although DFT can handle small to middle-size aggregates, the resulting properties show significant dependence over the functional, and it is only recently that the © 2016 American Chemical Society

Received: September 23, 2016 Published: September 26, 2016 8469

DOI: 10.1021/acs.jpca.6b09292 J. Phys. Chem. A 2016, 120, 8469−8483

Article

The Journal of Physical Chemistry A

which converges properly to the bulk cohesive energy. Conclusions and perspectives are provided in the last section.

computers, DFT is not yet efficient enough to be directly incorporated into global and extensive exploration of potential energy surfaces (PES). Density functional based tight-binding (DFTB) is a quantum mechanical method that is derived from DFT via several approximations.74−79 It is a quantum formalism that provides a balance between accuracy, computational efficiency, and the ability to explicitly take into account the electronic structure. It has proved to be particularly efficient for the description of complex molecular systems.80−84 Furthermore, DFTB has previously been applied to investigate the structural and dynamical properties of small gold anionic clusters, either free or supported on films and substrates.84−87 DFTB is a parametrized method and its accuracy may depend significantly on the quality of the parametrization. In particular, it requires the use of parametrized integrals, generally referred to as the Slater−Koster integrals, constructed for each atom pair of the chemical system of interest. These Slater−Koster integrals can be developed using various methodologies and reference DFT functionals and are not uniquely defined. In the particular case of silver and gold, several SCC-DFTB parameter sets have been proposed in the literature. For gold, Koskinen and co-workers developed a parametrization achieving a near DFT-like accuracy for the description of small anionic gold clusters AuN− up to N = 14.85 In a latter study by Fihey et al., another set of Slater− Koster integrals was developed for Au−Au, as well as integrals for all the Au−X (H, C, S, N, O) atomic pairs, to describe gold− thiolates hybrid systems.88 For silver, Szűcs et al. simulated the tunnelling current between a silver tip and a sulfur-passivated GaAs substrate using the SCC-DFTB formalism.89 The Ag−Ag parametrization used provided satisfactory results, nevertheless, and to our knowledge, it has not been further benchmarked for the description of other silver systems. The SCC-DFTB parameters should be, in principle, transferable to various chemical systems. However, it has been demonstrated by Wahiduzzaman et al. in the particular case of carbon that it is not always possible to develop a SCC-DFTB potential that is transferable to chemically different phases of the same element.90,91 Consequently, the Au−Au and Ag−Ag parametrization mentioned above are not necessarily accurate for the description of clusters, nanoparticles, surfaces, and bulks. The main goal of the present paper is to check and possibly improve the reliability and transferability of the available DFTB parameters for silver and gold, examining the whole domain from small to large clusters and for the bulk. In the smalland medium-size regime, we will check the ability of DFTB to reproduce, with a unique parametrization, details of the geometrical and energetic properties and electronic structures to account for the important differences between silver and gold aggregates in their neutrals, cations, and anions. At the same time we will check its ability to reproduce bulk structural and energetic properties and the convergence of large cluster energetics toward the bulk limit. The paper is organized as follows: First, the methodology and the SCC-DFTB parameters discussed along this paper are presented. Then, a comparison between DFT and SCC-DFTB results for small gold and silver clusters (N ≤ 13) is given, testing various SCC-DFTB parametrization sets with a particular attention to the 2D−3D transition. Bulk properties (structure, relative stabilities, and elastic constants) are then presented. Finally, medium-size clusters (N = 20 and 55) are investigated and we propose a two-parameter analytical extrapolation of the cohesive energy of large clusters (N = 147, 309, and 561),

■

METHODOLOGY AND COMPUTATIONAL DETAILS Several reviews on the DFTB and SCC-DFTB methods can be found in the literature.74−79 SCC-DFTB mostly differs from the Kohn−Sham DFT expressed in a local basis set by the following approximations: (i) the DFT energy is expanded up to the second-order with respect to electronic density fluctuations around a given reference density, (ii) all three centers interaction integrals are neglected, (iii) the molecular orbitals are expressed in a minimal atomic basis set, (iv) the short distance repulsive potential Erep is expressed as a function of two body interactions, and (v) the second-order term in the DFT energy expansion is expressed as a function of atomic charges and a Γ matrix, taking into account the charge density fluctuations. With these approximations, the total SCC-DFTB energy reads atoms

E SCC ‐ DFTB =

∑ α ,β≠α

rep Eαβ +

∑ ni⟨φi|Ĥ i

0

|φi⟩ +

1 2

atoms

∑

Γαβqαqβ

α ,β

(1)

̂0

where H is the Kohn−Sham operator at the reference density. The nondiagonal elements of the overlap and Ĥ 0 matrices in the atomic basis, as well as Erep αβ are interpolated from DFT calculations on atomic dimers and expressed as a function of interatomic distances. The diagonal elements of the Hamiltonian matrix are the orbital energies of the isolated atoms. Finally, Γαβ terms are derived from Hubbard parameters. ni is the occupation number of molecular orbital φi, and the qα are the atomic charges (Mulliken). In the following, the SCC-DFTB approach will simply be referred to as DFTB. Calculations involving clusters have been performed with the deMonNano code.92 The strategy followed for the structural study of small clusters in the range N ≤ 13 was 2-fold. (i) We achieved a global optimization search using parallel tempering canonical molecular dynamics (PTMD) followed by quenching. Precisely, the PTMD scheme involved simultaneous non-SCC MD runs for 60 temperatures equally distributed in the range 0−5000 K, allowing configurational exchanges between trajectories. Fragmentation was hindered by using a spherical box of radius 10 bohrs. The MD time step was set to 3 fs, and exchanges were attempted using a Boltzmann energy criterion every 10 timesteps. A chain of 5 Nosé−Hoover thermostats with a unique frequency 80 cm−1 was used. The lengths of the trajectories were 3 × 105 fs. From each of the three respective MD runs at 83 K (low T), 1000 K (medium T), and 3000 K (high T) of the PTMD process, 125 configurations equally spaced in time along the simulation were selected and optimized locally with the SCC-DFTB conjugated gradient scheme, providing the results for the final stable structures. This generated a bunch of low energy isomers for each size and was repeated for neutrals, anions, and cations. (ii) In addition, we have also checked some extra geometrical configurations found by previous authors and not present in the final yields of our own global optimization scheme. Actually, the global minimum was always found with the above PTMD+quenching scheme and in a few cases, new candidates were found with respect to the literature data. 8470


Article

The Journal of Physical Chemistry A The quality of the DFTB results will be probed by benchmarking the corresponding properties against reference calculations: (i) the evolution of the binding energy per atom as a function of the size defined as E b(N ) = (NE1 − EN )/N

1 (EN + 1 + EN − 1 − 2EN ) (3) 2 This quantity is relevant because it allows us to evidence particularly stable structures, also known as magic numbers structures. (iii) the size corresponding to the 2D−3D structural transition. For bulk materials, our computational investigations were conducted using the DFTB+ simulation package.93,94 The periodic simulation cell was a N = 4 atoms fcc Bravais cubic cell. A 40 × 40 × 40 Monkhorst−Pack mesh of k-points was used to sample the Brillouin zone. For computing the cohesive energies of the hcp bulk, a periodic rhombohedral simulation cell with N = 2 atoms was used together with a 40 × 40 × 40 Monkhorst−Pack mesh of k-points. To evaluate the performance of the DFTB parameters regarding the bulk properties, we calculated the equilibrium lattice parameter a0, the binding energy Eb, and the three independent cubic elastic constants B0, C′, and C44 of the Au and Ag fcc structures and finally checked the stability of the fcc structure against the hcp structure. The bulk modulus B0 represents the response of the crystal to an isotropic stress (homogeneous pressure P) whereas C′ and C44 its response to plane shear stresses on the (010), (110) plane C in the [100], [11̅0] direction, respectively. Their ratio z = C44′ (the Zener ratio) is thus a measure of the anisotropy of the cubic crystal. B0, C′, and C44 must be positive for the crystal stability against any deformation. A nice agreement between calculated and experimental values of a0, Eb, B0, C′, C44, and z is thus a good criterion of the ability of the DFTB potential to describe the equilibrium cubic crystalline structure of the material under consideration and of its anisotropy. The fcc and hcp close-packed structures have the same first and second neighbor shells of atoms and are thus very close in energy. A proper description of the relative stability of these two structures is thus a stringent test of the accuracy of any interaction model to describe subtle energy changes in the bulk properties of close packed structures. a0 and B0 at T = 0 K and P = 0 Pa were deduced from a fit of the T = 0 K energy versus volume curve using the Murnaghan equation of state (Supporting Information).95 C′ and C44 were obtained through a sixth-order polynomial fit of the elastic energy as a function of the corresponding applied shear strain δ (Supporting Information). The binding energy Eb per atom for both the fcc and hcp structures, defined as the energy gain in the crystal formation from the isolated atoms, was calculated using Δ 2 E (N ) =

E0 N

3B0 + 4C′ 3

(5)

C12 =

3B0 − 2C′ 3

(6)

Note that, to check the validity of our results, the C11 constant has also been computed through a sixth-order polynomial fit of the elastic energy as a function of a uniaxial tensile/compressive strain applied in the [100] direction. In both cluster and bulk calculations, the basis set includes s, p, and d shells, as standardly used in the description of noble and transition metal systems which are likely to undergo electronic valence transitions.

(2)

where EN is the energy of a cluster of N atoms. (ii) the second-order derivative of the total energy

E b = Eatom −

C11 =

■

SMALL CLUSTERS Gold Clusters. The reference data chosen for AuN clusters in the range N = 1−20 are taken from DFT calculations performed with the functionals LDA,96 GGA,29 revTPSS, and RPA.47 In the present study, we will benchmark the parameters from Fihey et al.88 (hereafter referred to as the first set of parameters, or DFTBα) and those from Koskinen et al.85 (referred to as the second set of parameters or DFTBβ). With these DFTB parameters, we optimized all the structures from Johansson et al.47 and Fernández et al.,29 in addition to the structures obtained from global optimization search. As will be discussed in the following, the DFTBα set provides reasonable results for bulk properties but not for clusters in term of the 2D−3D transition whereas the opposite is observed for the DFTBβ set. We noticed that a major difference between these two parameters relies on the energy attributed to the p orbital (−0.0279 Ha in the first case, −0.0096 Ha in the second case). Therefore, we introduced a third set of parameters consisting of DFTBα parameters with a significant shift of the p orbital energy to a value of −0.00001 Ha). This will be referred to as the third set, or DFTBχ. Figure 1 compares the binding energies obtained for the most stable clusters with different levels of theory. We note that revTPSS and RPA calculations on gold clusters are not presented in this figure as the authors only reported relative energies.47 Although the values reported at the DFT level are in fair agreement for the dimer, the binding energies obtained for larger sizes with the LDA functional exceed those published using GGA. At the DFTB level, the general trend for the binding energy increase is pretty well reproduced by the three sets of parameters. We notice, however, that although the first and third sets of parameters provide dimer binding energies in agreement with ab initio data, these are overestimated when the second set is used. For the larger sizes, the binding energies of the first and second sets lie between the two DFT values. Finally, the third set of parameters reproduces pretty well the GGA results. The cluster size corresponding to the 2D−3D transition in gold clusters has been widely investigated at the DFT level and attributed to relativistic effects.97 The most accurate calculations at the moment performed at the revTPSS and RPA level indicate that 2D and 3D structures are almost isoenergetic for clusters of 11 units.47 Figure 2a reports the binding energy of the most stable 2D and 3D structures of AuN obtained at the DFTB level. It can be seen that the transition occurs much too early with the first set of parameters (3D structures are clearly more stable for sizes above 6 units). With the DFTBβ parameter set, this transition is shifted toward larger sizes in

(4)

where Eatom is the energy of the isolated atom in vacuum and E0 is the energy of the bulk unit cell with N atoms. Finally, from the B0, C′, and C44 values, we deduced the z ratio and the elastic constants C11 and C12: 8471


Article


note that the 2D structures remain competitive even for sizes of 13 units. Interestingly, with the third set of parameters, the transition seems to occur within the correct size range with a 2D−3D competition for sizes from 8 to 11 units. Clearly, the shift of the p orbital (the only difference between the first and third sets of parameters) to higher energy values is responsible for preserving the planar geometry at large sizes. This increases the gap between the d and p orbitals and results in a larger population of d orbitals with respect to p orbitals, the first ones favoring planar structure and the latter ones 3D geometries. We note that DFT results also correlate the conservation of AuN planar geometries up to larger sizes with a larger d−p gap.29 Let us also mention that the choice of the atomic configuration to extract p orbitals (virtual) versus d and s orbitals (occupied) remains an open question in DFTB parametrization. To further document the validity of the here proposed DFTBχ parametrization, we also report a study of singly charged clusters AuN+ and AuN− (see illustrations and Cartesian coordinates in Supporting Information). Figure 2b shows that the 2D−3D transition occurs for N = 8, 9 in the case of cations, versus N = 12, 13 in the case of anions, which is consistent with recent DFT studies using expectedly accurate functionals,47 namely, the transition size increases from cations (N = 8) to neutrals (N = 11) and then anions (N = 12). The most stable structures obtained with the DFTBχ parameters are reported for neutrals in Figure 3, and the corresponding binding energies are given in Table 1.

Figure 1. Evolution of the DFTB binding energy of gold clusters as a function of the cluster size computed with three different DFTB parameter sets (DFTBα,β,χ, see text). DFT results have been added for comparison: LDA results from ref 96 and GGA results from ref 29.

Figure 3. Lowest energy geometries of AuN clusters (DFTBχ parameters). First row: 2D structures (N ≤ 7). Second and third rows: 2D and 3D structures (8 ≤ N ≤ 11). Fourth row: 3D structures (N = 12 and 13). (Associated Cartesian coordinates are given in the Supporting Information.)

Figure 2. Evolution of the DFTB binding energy of gold clusters as a function of the cluster size for the most stable 2D and 3D structures. (a) Neutral clusters computed with three different DFTB parameter sets (DFTBα,β,χ, see text). (b) Cationic and anionic clusters computed with the DFTBχ set.

Figure 4 reports the second energy differences Δ2E for DFT and DFTB results. It can be seen that all DFTB parametrizations yield the odd−even alternation of the second energy derivative with size, characterizing the most simple metal clusters. The hexamer appears as an electronic magic number at the DFT level, as a result of the planar geometry pattern

better agreement with DFT. We notice that the correct 2D−3D transition was already reported for these parameters in the context of negative and neutral gold clusters.85,87 We, however, 8472


Article

The Journal of Physical Chemistry A Table 1. Binding Energies, Eb (eV/at) of the Most Stable 2D and 3D Neutral Clusters with the DFTBγ (Silver) and the DFTBχ Parameters (Gold) size

AgN(2D)

2 3 4 5 6 7 8 9 10 11 12 13

1.03 1.09 1.37 1.46 1.62 1.59 1.65 1.69 1.78 1.76 1.83 1.81

AgN(3D)

AuN(2D)

AuN(3D)

1.42 1.64 1.74 1.86 1.84 1.90 1.92 1.96 1.98

1.59 1.70 2.10 2.27 2.55 2.51 2.63 2.63 2.72 2.72 2.80 2.79

2.12 2.35 2.49 2.62 2.63 2.72 2.73 2.81 2.82

these two sets will be benchmarked against DFT calculations reported in the literature that were performed using TPSS,32 PBE,29 M06,49 B97,49 and PW9149 functionals as well as against CI24 and CCSD(T).49 The latter can be considered as the highest level of calculation available for these systems. In addition to the structures obtained from global optimization search, the most stable 2D and 3D structures reported by Fernández et al.29 and Chen et al.27 were optimized within the DFTB framework. As can be seen from Figure 5, TPSS, M06, B97, and PW91 functionals, CI and CCSD(T) calculations give consistent

Figure 5. Evolution of the DFTB binding energy of silver clusters as a function of the cluster size computed with two different DFTB parameter sets (DFTBϕ,γ, see text). DFT (TPSS,32 PBE,29 M06,27 B97,27 PW9127), CCSD(T),27 and CI24 results have been added for comparison.

binding energies within the investigated range, whereas the PBE functional, although providing the correct trend for the size evolution, seems to significantly overestimate the binding energy versus the previous methods. Interestingly, the same observation can be made for the DFTB results with the first set of parameters, the binding energy overestimation being even larger. The modified parameters strongly reduce this overestimation, the binding energies being even closer to the CI ones. In the literature, the 2D−3D transition is usually reported to lie between 6 and 7 units.27 At the DFTB level, it can be seen from Figure 6 that, with the DFTBϕ set of parameters, the transition occurs too early, namely, below 5 units. On the contrary, the transition occurs at 6 units with the modified DFTBγ parameters. The corresponding optimized structures and binding energies are reported in Figure 7 and Table 1. For the second energy derivatives (Figure 8), all DFT calculations predict particular stabilities for clusters containing an even number of units. The same is observed at the CI level apart from the heptamer that displays a value of Δ2E larger than the hexamer. The latter singularity is not present in the recent CCSD(T) results. At the DFTB level with the first set of parameters the odd−even behavior is respected except in the 5−8 size range with Δ2E(5) < Δ2E(6) < Δ2E(7) < Δ2E(8).

Figure 4. Second energy derivative as a function of cluster size for Au clusters. The ab initio methods are the same as referenced on Figure 1 and DFTB is performed with three sets of parameters (α, β, and χ on top, middle, and bottom patterns, respectively).

(possibly interpreted in the 2D quantum dot model, see ref 41). This feature is well reproduced at the DFTB level with the second and third parameters whereas it is less pronounced with the first set of parameters. This can be correlated with the too early appearance of 3D geometries for the heptamer with DFTBα. Silver Clusters. The DFTB parameters tested for silver belong to a set originally developed in the context of surface modeling by Szűcs et al.89 (hybride set from www.dftb.org Web site). This latter will be referred to as DFTBϕ. As will be shown in the following, it leads to an overestimation of the silver cluster binding energies and a too early 2D−3D transition. To reach reasonable binding energies, we introduced a slight modification to this parameter consisting of scaling the interatomic DFTB Hamiltonian matrix elements by a factor of 0.9. This modified parameter will be referred to as DFTBγ. The results of 8473


Article


Figure 8. Second energy derivatives as a function of size for AgN clusters. The ab initio methods are the same as those referenced in Figure 5. DFTB is performed with two sets of parameters (ϕ and γ in upper and lower panels, respectively).

The modified DFTBγ parameters present the odd−even pattern in agreement with both DFT and CCSD(T) results. We notice finally that the octamer appears as a magic number whatever the DFTB parametrization used. This is in agreement with all ab initio results. This size is characteristic of electronic magic number of spherical systems (with reference to the 3D spherical jellium model, see the discussion in ref 41). Vertical Ionization Potentials and Vertical Detachment Energies. Finally, we show in Figure 9 and Table 2 two electronic properties that can be almost directly compared with experimental data, namely, the vertical ionization potentials (VIPs, probing the geometries of neutrals) and the vertical detachment affinities (VDEs, probing the geometries of the anions). The optimized geometries of anions and cations cations are shown in the Supporting Information. In the case of silver (left), the VIPs obtained with DFTB show a global underestimation of 0.5 eV with respect to the experimental IPs of Jackschath et al.98 However, despite this shift, the size evolution of the IPs is in excellent agreement with the experimental one, and in particular the odd−even alternation is quite well reproduced by DFTB. An even better convergence is obtained for the VDEs that follow almost perfectly the experimental points of Ho et al.99 The situation is clearly somewhat deteriorated for the VIPs of gold. Although the DFTB results still exhibit an odd−even alternation, this alternation is much smoother than in the experiment.98 Also the DFTB IPs tend to decrease too fast at the larger sizes. This does not seem to be due to the DFTB geometries of neutrals, which are in general good correspondence (with a few exceptions) with those obtained in DFT/B3PW91 calculations.100 On the contrary to VIPs, the evolution of the VDEs is quite satisfactory with respect to the experimental determination of Taylor et al.,101 except for size n = 10, which provides an exception to the odd−even alternation of the VDEs. Let us mention here that the geometries obtained for gold anions are generally the same as those derived using DFT (BP86 functional) by Furche et al.,102 except for N = 5, 6, 7, and 13. For size N = 10, the present geometry is the same as that of Furche et al., who, however, found a much lower VDE from their second anionic isomer than from the lowest one.

Figure 6. Evolution of the DFTB binding energy of silver clusters as a function of the cluster size for the most stable 2D and 3D structures. (a) Neutral clusters computed with two different DFTB parameter sets (DFTBϕ,γ, see text). (b) Cationic and anionic clusters computed with the DFTBγ set. (See illustrations and Cartesian coordinates in the Supporting Information.)

Figure 7. Lowest energy geometries for silver clusters (DFTBγ parameters). First row: 2D structures (N ≤ 4). Second row: 2D and 3D structures (N = 5 and 6). Third and fourth rows: 3D structures (7 ≤ N ≤ 13). (Associated Cartesian coordinates are given in the Supporting Information.) 8474


Article


Figure 9. Size evolution of the vertical ionization potentials (VIP) and vertical detachments energies (VDE) for silver (left) and gold (right) clusters. DFTB is in the full line. Symbols with joining lines stand for the experimental IPs of silver and gold clusters,98 and the experimental VDEs of silver99 and gold101 clusters

Table 2. Vertical Ionization Potential (VIP, eV) and Vertical Electron Detachement Energies (VDE, eV) Obtained for Silver and Gold Cluters (with DFTBγ and DFTBχ Parameters, Respectively) Ag

■

Table 3. Experimental and DFTB Calculated Properties of the fcc Gold Crystala

VIP

VDE

VIP

VDE

2 3 4 5 6 7 8 9 10 11 12 13

6.90 5.46 6.19 5.74 6.42 5.80 6.35 5.01 5.62 5.43 5.77 5.31

1.22 2.42 1.70 2.23 2.19 2.28 1.45 2.29 1.95 2.41 2.00 2.57

9.12 7.91 7.56 7.05 7.42 6.45 6.78 6.58 6.75 5.88 6.50 6.46

1.29 3.39 2.63 3.35 1.76 3.03 1.97 2.81 2.97 2.99 2.16 2.84

DFTBβ

DFTBχ

GGA105

4.0782 2.884106 3.69107 160−180108−110 15.98110 45.44110 201.63110 169.67110 2.84

4.082 2.886 3.52 162 19 (12b) 14 172 (187c) 149d 0.74 3.4

4.401 3.112

4.084 2.889 3.39 162 20 (11b) 13 170 (189c) 149d 0.65 2.4

4.197 2.968 3.21 137.6 12.6b 44.1 154.4 129.2 3.50 1.9

1.635104

1.736

1.736

1.655

2.837

2.838

2.952

exp

a0 d1 Eb B0 C′ C44 C11 C12 z ΔE

106

Au

size

DFTBα

property

hcp

( ac )

hcp (a)

2.96

104

a

Lattice constant a0 (Å). Binding energy Eb (eV/at). First neighbours fcc distance d1 (Å). Elastic constants (GPa): bulk modulus B0, shear elastic moduli C′ and C44, cubic elastic constants C11 and C12. Zener ratio z = C44/C′. ΔE = Eb(fcc) − Eb(hcp) (meV/at). c ratio and a a lattice parameter (Å) of the hcp cell. GGA data from the literature were added for comparison. b Calculated with the formula 3B + 4C ′ d C −C C′ = 11 2 12 cCalculated with the formula C11 = 0 3 Calculated

BULK PROPERTIES Gold Bulk. Experimental and theoretical bulk properties are reported in Table 3. We see that the lattice constant of fcc gold computed with the DFTBβ parameters differs by 8% from the experimental value. By contrast, the lattice constant computed with the DFTBα and the DFTBχ parameters match very nicely (+0.09% and +0.14%, respectively) the experimental one, even better than DFT (+2.9%). We thus exclude the DFTBβ parameters from subsequent tests. The binding energy of the solid is underestimated by all methods. However, despite being too small, the DFTB computed values remain close to experimental ones (−4.6% for DFTBα and −8.1% for DFTBχ), which is not the case for DFT (−13.0%). Concerning the elastic constants, one must keep in mind that, even the most accurate DFT quantum-mechanical calculations provide results that are typically off by 15% with respect to the experimental values103 and this error is often much larger for small elastic constants. The bulk modulus, which experimental value is 170 ± 10 (±6%) GPa, is surprisingly rather well reproduced by both DFTBα and DFTBγ models (−5%), even better than with DFT (−19%). Regarding the shear elastic moduli C′ and C44, both DFTB models give similar poor results. Not only do they differ significantly from the experimental ones (respectively ≈+22% and −70%) but also, because they underestimate C44 and overestimate C′, they lead to a highly underestimated Zener ratio (≈0.7 compared to 2.84) and thus a very bad description of the crystalline

with the formula C12 =

3B 0 − 2C ′ 3

anisotropy. This is not the case for DFT, which slightly underestimates both elastic moduli (−21 and −3%) and thus leads to a Zener ratio (3.50) closer to the experimental one. The performance of the DFTB parameters regarding the relative stability of the fcc and hcp close-packed structures is similar to that of DFT. Concerning the first neighbors distance (respectively d1 (fcc) and a (hcp)) and the c ratio, both DFTB a

models present the same tendencies as DFT (a < d1 and c larger a than the theoretical compact ratio (1.633)). Compared to experimental results, c is overestimated and the first neighbor a variation between fcc and hcp structures is opposite. However, the comparison with the available experimental data104 is not straightforward, because these values were obtained on nanoparticles. As a partial conclusion, the performance of the DFTBα and DFTBχ parameters are comparable and very good regarding the lattice constant, the bulk modulus, the cohesive energy, and the stability of phases and only acceptable regarding the shear elastic constants and the anisotropy. Silver Bulk. Table 4 gathers the different bulk properties computed to evaluate the performance of the DFTBϕ and DFTBγ 8475


Article


Clusters N = 20 and 55. We investigated larger clusters with sizes N = 20 and 55. Those are particularly sensitive tests in the medium-size range because vibrational or photoelectrons spectra were recorded and comparison with DFT calculations were carried out. The structures of four Ag20 isomers are reported in Figure 10. The three lowest energy ones (labeled a, b, and c) are partly

Table 4. Experimental and DFTB Calculated Properties of the fcc Silver Crystala property

106

a0 d1 Eb B0 C′ C44 C11 C12 z ΔE

hcp

exp 4.0857 2.889106 2.95113 108.72110 17.08110 51.09110 131.49110 97.33110 2.99

( ac )

DFTBΦ

DFTBγ

GGA

4.255 3.009 4.07 105

4.067 2.876 3.02 124 24 (30b) 65 167 (156c) 108d 2.71 0.6

4.159111−4.069112 2.941111−2.877112 2.49111−3.445112 83.3111−116112 21.05b 58.1114 161.2114 119.1114 2.76 3112

1.633

1.656112

2.877

2.866112

1.63106

hcp(a) a

Lattice constant a0 (Å). First neighbours fcc distance d1 (Å). Cohesive energy Eb (eV/at). Elastic constants (GPa): bulk modulus B0, shear elastic moduli C′and C44, cubic elastic constants C11 and C12. Zener ratio z = C44/C′. ΔE = Eb(fcc) − Eb(hcp) (meV/at). c ratio and a a lattice parameter (Å) of the hcp cell. GGA data from the literature b were added for comparison. Calculated with the formula 3B + 4C ′ d C −C C′ = 11 2 12 cCalculated with the formula C11 = 0 3 Calculated with the formula C12 =

3B 0 − 2C ′ 3

parameters, GGA data from the literature were added for comparison. The lattice constant of fcc silver computed with the DFTBϕ parameters differs from the experimental value by more than 4%. Its performance is thus nonconvincing for both clusters and bulk. The lattice constant obtained with the DFTBγ model is in fair agreement with the experimental one (−0.46%). The binding energy is slightly overestimated by the DFTBγ model (+2.3%) whereas it is poorly estimated by DFT calculations (≈16%). Concerning the bulk modulus, although it is somehow significantly overestimated (+14%) with the DFTBγ model, it is incorrectly estimated by quite large amounts (−23%111 and +7%112) by DFT. Regarding the shear elastic moduli C′ and C44, the results of the DFTBγ parametrization are acceptable. Both the C′ and the C44 calculated values overestimate the experimental ones (+41% and +27%, respectively), leading to a good Zener ratio (2.71 versus 2.99) and thus to a satisfactory description of the anisotropy. These results are reasonnably similar to the DFT results, which also overestimate both shear constants but to a lesser extent (respectively, 23% and 14%), leading to a similar Zener ratio (2.76). The DFTB computed relative stability of the fcc and hcp close-packed structures is correct. In hcp bulk, the c ratio is very close to experiments a (+0.2%), even closer than the DFT one (+1.6%). The first neighbors distance in hcp bulk does not vary significantly from fcc bulk in both DFTB and DFT calculations. The overall performance of the DFTBγ model relative to the description of bulk properties is thus good, nearly comparable to the DFT one.

Figure 10. Low-energy isomers of Ag20.

ordered structures based on pentagonal patterns. Structure c in particular derives from a double 19-atom double icosahedron from which one apex and two atoms on the side have been subtracted. The pyramidal Td structure lies at higher energy (Table 5). Depending on the conditions of the calculations and in particular the basis set extension and the functionals, DFT calculations did not yield the same results for the structure of Ag20.27,32 Alternatively, Au20 was experimentally investigated via vibrational spectroscopy42 and the resulting data were shown to correlate with DFT vibrational frequencies originating from the pyramidal Td structure.115 Such finding is also obtained in the present DFTB framework as reported in Table 5 and Figure 11 (the labels a, b, and c correspond to the lowest energy structures). Moreover, from correlation of photoelectron spectroscopy data and DFT calculations, Au20− was also proven to have a Td symmetry.115 The present DFTB results are thus consistent with the literature. It is also interesting to see that although small differences in the structures are obtained for various charge situations, the isomer ordering is conserved in each case, constantly different in gold and silver, whatever the charge state. This shows that the addition or removal of a single electron at that size, N = 20, does not fully alter the structure of those clusters. We have also examined the larger clusters of 55 units, as well as their anions and cations. The results are shown in Table 6 and in Figure 12. In the present study, we started the DFTB optimizations of the neutral, cationic, and anionic clusters with the same structures tested by Häkkinen et al.116 in their LDA study of the Ag55− and Au55− clusters: the symmetrical cuboctahedral, decahedral, and icosahedral structures and three less

■

LARGER CLUSTERS AND CONVERGENCE TO THE BULK From the two previous sections, we showed that only DFTBχ for gold and DFTBγ for silver are able to reproduce reasonably both cluster and bulk properties. In the following, only these sets of parameters will be used. 8476


Article

The Journal of Physical Chemistry A Table 5. Energy Ordering (eV) of Various Isomers of Ag20 and Au20 at Different Charge States Obtained with the Modified DFTB Parameters (DFTBχ and DFTBγ for Gold and Silver, Respectively)a isomer (a) (b) (c) Td

ΔE 0 0.342 0.599 0.972

IP Ag20 5.522

EA

ΔE

0 0.234 0.487 1.297

(a) (b) (c) Td

0 0.245 0.452 1.043

EA

isomer

ΔE

Au20 1.837

1.180 1.766 1.803 0

Ag20+ (a) (b) (c) Td

Table 6. Energy Ordering (eV) of Various Isomers of Ag55 and Au55 at Different Charge States Obtained with the Modified DFTB Parameters (DFTBχ and DFTBγ for Gold and Silver, Respectively)a

6.421 Au20+

1.333

0.452 0.677 1.207 0 Ag20−

IP

EA

ΔE

Au20− 1.143 1.253 1.281 0

cubo deca G M S−C ico

2.605 1.855 1.042 1.795 2.412 0


2.735 1.777 1.243 1.936 2.239 0


2.742 2.300 1.067 2.045 2.810 0

5.024 Ag55+

3.057

1.620 1.313 0.169 0.006 0 1.204

Ionization potentials (IP) and electron affinities (EA) are reported in eV (all structures are relaxed).

EA

4.968

2.827

Au55+ 2.086 1.281 0.426 0.376 0 1.448

Ag55−

a

IP Au55

Ag55

Au55− 1.421 1.180 0.220 0 0.342 1.265

a

Ionization potentials (IP) and electron affinities (EA) are reported in eV (all structures are relaxed).

cuboctahedron, and decahedron have much higher energies than the less ordered structures (Table 6): the Sutton−Chen, Morse, and Glue structures are always at energies lower than 0.5 eV from the ground state structure, whereas the icosahedron, decahedron and cuboctahedron lie at energies higher than ≈1.2 eV. Let us mention that the exact ordering of these low energy isomers (G, M, S−C) is different from that of LDA calculations. Nevertheless, the present DFTB results are globally consistent with the Au55− LDA calculations of Häkkinen et al.116 who showed that, contrarily to other metallic 55-mers, the low-energy structures of Au55 are not the symmetric icosahedral structure, but less symmetric structures. The same authors also determined the theoretical photoelectron spectra corresponding to these various isomers and showed that these less ordered structures yielded spectra in good agreement with experiments. Convergence of Clusters Properties to the Bulk. Several investigations were concerned with convergence to the bulk and structural transitions in noble metal clusters,67−72,120,121 most of them using various semiempirical many-body potentials (MBP). Baletto et al. specifically addressed noble metal clusters with size ranging up to ∼2000 atoms, focusing on the transitions between non-fcc structures (icosahedron (Ih) or Marks decahedron (Dh)) toward fcc structures (truncated octahedron (TO)). They showed that an Ih to Dh transition occurs for N ≤ 400 in the case of silver and N < 100 for gold, and the transition from Ih to TO at N ≈ 400 for silver and

Benchmarking Density Functional Based Tight ... - ACS Publications

Recommend Documents