Calculations Manual

From Wiki
Jump to: navigation, search

Overview

This manual explains our calculation methods, adjustments made to our calculated energies, and errors in our calculations.

Total Energy Calculation Details

We use density functional theory as implemented in the Vienna Ab Initio Simulation Package (VASP) software[1] to evaluate the total energy of compounds. For the exchange-correlational functional, we employ a mix of Generalized Gradient Approximation (GGA) and GGA+U, as described later in this manual. We use the Projector Augmented Wave (PAW) method for modeling core electrons with an energy cutoff of 520 eV. This cutoff corresponds to 1.3 times the highest cutoff recommended among all the pseudopotentials we use (more details can be found in the Pseudopotentials Choice manual). All calculations are performed at 0K and 0atm. All computations are performed with spin polarization on and with magnetic ions in a high-spin ferromagnetic initialization (the system can of course relax to a low spin state during the DFT relaxation). We used a k-point mesh of 1000/(number of atoms in the cell) for all computations, the Monkhorst-Pack method for the k-point choices (but \Gamma-centered for hexagonal cells) and the tetrahedron method to perform the k-point integration. Pymatgen could change those default parameters if they are not adequate for the computation (e.g., switch to another k-point integration scheme). Some details of our calculation method can be found in ref [2]; however, the Materials Project has updated many parameters as documented throughout the Calculations Wiki.

Crystal structures

We use input structures from the Inorganic Crystal Structure Database (ICSD),[3] and relax all cell and atomic positions in our calculation two times in consecutive runs. When multiple crystal structures are present for a single chemical composition, we attempt to evaluate all unique structures as determined by an affine mapping technique.[4]

Total energy convergence

We currently employ a k-point mesh of 1000/number of atoms. We performed a convergence test of total energy with respect to k point density and convergence energy difference for a subset of chemically diverse compounds for a previous parameter set, which employed a smaller k-point mesh of 500/atom. Using 500/atom k-point mesh, the numerical convergence for most compounds tested was within 5 meV/atom, and 96% of compounds tested were converged to within 15 meV/atom. Results for the new parameter set will be better due to the higher k-point mesh employed. Convergence will depend on chemical system; for example, oxides were generally converged to less than 1 meV/atom.[2]

Structure convergence

The energy difference for ionic convergence is set to 0.0005 * natoms in the cell. Data on expected accuracy on cell volumes can be found in a previous paper.[2] We have found these parameters to yield well-converged structures in most instances; however if the structures are to be used for further calculations that require strictly converged atomic positions and cell parameters (e.g. elastic constants, phonon modes, etc.) we recommend that users re-optimize the structures with tighter cutoffs or in force convergence mode.

Total Energy Adjustments

To better model energies across diverse chemical spaces, we apply several adjustments to the total energy. These adjustments are described below.

Default corrections

Gases, liquids, and elements

Our total energy calculations are for periodic solids at 0K. For elements that are gaseous in their standard state, our raw calculations do not represent the same phase as standard experimental data. In addition, there may be GGA errors associated with reducing gas phases in gas-to-solid reactions. Rather than calculating the liquid/gas energies directly, we adjust the energies of several elements that are liquid or gaseous at room temperature using Wang's method.[5] The best-tested fit is for oxygen gas reacting to form to oxides. We have adjusted energies of the following compounds:

  • O2, N2, Cl2, F2, H2

The adjusted energies outside of O2 are not as well-tested and calculations involving these elements should be taken with greater caution. For example, it is assumed in our elemental gas adjustments that reactions to the solid will involve reduction of the gas.

GGA+U calculations and adjustments to match GGA energies

Some compounds are better modeled with a U correction term to the density functional theory Hamiltonian while others are better modeled without (i.e., straight GGA). Energies from calculations with the +U correction are not directly comparable to those without. To obtain better accuracy across chemical systems, we use GGA when appropriate, GGA+U otherwise, and mix energies from the two calculation methodologies by adding an energy correction term to the GGA+U calculations to make them comparable to the GGA calculations. The idea behind this approach is to split reactions into sub-reactions that are well-modeled by GGA, well-modeled by GGA+U, or a binary formation reaction that can be estimated from known experimental data. More details on this method can be found in ref [6].

Band Structure and DOS computations details

We used the relaxed structure from the total energy run to perform using the same GGA functional (and U if any) a DOS and band structure run. This data is available for many entries (a few thousands as for now) in the materials explorer. We ran first a static run with a uniform (monkhorst pack or gamma centered for hexagonal systems) k-point grid. From this run the DOS and the charge density was extracted. The charge density was used to perform then a band structure computation for k-point along the symmetry lines of the Brillouin zone. The symmetry lines have been determined following the methodology from Curtarolo et al.[7] using the aconvasp software.

The band structure data is displayed with full lines for up spin and dashed lines for down spin. For insulators, the band gap is computed according to the band structure and the nature of the gap (direct or undirect) as well as the k-points involved in the band gap transition are displayed. The VBM and CBMs are displayed for insulators as well by red and white dots. The web site does not display the 10% highest bands as those higher energies bands tend to converge more slowly.

The DOS indicates the full DOS by default but projections along atoms or orbitals are available too through a scroll down menu. Please note that the DOS data and band structure can be slightly different as the kpoint grid is not the same. For instance, the k-point grid for the DOS might not include Gamma while the band structure does.

Accuracy of Total Energies

To estimate the accuracy of our total energy calculations, we compute reaction data and compare against experimental data. Note that this data set was compiled using a lower k-point mesh and pseudopotentials with fewer electrons than the current Materials Project parameter set.

Estimating errors in calculated reaction energies

The accuracy of calculated reaction energies depends on the chemical system investigated. In general, GGA calculations have similar errors among chemically similar systems. Hence, reaction energies between chemically similar systems (e.g., a reaction where the reactants and products are all oxides, such as MgO + Al2O3 -> MgAl2O4) tend to have smaller errors than reactions between chemically dissimilar systems (e.g., between metals and insulators).

Figure 1: Errors in Calculated Formation Energies for 413 binaries in the Kubaschewski Tables. Energies are normalized to per mol atom.

To provide a quantitative indicator of the error we may expect from the reaction calculator, we have computed the reaction energies of File:KtablesBinaries.txt in the Kubaschewski Tables formed with Group V, VI and VII anions. Figure 1 shows the errors in the calculated formation energies (compared to the experimental values) for these compounds. The mean absolute error (MAE) is around 14 kJ mol-1. 75% of the calculated formation energies are within 20 kJ mol-1. We also found that compounds of certain elements tend to have larger errors. For example, Bi, Co, Pb, Eu, U, Tl and W compounds often have errors larger than 20 kJ mol-1.

It should be noted that while an MAE of 14 kJ mol-1 is significantly higher than the desired chemical accuracy of 4 kJ mol-1, it compares fairly well with the performance of most quantum chemistry calculations[8] Other than the most computationally expensive model chemistries such as G1-G3 and CBS, the reaction energy errors of most computational chemistry model chemistries are well above 10 kJ mol-1.

For oxidation of the elements into binary compounds, an average error of ~4% or 33 kJ/mol-O2 is typical.[9] For conventional ternary oxide formation from the elements, we have found a mean relative absolute error of about 2%.[6]

Sources of error

The largest contribution to the error comes from the inability of the GGA to fully describe electronic exchange and correlation effects. In addition, there is some error associated with neglecting zero-point effects and with comparing 0K, 0atm computations with room-temperature enthalpy experiments. The latter effect was estimated to contribute less than 0.03 eV/atom by Lany.[10] The stability of antiferromagnetic compounds may be underestimated, as the majority of our calculations are performed ferromagnetically only. The effect of magnetism may be small (under 10 meV/atom) or large (100 meV/atom or greater), depending on the compound. For compounds with heavy elements, relativistic effects may lead to greater-than-expected errors.

GGA errors on reaction energies between chemically similar compounds

We recently conducted a more in-depth study comparing GGA (+U) reaction energies of ternary oxides from binary oxides on 135 compounds. [11]

The main conclusions are:

  • The error in reaction energies for the binary oxide to ternary oxides reaction energies are an order of magnitude lower than for the more often reported formation energies from the element. An error intrinsic to GGA (+U) is estimated to follow a normal distribution centered in zero (no systematic underestimation or overestimation) and with a standard deviation around 24 meV/at.
  • When looking at phase stability (and for instance assessing if a phase is stable or not), the relevant reaction energies are most of the time not the formation energies from the elements but reaction energies from chemically similar compounds (e.g., two oxides forming a third oxide). Large cancelation of errors explain this observation.
  • The +U is necessary for accurate description of the energetics evene when reactions do not involve change in formal oxidation states

Accuracy of Calculated Volumes

A discussion of errors in calculated volumes can be found in the Volume Change Error manual.

Accuracy of Band Structures

Note: the term 'band gap' in this section generally refers to the fundamental gap, not the optical gap. The difference between these quantities is reported to be small in semiconductors but significant in insulators. [12]

Figure 2: Experimental versus computed band gaps for 237 compounds in an internal test. The computed gaps are underestimated by an average factor of 1.6, and the residual error even after accounting for this shift is significant (MAE of 0.6 eV). We thank M. Chan for her assistance in compiling this data.

Density functional theory is formulated to calculate ground state properties. Although the band structure involves excitations of electrons to unoccupied states, the Kohn-Sham energies used in solving the DFT equations are often interpreted to correspond to electron energy levels in the solid.

The correspondence between the Kohm-Sham eigenvalues computed by DFT and true electron energies is theoretically valid only for the highest occupied electron state. The Kohn-Sham energy of this state matches the first ionization energy of the material, given an exact exchange-correlation functional. However, for other energies, there is no guarantee that Kohn-Sham eigenvalues will correspond to physical observables.

Despite the lack of a rigorous theoretical basis, the DFT band structure does provide useful information. In general, band dispersions predicted by DFT are reported to match experimental observations; one small test of band dispersion accuracy found that errors ranged from 0.1 to about 0.4 eV.[13] However, predicted band gaps are usually severely underestimated. Therefore, a common way to interpret DFT band structures is to apply a 'scissor' operation whereby the conduction bands are shifted in energy by a constant amount so that the band gap matches known experimental observations.

Band gaps

In general, band gaps computed with common exchange-correlation functionals such as the LDA and GGA are severely underestimated. Typically the disagreement is reported to be around 50% in the literature. Some internal testing by the Materials Project supports these statements; typically, we find that band gaps are underestimated by about 40% (Figure 2). We additionally find that several known insulators are predicted to be metallic.

Origin of band gap error and improving accuracy

The errors in DFT band gaps obtained from calculations can be attributed to two sources: 1. Approximations employed to the exchange correlation functional 2. A derivative discontinuity term, originating from the true density functional being discontinuous with the total number of electrons in the system.

Of these contributions, (2) is generally regarded to be the larger and more important contribution to the error. It can be partly addressed by a variety of techniques such as the GW approximation but typically at high computational cost.

Strategies to improve band gap prediction at moderate to low computational cost now been developed by several groups, including Chan and Ceder (delta-sol),[14] Heyd et al. (hybrid functionals) [15], and Setyawan et al. (empirical fits) [16]. (These references also contain additional data regarding the accuracy of DFT band gaps.) The Materials Project may employ such methods in the future in order to more quantitatively predict band gaps. For the moment, computed band gaps should be interpreted with caution.

Citation

To cite the calculation methodology, please reference the following works:

  1. A. Jain, G. Hautier, C. Moore, S.P. Ong, C.C. Fischer, T. Mueller, K.A. Persson, G. Ceder., A High-Throughput Infrastructure for Density Functional Theory Calculations, Computational Materials Science, vol. 50, 2011, pp. 2295-2310.
  2. A. Jain, G. Hautier, S.P. Ong, C. Moore, C.C. Fischer, K.A. Persson, G. Ceder, Accurate Formation Enthalpies by Mixing GGA and GGA+U calculations, Physical Review B, vol. 84, 2011, p. 045115.

References

  1. Kresse, G. & Furthmuller, J., 1996. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical Review B, 54, pp.11169-11186.
  2. 2.0 2.1 2.2 A. Jain, G. Hautier, C. Moore, S.P. Ong, C.C. Fischer, T. Mueller, K.A. Persson, G. Ceder., A High-Throughput Infrastructure for Density Functional Theory Calculations, Computational Materials Science. vol. 50 (2011) 2295-2310.
  3. G. Bergerhoff, The inorganic crystal-structure data-base, Journal Of Chemical Information and Computer Sciences. 23 (1983) 66-69.
  4. R. Hundt, J.C. Schön, M. Jansen, CMPZ - an algorithm for the efficient comparison of periodic structures, Journal Of Applied Crystallography. 39 (2006) 6-16.
  5. L. Wang, T. Maxisch, G. Ceder, Oxidation energies of transition metal oxides within the GGA+U framework, Physical Review B. 73 (2006) 1-6.
  6. 6.0 6.1 A. Jain, G. Hautier, S.P. Ong, C. Moore, C.C. Fischer, K.A. Persson, G. Ceder, Formation Enthalpies by Mixing GGA and GGA+U calculations, Physical Review B, vol. 84 (2011), 045115.
  7. Setyawan, W.; Curtarolo, S. Computational Materials Science 2010, 49, 299-312.
  8. J.B. Foresman, A.E. Frisch, Exploring Chemistry With Electronic Structure Methods: A Guide to Using Gaussian, Gaussian. (1996).
  9. A. Jain, S.-a Seyed-Reihani, C.C. Fischer, D.J. Couling, G. Ceder, W.H. Green, Ab initio screening of metal sorbents for elemental mercury capture in syngas streams, Chemical Engineering Science. 65 (2010) 3025-3033.
  10. S. Lany, Semiconductor thermochemistry in density functional calculations, Physical Review B. 78 (2008) 1-8.
  11. G. Hautier, S.P. Ong, A. Jain, C. J. Moore, G. Ceder, Accuracy of density functional theory in predicting formation energies of ternary oxides from binary oxides and its implication on phase stability, Physical Review B, 85 (2012), 155208
  12. E.N. Brothers, A.F. Izmaylov, J.O. Normand, V. Barone, G.E. Scuseria, Accurate solid-state band gaps via screened hybrid electronic structure calculations., The Journal of Chemical Physics. 129 (2008)
  13. R. Godby, M. Schluter, L.J. Sham, Self-energy operators and exchange-correlation potentials in semiconductors, Physical Review B. 37 (1988).
  14. M. Chan, G. Ceder, Efficient Band Gap Predictions for Solids, Physical Review Letters 19 (2010)
  15. J. Heyd, J.E. Peralta, G.E. Scuseria, R.L. Martin, Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional, Journal of Chemical Physics 123 (2005)
  16. W. Setyawan, R.M. Gaume, S. Lam, R. Feigelson, S. Curtarolo, High-throughput combinatorial database of electronic band structures for inorganic scintillator materials., ACS Combinatorial Science. (2011).

Authors

  • Anubhav Jain
  • Shyue Ping Ong
  • Geoffroy Hautier
  • Charles Moore