Article pubs.acs.org/JCTC
Terms of Use
Bulk Properties of Transition Metals: A Challenge for the Design of Universal Density Functionals Patanachai Janthon,†,‡,§ Sijie (Andy) Luo,∥ Sergey M. Kozlov,† Francesc Viñes,† Jumras Limtrakul,‡,§ Donald G. Truhlar,∥ and Francesc Illas*,† †
Departament de Química Física & Institut de Química Teòrica i Computacional (IQTCUB), Universitat de Barcelona, c/Martí i Franquès 1, 08028, Barcelona, Spain ‡ Department of Chemistry and NANOTEC Center for Nanoscale Materials Design for Green Nanotechnology, Kasetsart University, Bangkok 10900, Thailand § PTT Group Frontier Research Center, PTT Public Company Limited, 555 Vibhavadi Rangsit Road, Chatuchak, Bangkok 10900, Thailand ∥ Department of Chemistry, Chemical Theory Center, and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *
ABSTRACT: Systematic evaluation of the accuracy of exchangecorrelation functionals is essential to guide scientists in their choice of an optimal method for a given problem when using density functional theory. In this work, accuracy of one Generalized Gradient Approximation (GGA) functional, three meta-GGA functionals, one Nonseparable Gradient Approximation (NGA) functional, one metaNGA, and three hybrid GGA functionals was evaluated for calculations of the closest interatomic distances, cohesive energies, and bulk moduli of all 3d, 4d, and 5d bulk transition metals that have face centered cubic (fcc), hexagonal closed packed (hcp), or body centered cubic (bcc) structures (a total of 27 cases). Our results show that including the extra elements of kinetic energy density and Hartree−Fock exchange energy density into gradient approximation density functionals does not usually improve them. Nevertheless, the accuracies of the Tao−Perdew−Staroverov−Scuseria (TPSS) and M06-L meta-GGAs and the MN12-L metaNGA approach the accuracy of the Perdew−Burke−Ernzerhof (PBE) GGA, so usage of these functionals may be advisable for systems containing both solid-state transition metals and molecular species. The N12 NGA functional is also shown to be almost as accurate as PBE for bulk transition metals, and thus it could be a good choice for studies of catalysis given its proven good performance for molecular species. Approximation3,4 (GGA) or Nonseparable Gradient Approximation5 (NGA), which constitute the second rung; further addition of kinetic energy density leads to meta functionals, such as meta-GGAs and meta-NGAs, which form the third rung, and, finally, including Hartree−Fock (HF) exchange leads to hybrid functionals, in particular hybrid gradient approximations and hybrid meta approximations, which form the fourth rung. Through the third rung, the exchange-correlation energy density depends only on local variables, but hybrid functionals are nonlocal. One tentative conclusion about xc functionals that has been advanced is that adding elements from the third and fourth rungs has not (with some exceptions) led to better performance for molecules containing metals,6,7 and various studies have been carried out with gradient approximations that provide further experience related to this issue.7−9 Furthermore, although hybrid functionals are justifiably the most popular
1. INTRODUCTION Kohn−Sham Density Functional Theory (KS-DFT) has become a workhorse for treating complex problems in both the gaseous and condensed phases. The accuracy of KS-DFT rests entirely on the accuracy of one’s approximation to the exchange-correlation (xc) functional. Functional development has considered a wide range of application targets. In particular, there are numerous validation studies for atoms and molecules, many more than for condensed-phase quantum chemistry and solid-state physics. A particularly important area that still needs further study in its own right, in addition to its importance for studying heterogeneous catalysis and electrochemistry,1 is the accuracy of available xc functionals for transition metal solids. One route to improving xc functionals is to add more elements, which takes one to a higher rung on the Jacob’s ladder2 of xc functionals. The first rung is the Local Spin Density Approximation (LSDA), where the xc functional depends on only local spin densities. Including dependences not only on spin densities but also on their gradients yields gradient approximations, such as the Generalized Gradient © 2014 American Chemical Society
Received: June 20, 2014 Published: August 6, 2014 3832
dx.doi.org/10.1021/ct500532v | J. Chem. Theory Comput. 2014, 10, 3832−3839
Journal of Chemical Theory and Computation
Article
kind of functionals in molecular chemistry10 (because of their good performance), for extended systems they have increased computational demands as compared to local functionals due to the long range of the exchange in real space programs and to the requirement for dense Brillouin zone sampling in plane wave programs.11 Hybrid functionals also have the disadvantage of bringing in HF static correlation error, which is important for many molecular and solid-state transition metal compounds.6−9,12,13 In fact, when one tries to estimate the optimal percentage of HF exchange in a hybrid functional suitable for transition metal solids, the resulting fraction is close to zero,14,15 leading to the conclusion that the inclusion of HF exchange is detrimental for the accuracy of xc functionals for solid transition metals. The problems with Hartree−Fock exchange for metals are in fact well-known in that Hartree− Fock exchange for free electrons (which constitute the zerothorder model of a metal) leads to an unphysical density of states such that there is a singularity in the electron group velocity at the Fermi level, although this defect is eliminated if the Hartree−Fock exchange is screened by correlation effects.16 As a result of the above kinds of issues with including Hartree−Fock exchangecost, performance, and an incorrect treatment of the density of stateslocal functionals have remained the most popular choice for solid-state transition metal investigations even though hybrid functionals by construction reduce self-interaction error. A similar situation pertains even when the choice of xc functionals is restricted to local ones, where, although meta functional also has the potential to reduce self-interaction error,17 gradient approximations have been preferred to meta functionals partly for the greater simplicity of popular GGAs, but also because meta functionals have not always seemed to significantly improve the accuracy. However, this situation seems unsatisfactory because both hybrid functionals and meta functionals have better performance than GGAs on molecules composed of atoms lighter than the transition metals.18−21 Thus, there is a strong interest in certain fields, such as heterogeneous catalysis, nanotechnology, and materials chemistry, to use meta or hybrid functionals because these research lines often deal with the interaction between transition metal surfaces and light maingroup compounds. Meta and hybrid functionals also yield more accurate band gaps,15,22,25 which may be important for metaloxide systems that are ubiquitous as supports in catalysis, photocatalysis, and nanotechnology. In light of these practical considerations, there is a need for a more systematic investigation of the accuracy of meta and hybrid functionals for various properties of transition metal solids. In this article, we assess the accuracy of several meta and hybrid functionals for bulk transition metals. Namely, we consider shortest interatomic distances, δ, cohesive energies, Ecoh, and bulk moduli, B0, of all 27 bulk transition metals that have face centered cubic (fcc), hexagonal closed packed (hcp), or body centered cubic (bcc) structures. We present new results for the following functionals: GGA: SOGGA11;26 NGA: N12; 5 meta-GGAs: Tao−Perdew−Staroverov−Scuseria (TPSS), 27 revised Tao−Perdew−Staroverov−Scuseria (revTPSS),28 and meta M06-L;29 meta-NGA: MN12-L;30 hybrid gradient approximations with unscreened Hartree− Fock exchange: Perdew−Burke−Ernzerhof (PBE0),31,32 and Lee−Yang−Parr B3LYP;32 and hybrid gradient approximation with screened Hartree−Fock exchange: Heyd−Scuseria− Ernzerhof (HSE06).33
We compare the new results obtained here to the results obtained previously34 with PBE, a GGA functional that performed best in a previous study of two LSDA functionals [the Perdew−Zunger functional (PZ81)35 based on Ceperley− Alder (CA) numerical results36 and the Vosko−Wilk−Nusair functional (VWN, usually denoted also as VWN5)]37 and four gradient approximations [PBE,38 PBEsol,39 RPBE,40 and PW9141]. Note that the performance of the PW91 GGA was found to be close to that of PBE.
2. COMPUTATIONAL DETAILS The computational procedures were kept as close as possible to those used in a previous work by Janthon et al.34 Calculations were performed using a locally modified version of the Vienna Ab Initio Simulation Package (VASP).42 The atomic cores were described by the Projector Augmented Wave (PAW) method,43 using the potentials recommended in the documentation (see Table S1 in the Supporting Information for full specification) which take account of kinetic energy density dependencies when present. An optimized Monkhorst−Pack44 k-points grid of 7 × 7 × 7 was used in all bulk calculations, which was found to be sufficient for accurate total energy calculations with the smallest unit cell. A kinetic energy cutoff of 415 eV for the plane-wave basis set was employed throughout, guaranteeing variations of total energy below 1 meV with respect to more complete basis sets. Scalar relativistic effects in the core region are included in the PAW potentials. Previous calculations have shown that relativistic effects on the valence electrons of heavy transition metals lead to negligible deviations from this standard PAW approximation, i.e., interatomic distance changes in the 0.002− 0.005 Å range and differences in bulk moduli of 3−5 GPa.45 The electronic structure calculations were not spin polarized, with the exception of the calculations on the ferromagnetic Fe, Ni, and Co bulk systems and on isolated metal atoms. Although it is a digression, we note in passing that some of the calculated magnetic moments agree nicely with experimental values. For instance, the PBE values of 2.21, 1.55, and 0.59 μB agree well with experimental values of 2.2, 1.7, and 0.6 μB for Fe, Co, and Ni, respectively.46 Values predicted by TPSS are very close to the PBE ones with differences smaller than 0.02 μB, whereas HSE06 overestimates the magnetic moments with predicted values of 2.89, 1.85, and 0.82 μB for Fe, Co, and Ni, respectively. Note also that, even though a nearly degenerate antiferromagnetic solution for the bulk structure of Cr with a magnetic moment of of 0.92 μB has been reported from calculations at the PW91 level,47 the present calculations for Cr correspond to a nonmagnetic solution. Optimizations were performed using the tetrahedron smearing method of Blöchl et al.48 with an energy width of 0.2 eV to speed up convergence with the final energies extrapolated to zero smearing. In bulk calculations, ionic positions and cell volumes were optimized using the conjugate gradient algorithm until pressures and total energies were converged within 0.01 GPa and 10 meV, respectively. Calculated cohesive energies are sensitive to the treatment adopted for calculating the energies of the atoms. The isolated atoms were placed in a large unit cell with a broken symmetry of 9 × 10 × 11 Å dimensions to ensure proper occupancy of orbitals. When needed, orbital occupancy was imposed to match that of the isolated atom; further details concerning atomic configuration and energies are found in a previous article.34 The large unit cell allows one to sufficiently suppress 3833
dx.doi.org/10.1021/ct500532v | J. Chem. Theory Comput. 2014, 10, 3832−3839
Journal of Chemical Theory and Computation
Article
Table 1. Calculated and Experimental Shortest Interatomic Distances (δ, in Å) for Bulk Transition Metals GGA Sc Ti V Cr Fe Co Ni Cu Zn Y Zr Nb Mo Tc Ru Rh Pd Ag Cd Hf Ta W Re Os Ir Pt Au a
NGA
meta-GGA
meta-NGA
hybrid
CS
PBE
SOGGA11
N12
TPSS
revTPSS
M06-L
MN12-L
PBE0
HSE06
B3LYP
exptl.
exptl. corr.a
hcp hcp bcc bcc bcc hcp fcc fcc hcp hcp hcp bcc bcc hcp hcp fcc fcc fcc hcp hcp bcc bcc hcp hcp fcc fcc fcc
3.214 2.884 2.577 2.459 2.453 2.470 2.489 2.567 2.644 3.543 3.197 2.873 2.749 2.722 2.658 2.717 2.794 2.941 2.990 3.139 2.875 2.751 2.755 2.694 2.741 2.811 2.950
3.066 2.806 2.561 2.479 2.498 2.496 2.506 2.502 2.558 3.419 3.144 2.882 2.716 2.707 2.643 2.695 2.756 2.861 2.941 3.067 2.875 2.761 2.756 2.700 2.747 2.809 2.921
3.194 2.864 2.582 2.467 2.457 2.438 2.455 2.551 2.658 3.511 3.166 2.871 2.692 2.685 2.627 2.683 2.766 2.901 3.018 3.110 2.854 2.746 2.723 2.673 2.712 2.780 2.912
3.198 2.862 2.559 2.442 2.431 2.444 2.453 2.517 2.522 3.522 3.179 2.862 2.734 2.700 2.637 2.689 2.761 2.890 2.920 3.114 2.857 2.736 2.742 2.685 2.727 2.789 2.912
3.192 2.855 2.554 2.435 2.424 2.438 2.443 2.499 2.465 3.518 3.171 2.852 2.729 2.689 2.624 2.674 2.742 2.867 2.880 3.100 2.844 2.724 2.730 2.673 2.710 2.769 2.888
3.149 2.859 2.581 2.457 2.484 2.472 2.419 2.476 2.490 3.562 3.218 2.888 2.748 2.706 2.649 2.702 2.794 2.940 2.883 3.151 2.886 2.744 2.742 2.688 2.733 2.806 2.945
3.013 2.850 2.559 2.428 2.598 2.489 2.370 2.485 2.552 3.615 3.195 2.871 2.726 2.689 2.623 2.676 2.755 2.899 3.002 3.099 2.839 2.718 2.724 2.656 2.705 2.771 2.913
3.238 2.862 2.548 2.429 2.502 2.452 2.471 2.563 2.638 3.554 3.187 2.860 2.724 2.679 2.608 2.674 2.768 2.924 2.964 3.123 2.864 2.725 2.721 2.660 2.712 2.778 2.922
3.252 2.865 2.548 2.428 2.522 2.482 2.480 2.564 2.640 3.559 3.185 2.857 2.721 2.677 2.610 2.672 2.765 2.925 2.971 3.151 2.861 2.724 2.720 2.661 2.711 2.777 2.919
3.230 2.883 2.573 2.457 2.517 2.498 2.510 2.606 2.781 3.565 3.213 2.888 2.755 2.717 2.649 2.723 2.827 2.992 3.139 3.148 2.886 2.761 2.759 2.698 2.759 2.835 2.992
3.254 2.897 2.613 2.498 2.460 2.497 2.493 2.556 2.665 3.556 3.179 2.859 2.725 2.710 2.650 2.539 2.753 2.889 2.979 3.131 2.860 2.741 2.567 2.675 2.715 2.772 2.879
3.244 2.889 2.606 2.485 2.450 2.488 2.484 2.544 2.645 3.548 3.174 2.854 2.721 2.705 2.642 2.532 2.745 2.877 2.959 3.126 2.856 2.738 2.562 2.671 2.710 2.766 2.870
ZPE and finite temperature corrections to experimental values adapted from Lejaeghere et al.52
the volumes at the equilibrium geometry and at geometries with ±0.05 and ±0.10 Å variations of the lattice constants. A more accurate procedure would be the adjustment of energy points to a Murnaghan equation of state.49 Bulk moduli obtained at the PBE level with a similar calculation setup but through the adjustment to the Murnaghan equation are ∼5 GPa lower than those calculated here.11 In the present work, we consider the 27 transition metals that exhibit close-packed structures, which may display fcc, hcp, or that exhibit a bcc structure. We specifically exclude La and solid Hg, which have hexagonal and rhombohedral unit cells, respectively, and Mn, which has a cubic unit cell containing 58 atoms, although it is similar to a bcc structure. The reason for these exclusions is that preliminary work showed that metals exhibiting these structures can exhibit different trends than the typical metals, and specialized consideration is required to do them justice. In the following, as in the previous work,34 in addition to cohesive energy and bulk modulus, the shortest interatomic distance within a crystal cell, δ, will be compared to the experimental values. We present results for individual metals, and to more clearly see some trends, we also present mean errors, in particular mean signed error (MSE), mean absolute error (MAE), and mean absolute percentage error (MAPE).
spurious interactions between periodic images (
