@Ab-Initio Simulations of Materials Using VASP

Ab-Initio Simulations of Materials Using VASP: Density-Functional Theory and Beyond JÜRGEN HAFNER Faculty of Physics and Center for Computational Materials Science, Universität Wien, Sensengasse 8, A–1090 Wien, Austria Received 22 February 2008; Accepted 9 May 2008 DOI 10.1002/jcc.21057 Published online 11 July 2008 in Wiley InterScience (www.interscience.wiley.com). Abstract: During the past decade, computer simulations based on a quantum-mechanical description of the interactions between electrons and between electrons and atomic nuclei have developed an increasingly important impact on solid-state physics and chemistry and on materials science—promoting not only a deeper understanding, but also the possibility to contribute significantly tomaterials design for future technologies. This development is basedon two important columns: (i) The improved description of electronic many-body effects within density-functional theory (DFT) and the upcoming post- DFTmethods. (ii) The implementation of the new functionals andmany-body techniqueswithin highly efficient, stable, and versatile computer codes, which allow to exploit the potential of modern computer architectures. In this review, I discuss the implementation of various DFT functionals [local-density approximation (LDA), generalized gradient approximation (GGA), meta-GGA, hybrid functional mixing DFT, and exact (Hartree-Fock) exchange] and post-DFT approaches [DFT + U for strong electronic correlations in narrow bands, many-body perturbation theory (GW) for quasiparticle spectra, dynamical correlation effects via the adiabatic-connection fluctuation-dissipation theorem (AC-FDT)] in the Vienna ab initio simulation package VASP. VASP is a plane-wave all-electron code using the projector-augmented wave method to describe the electron-core interaction. The code uses fast iterative techniques for the diagonalization of the DFT Hamiltonian and allows to perform total-energy calculations and structural optimizations for systems with thousands of atoms and ab initio molecular dynamics simulations for ensembles with a few hundred atoms extending over several tens of ps. Applications in many different areas (structure and phase stability, mechanical and dynamical properties, liquids, glasses and quasicrystals, magnetism and magnetic nanostructures, semiconductors and insulators, surfaces, interfaces and thin films, chemical reactions, and catalysis) are reviewed. © 2008 Wiley Periodicals, Inc. J Comput Chem 29: 2044–2078, 2008 Key words: density-functional theory; plane-wave basis; pseudopotentials; projector-augmented-waves; hybrid functionals; many-body perturbation theory; solid state physics; solid state chemistry; materials science; surface science; catalysis Introduction The last two decades have witnessed tremendous progress in the development of methods for ab initio calculations of materials prop- erties and for simulations of processes in materials. The cornerstone of this development was laid by density-functional theory (DFT), which casts the intractable complexity of the electron–electron inter- actions in many-electron systems into an effective one-electron potential, which is a functional of the electron density only.1–4 Although the formof this functionalwhichwouldmake the reformu- lation of the many-electron Schrödinger equation (the Kohn–Sham equations) exact is not known, starting with the pioneering work of Perdew, Becke, and coworkers a hierarchy of approximate func- tionals has been developed, which allow to predict many properties of solids with increasing accuracy.5–8 However, it cannot be over- looked that in the quest for better density-functionals, different cultures have developed in quantum chemistry and in condensed- matter physics. Quantum chemists have realized very early the need to correctly reproduce the atomic one- and two-electron densities in atoms and molecules, but in addition, to maximize accuracy a num- ber of empirical parameters are often introduced in a functional6,9–14 The resulting semiempirical approximations may be very accu- rate within their “training sets” (i.e., for systems similar to those in the data sets used for the optimization of the adjustable para- meters) but often fail when applied to different situations because the optimized parameters are not transferable from one system (e.g., small molecules) to another (e.g., solids or solid surfaces). A striking example is provided by the B3LYP functional,15, 16 which Correspondence to: J. Hafner; e-mail: juergen.hafner@univie.ac.at © 2008 Wiley Periodicals, Inc. Ab-Initio Simulations of Materials Using VASP 2045 is the most popular functional in quantum chemistry. This func- tional contains a total of eight empirical parameters, which allow to achieve a very high accuracy for almost all properties of small molecules. However, the B3LYP functional fails to reproduce the correct exchange-correlation energy of the homogeneous electron gas. Consequently, the accuracy of B3LYP predictions deteriorates rapidly with increasing molecular size,17,18 and the functional fails quite badly for metallic solids.19 The solid-state physics commu- nity, under the leadership of John Perdew, has followed another strategy: Exchange-correlation functionals may also be designed to satisfy as many exact constraints (known, e.g., from many-body theory) as possible. Parameters introduced into the functional are nonempirical in the sense that they are determined by these con- straints instead of being adjusted to reproduce experimental data. Nonempirical functionals that satisfy a sufficient number of con- straints20–23 aremore likely to be transferable and to perform equally well in different tasks. A particular challenge in the construction of nonempirical functionals is to respect the two paradigms prevalent in quantum chemistry and in solid state physics: the one- and two- electron densities in molecular quantum chemistry and the slowly varying electron densities in condensed matter physics. The rel- evance of the electron-gas limit to atoms, molecules, and solids and the difficulties to construct an exchange-correlation functional whose gradient-dependence allows to achieve good accuracy for atomization energies on one, and lattice parameters or surface ener- gies on the other hand have very recently been discussed by Perdew and coworkers.24,25 A substantial part of this review will be devoted to discuss the implementation of these functionals in a code using a plane-wave basis-set and their performance in describing variousmaterials prop- erties. However, even with the radical simplification introduced by DFT, ab initio calculations for solids were restricted, for a long time, to simple systems with small unit cells. The rapid development of electronic-structure theory of solids during the last decades was triggered by a seminal paper pub- lished in 1985 by Car and Parrinello,26 in which they proposed to solve the equations of motions of the coupled many-atom, many- electron system via a dynamical simulated annealing strategy. The Car-Parrinello (CP) method was designed to replace the traditional approach consisting of the iterative selfconsistent solution of the Kohn–Sham equations for the electrons, the calculation of the forces acting on the atoms via the Hellmann–Feynman theorem, and the integration of the Newtonian equations of motions of the ions—the procedure having to be repeated after each ionic integration step until the ground state of the many-atom, many-electron system had been reached. In addition, the CP paper introduced several other important innovations. One which is very important for the method- ology to be described below is the use of Fast Fourier Transforms to switch between real-space and momentum-space representations of the wave function because different parts of the calculation can be done most efficiently in one space or another: The kinetic energy has a diagonal representation in momentum space while the poten- tial energy is diagonal in real space. The second step forward was based on the observation that it is inefficient to do one part of the calculation (the solution of the Kohn-Sham equations) with very high accuracy while the other part (the determination of the equilibrium ionic configuration) is still far from convergence. This led immediately to the bold idea—already referred to above—that the total energy of a system could be minimized simultaneously with respect to both the electronic and ionic degrees of freedom. It was only after the publication of the CP paper that the full poten- tial of DFT has been exploited: what had been so far a technique used by a small community of solid-state theorists was transformed into a powerful tool for materials research, with applications in many different areas such as structural materials, catalysis and surface science, nanomaterials, biomaterials, and geophysics (for a recent review of the impact of DFT on materials research see ref. 27). Although much of the recent development was undoubtedly triggered by the Car-Parrinello paper, it is a bit ironical that the development of modern DFT calculations is characterized by a rather quick return to the more traditional approach. The reason is twofold: First, the CP approach of a dynamical updat- ing of the electronic degrees of freedom requires electrons and ions to be effectively decoupled such that, once the electronic ground-state has been reached, the system remains close to the adiabatic Born-Oppenheimer surface. This condition is met with good accuracy for insulators and wide-gap semiconductors but violated for metals and narrow-gap materials. The second rea- son is that the minimization of the total energy does not allow an efficient control of charge-density fluctuations during the itera- tive process—for metallic systems such fluctuations (often referred to as “charge-sloshing”) may even prevent a convergence of this process. Modern DFT calculations for solids are determined by sev- eral technical choices: (i) The choice of a basis set to expand the Kohn-Sham eigenfunctions. Essentially, the choice is between plane waves and localized basis functions. (ii) The interactions between the ionic core and the valence electrons can be described either by a full-potential approach or by a pseudopotential elim- inating the need to account for the complex nodal character of the valence orbitals. (iii) The method adopted for the determina- tion of the eigenstates of the Kohn-Sham Hamiltonian. (iv) The description of the electron-electron interactions by choosing an exchange-correlation functional within the hierarchy of functionals proposed within DFT.3 In those cases where DFT alone does not provide an adequate solution (strong electronic correlations, excited eigenstates, . . .) post-DFT corrections such as many-body perturba- tion theory28,29 or dynamical mean field theory30,31 may be used to improve the DFT predictions. In the present review, I shall concentrate on the theoretical back- ground of the Vienna ab initio simulation package VASP developed by Georg Kresse and his coworkers32–35 and on applications of this code in key areas ofmodern solid-state physics and chemistry.VASP is a plane-wave code for ab-initio density-functional calculations. It attempts to match the accuracy of the most advanced all-electron codes by using a projector-augmented-wave approach (PAW)36 for describing the electron-ion interaction. A stable and accurate solu- tion of the Kohn-Sham equations, as well as a favorable scaling of the computational effort with system size, are achieved by adopting iterative diagonalization techniques and optimized charge-mixing routines. Different levels of exchange-correlation functionals and different post-DFT approaches have been implemented. A vari- ety of routines added to the basic DFT solver allow to calculate a wide variety of materials properties. I begin by reviewing the basic methodology. Journal of Computational Chemistry DOI 10.1002/jcc 2046 Hafner • Vol. 29, No. 13 • Journal of Computational Chemistry Figure 1. Convergence of the relative bond-lengths of Cl2, ClF, and HCl molecules calculated using various local basis sets [aug-cc-pVXZ with X = 2 (double), 3 (triple), 4 (quadruple), 5 (quintuple)] rela- tive to the plane-wave results. Cf. text. After Kresse et al.37 [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.] The Vienna Ab Initio Simulation Package VASP — Basic Methodology Why Plane Waves? Modern electronic structure methods fall into two broad classes, depending on the choice of the basis set for the expansion of the valence orbitals, charge densities and potentials: plane-wave meth- ods or methods using some kind of localized basis functions, e.g., Gaussian-type orbitals. The use of a plane-wave basis has several immediate advantages: (i) It is easy to change from a real-space representation (where the potential energy V has a diagonal repre- sentation) via a Fast Fourier Transform to momentum-space where the kinetic energy T is diagonal. (ii) The control of basis-set con- vergence is almost trivial; it is sufficient to monitor the eigenvalues and total energies as a function of the cut-off energy, i.e., the high- est kinetic energy of a plane-wave within the chosen basis set. (iii) The Hellmann-Feynman forces acting on the atoms and the stresses on the unit cell may be calculated straightforwardly in terms of the expectation value of the Hamiltonian with respect to the ionic coordinates. (iv) Basis-set superposition errors that have to be carefully controlled in calculations based on local basis sets are avoided. On the other hand, a set of local Gaussian basis functions allows an analytic integration of the 1/r singularity of the Coulomb potential—this is instrumental to a fast calculation of exact (Hartree- Fock) exchange. The treatment of exact exchange is more difficult with a plane-wave basis. We shall return to this point below. A very important point is that a reasonable convergence of a plane-wave expansion can be achieved only if the nodal character of the valence orbitals is eliminated, i.e., if the ion–electron inter- action is described by some kind of pseudopotential. This raises the question of the accuracy and transferability of pseudopotentials, and the necessity linearize the valence-core exchange-correlation inter- actions. I shall discuss below how these problems can be solved within the projector-augmented wave method.35,36 The pseudopotential and related methods are sometimes regarded as unnecessary approximations by quantum-chemists, while the plane-wave community considers local basis set results with some suspicion because of basis-set completness and basis- set superposition errors. Therefore, a demonstration that both approaches lead to perfectly converged results is very important. Kresse et al.37 have reported benchmark results on optimized geometries and atomization energies of molecules calculated with VASP (using a plane-wave basis set and the projector-augmented wave method) and GAUSSIAN03 (G03)38 using large local basis sets. An illustration of the results is given in Figure 1 for the diffi- cult (for plane waves) case of diatomic molecules containing Cl. It is evident that a very large local basis set [an augmented correlation- consistent polarized valence quintuple-zeta (aug-cc-pV5Z) basis set] is required tomatch the converged plane-wave results.Generally bond lengths of small molecules agree within 0.1%. The necessity to use a large basis set also has a consequence on the computational cost. While for the small molecules the plane-wave VASP and local basis-set quadruple-zeta calculations perform roughly equally, the cost increases dramatically for the quintuple-zeta basis necessary to achieve full convergence. In fairness, it must be pointed out that the G03 calculations are all-electron calculations, while the VASP- PAW calculations are full-potential valence-only calculations. As a consequence, G03 calculations will be very expensive for larger systems for which plane-wave calculations are suited very well. The influence of the frozen-core approximation used in the PAWmethod will be discussed below. Very recently, a comparative investigation of the performance of plane-wave (VASP) and local-basis set methods (using the GAUS- SIAN and SIESTA packages40) in structural studies of small gold clusters has been presented by Gruber et al.39 For a wide class of relatively compact cluster structures the authors found excel- lent agreement between between the binding energies calculated using both methods, while planar structures where found to have a somewhat reduced stability in the local-basis set calculations. This differencewas attributed to the fact that the quality of the plane-wave basis set is independent of the topology of the systemwhile the qual- ity of a basis composed of atom-centred local orbitals depends on the relative atomic positions (a situation which is evidently reminiscent of the basis-set superposition error). It was concluded that the rela- tively lower binding energy of planar clusters provided by SIESTA and GAUSSIAN03 could be a consequence of a lower “effective quality” of the basis set for systems that are more extended in one or two dimensions compared with more compact structures Potentials, pseudopotentials Pseudopotentials have been introduced to avoid the need for an explicit treatment of the strongly bound and chemically inert core electrons. They are a necessary ingredient of all plane-wave meth- ods, but they can also be used in local-basis set methods to reduce the computational effort. The theory of pseudopotentials is mature, but the practice of constructing accurate, transferable, and effi- cient pseudopotentials is far from straightforward. Methods for generating pseudopotentials include the ‘norm-conserving’ pseu- dopotentials41,42 (the “norm-conservation” criterion applied to the node-less pseudo wave functions ensures that not only the loga- rithmic derivative of the exact and pseudo-wavefunctions, but also Journal of Computational Chemistry DOI 10.1002/jcc Ab-Initio Simulations of Materials Using VASP 2047 their derivatives with respect to the energy agree at the chosen reference energy and cut-off radius) and the “ultrasoft” pseudopo- tentials (where the norm-conservation criterion is dropped, but the logarithmic derivatives are matched at two or more reference energies spanning the entire range of eigenvalues of the valence electrons).41,43 The ultrasoft pseudopotentials have the merit to make calculations for first-row elements and for systems with d− or f−electrons feasible at tractable effort. The criterion for the qual- ity of a pseudopotential is not how well it matches experiment, but how well it reproduces the results of accurate all-electron calcu- lations. A certain drawback of pseudopotential calculations is that because of the nonlinearity of the exchange interaction between valence and core electrons, elaborate non-linear core corrections44 are required for all systems where the overlap between valence- and core-electron densities is not completely negligible. This deficiency may be removed by using the projector-augmented wave method. Projector-augmented waves The projector-augmented wave (PAW) method originally intro- duced by Blöchl36 represents an attempt to achieve simultaneously the computational efficiency of the pseudopotential method as well as the accuracy of the full-potential linearized augmented-plane- wave (FLAPW) method,45 which is commonly regarded as the benchmark for DFT calculations on solids. Unlike the pseudopo- tential approach, the PAW method accounts for the nodal features of the valence orbitals and ensures orthogonality between valence and corewave functions. In the PAWapproach, the all-electron (AE) valence wave functionsψAEn are reconstructed from