Understanding the difference in cohesive energies between alpha and beta tin in DFT calculations

The transition temperature between the low-temperature alpha phase of tin to beta tin is close to the room temperature (Tab =13C), and the difference in cohesive energy of the two phases at 0 K of about dEcoh=0.02 eV/atom is at the limit of the accuracy of DFT (density functional theory) with available exchange-correlation functionals. It is however critically important to model the relative phase energies correctly for any reasonable description of phenomena and technologies involving these phases, for example, the performance of tin electrodes in electrochemical batteries. Here, we show that several commonly used and converged DFT setups using the most practical and widely used PBE functional result in dEcoh of about 0.04 eV/atom, with different types of basis sets and with different models of core electrons (all-electron or pseudopotentials of different types), which leads to a significant overestimation of Tab. We show that this is due to the errors in relative positions of s and p -like bands, which, combined with different populations of these bands in alpha and beta Sn, leads to overstabilization of alpha tin. We show that this error can be effectively corrected by applying a Hubbard +U correction to s -like states, whereby correct cohesive energies of both alpha and beta Sn can be obtained with the same computational scheme. We quantify for the first time the effects of anharmonicity on dEcoh and find that it is negligible.


Introduction
Tin (Z=50) is part of the group IV materials. Similarly to its lighter counterparts (C, Si, and Ge), Sn exists in the diamondtype crystal structure, commonly called alpha Sn (α-Sn, space group 227; Fd-3m) or grey tin; this is the most stable lowtemperature phase, and it is a zero-gap semiconductor. 1 Unlike C, Si, and Ge, however, the diamond structure of Sn is only stable under 286 K (13°C), above which Sn transforms to beta Sn (β-Sn, space group 141; I41/amd), the metallic and room temperature phase of tin. The relatively low transition temperature T αβ is reflective of an only slight stabilization of α Sn over β Sn at 0 K (estimated at about ∆E coh =0.02 eV/atom), so that relatively small vibrational contributions are required to reverse the phase stability. The proximity of T αβ to room temperature means that the α−β transition can be of practical importance, as the relative phase stability can be easily changed by perturbations to the lattice by e.g. doping or strain. An example of a technology where this transition has been found important are electrochemical batteries. Tin has been shown to be a high capacity anode material for Li, Na, and Mg ion batteries. [2][3][4][5] Specifically in Li ion batteries, for which Snbased electrodes have been most studied, 6,7 the formation of α Sn upon lithiation of (initially β phase) Sn electrodes has been reported in experimental studies. 8,9 This beta-to-alpha phase transformation in battery electrodes is still not fully understood. To produce reliable ab initio models of this process or other processes which can affect the phase stability at room temperature, the transition temperature and therefore ∆E coh need to be reproduced accurately.
The practical approach to ab initio modelling of solids is DFT (density functional theory) with PBC (periodic boundary conditions), and GGA (generalized gradient approximation) 10 with the PBE (Perdew-Burke-Ernzerhof) 11 functional is the most practical and widely used approximation to the exchange-correlation functional. 12,13 Hybrid functionals 14,15 can be used with small simulation cells 16 but are increasingly unwieldy for larger supercells which are required to study doping, intercalation, and interfaces. A difference in cohesive energy on the order of 10 -2 eV/atom is at the limit of DFT accuracy with available exchange-correlation functionals. Previous DFT simulations reported ∆E coh of 0.022-0.047 eV/atom. 17,18 Below, we show that converged calculations using the PBE functional and different types of basis sets and core electron treatment invariably result in ∆E coh ≈0.04 eV/atom. This leads to an overestimation of T αβ by more than 100 K if the harmonic approximation for the phononic contribution is used and is therefore not suited for modelling near-room temperature behaviour of tin-based materials. The transition temperature is the temperature at which the sum of the electronic ( ) and of the vibrational energy ( ) and entropy (− ) contributions, is equal in the two phases. Eq. (1) is an approximation for the Gibbs free energy that neglects the pV term; this is appropriate for low pressure conditions including standard conditions. At the transition temperature, the difference in vibrational contributions to G counterbalances exactly the DFT energy difference, i.e. Δ ℎ = Δ( − ). This means that the error in T αβ can come either from the error in electronic energy or/and in the vibrational contribution, on the order of 0.02 eV/atom. The difference in vibrational contribution to G at the experimental transition temperature between alpha and beta Sn is around 0.024 eV based on GGA DFT and the harmonic approximation (see below).
With an accurate computational setup, each of the two terms of Eq. (1) should therefore be exact to within a few meV, and not within a few dozen meV as is the case with the standard setup using the PBE functional and the harmonic approximation. It is also clear that an accurate setup must be able to model correctly both α and β Sn simultaneously to model systems where the phase transition is possible. It means that fixes such as selecting the degree of exact exchange 15,19 or Hubbard (+U) corrections 20 that tune the cohesive energy and band structure with parameters different for the two phases are not useful.
In this paper, we investigate the origin of the error in T αβ and identify a computational scheme that fixes it. We analyse both the electronic and the phononic contributions to free energy. Specifically, for the first time, we analyse the effect of anharmonicity of vibrations of bulk tin on relative phase stability; it is shown to be unimportant. We show that the computed ∆E coh of about 0.04 eV is ubiquitous across DFT setups using localized or plane wave based basis sets and full ionic potential or pseudopotentials. We trace the origin of the error in the electronic energy to destabilization of the energy of the s band vs the p band and show that the correct ∆E coh can be obtained by applying a Hubbard +U correction to the atomic-part s orbitals, with which absolute and relative cohesive energies of both alpha and beta tin can be accurately modelled with one and the same computational setup.

Methods
DFT calculations were performed with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional. 11 Different codes were used employing (augmented) plane wave (VASP 21 and Elk 22 ) or atom-centered (SIESTA 23 and FHI-AIMS 24 ) basis sets, full ionic potential (FHI-AIMS and Elk) or pseudopotentials (VASP and SIESTA). When using pseudopotentials, the valence configuration was 5s 2 5p 2 ; the effect of the addition of the 4d electrons to the valence shell was found to be insignificant. Spin polarized calculations were used where spin polarization is important (i.e. tin atoms and dimers) and spin restricted where it is not (i.e. bulk phases).
The VASP 21 calculations were performed using a PAW (projector-augmented wave) pseudopotential 25 and a plane wave energy cutoff of 300 eV. Cubic 8 atom cells were used for both α and β Sn. 16×16×16 Г-centred Monkhorst-Pack meshes were used for the k-point sampling of the Brillouin zone integration for α and β Sn. A Gaussian smearing is used with σ=0.005. Atomic coordinates and cell vectors were optimized until the residual forces were below 0.01 eV/Å. To consider the vibrational contributions, we computed the Hessian matrices on cubic 64-atom cells (8×8×8 k-mesh) using the displacement method and a step of 0.01 Å. The Elk 22 calculations were performed using the full potential linearized augmented plane wave method (FP-LAPW). This approach implements a solution to an all-electron and a full-potential system. The 36 energetically lowest electrons are considered to be core electrons and the remaining 14 electrons are valence electrons. The core electrons are treated fully relativistically in the spincompensated Dirac equation for spherical potential, and the valence electrons are treated within the scalar-relativistic approximation where the spin-orbit coupling is included as a perturbation. The muffin-tin radius is R = 2.3 a.u. 26 The planewave cutoff K max was determined from R•K max = 10.0 a.u. The angular momentum cutoff for the APW functions and for the muffin density and potential was set to 10 and 9, respectively. Maximum length of the reciprocal vectors G max for the expansion of the interstitial density and potential was get to 14.0 a.u. Total energy convergence was set to <0.5×10 -6 a.u./atom, and the electron population function was broadened with the Fermi-Dirac distribution with width 0.001 a.u. Eight-and four-atom cells were used for α and β Sn, respectively, and 10x10x10 and 10x10x20 k points, respectively, provided converged results. The SIESTA 23 calculations were performed using a Troullier-Martins pseudopotential 27 (generated by us) and DZP and TZP (double-ζ polarized and triple-ζ polarized, respectively) basis sets. The DZP basis was tuned by using the PAO.EnergyShift parameter to reproduce both the absolute value of E coh (to the experimental value, i.e. ≈3.14 eV/atom) 28 and the difference between the phases, ∆E coh (to the converged plane wave values, i.e. 0.04 eV/atom). We then checked that the conclusions hold for a TZP basis made with the same value of PAO.EnergyShift and are therefore not polluted by a particular choice of the atom-centred basis. A 200 Ry cutoff was used for the expansion of the density with a bcc-type oversampling of the Fourier grid. Eight-and four-atom cells were used for α and β Sn, respectively, with 10x10x10 and 10x10x20 k points, respectively, for most calculations; we confirmed that doubling the number of k points in each direction does not change the result. Atomic coordinates and cell vectors were optimized until the forces on atoms were below 0.01 eV/Å and stresses below 0.1 GPa. A development version of SIESTA which includes the LDA+U functionality 29 with the Dudarev approximation 30 was used to apply Hubbard correction on s or p -like levels.
The FHI-AIMS 24 calculations were performed with basis sets and integration grids set to the "really_tight" settings to approach the basis set limit. The convergence criteria for the self-consistency cycle were 1×10 -7 eV for the energy and 1×10 -6 for the density. Unit cells were used (with 2 atoms per cell, i.e. non-rectilinear). A smearing of σ=0.1 eV was used. 15 kpoints per dimension provided converged Brillouin zone sampling (change in energy of either phase of less than 0.001 eV vs. 30 k-points). Atomic coordinates and cell vectors were optimized until the residual forces were below 1×10 -2 eV/Å.
The atomic ZORA 31 approximation was used to account for relativistic effects.
CCSD(T) 32 calculations were performed in Gaussian 09 (G09). 33 Calculations were performed on Sn atoms and dimers using the LANL2DZ 34,35 basis set, and the trends in the differences in s and p levels between CCSD(T) and DFT were also confirmed with calculations on atoms using a larger QZVPPD basis set with a (small-core) ECP28MDF ECP. 36 DFT calculations were also performed in Gaussian 09 on atoms and dimers for comparison with CCSD(T) calculations.
Contributions to the free arising from the vibrational energy and entropy were computed as: where ν i is the energy of one quantum in the i-th normal mode, T is the temperature, and k B is the Boltzmann constant.
The sum is over all N normal modes. Adding these vibrational contributions to the DFT ground state energy E DFT gives the free energy of the system, Eq. (1). Molecular dynamics (MD) simulations were performed in VASP with a time step of 2 fs. Cubic cells with 64 Sn atoms were used for both phases. The k-grid was reduced to 3x3x3 for alpha and 4x4x4 for beta Sn. The temperature was ramped from 0 K to 300 K over 2000 (3000) steps for β (α) Sn. A total of 15,000 steps were computed for α and 10,000 steps for β Sn.

Alpha and beta tin with different DFT setups
Crystal structures of α and β tin are shown in Fig. 1. Table 1 lists α Sn cohesive energies as well as the differences in cohesive energies between α and β Sn obtained with different DFT codes, basis set types and core electron treatments. A difference in cohesive energies of about 0.04 eV/atom is persistent among different computational setups and is different by about 0.02 eV/atom from experimental estimates of about 0.02 eV/atom. 37 This results in a significant overestimation of the transition temperature, e.g. value T αβ = 475 K is obtained in VASP using the harmonic approximation, Eq. (2). We have confirmed that using the quasiharmonic approximation, namely, fixing the cell volume according to the thermal expansion coefficient of Sn, does not significantly affect the result, changing the transition temperature by only 45 K (as computed in VASP). The error in ∆E coh and T αβ could come mainly from the electronic or vibrational contributions to the free energy, or both. In what follows, we therefore analyse each of these contributions.

Effects due to anharmonicity
The vibrational part of free energy as expressed in Eq. (2) neglects the effects of anharmonicity and coupling on vibrational levels. To estimate these effects, we performed MD simulations of both phases of tin to sample the potential energy surface up to potential energies relevant for T αβ . The collected potential values, V DFT , were compared to the harmonic uncoupled potential built on normal mode frequencies: where Q i are normal mode coordinates related to the Cartesian coordinates of Sn atoms via where Q is a vector of all normal mode coordinates, L is the matrix that diagonalizes the Hessian matrix, x is a vector of Cartesian coordinates of all atoms, and x eq are equilibrium positions. M is a diagonal matrix of atomic masses. We consider N=3N A , where N A is the number of atoms in the simulation cell, in which case all matrices are square of size NxN and vectors are of length N. This therefore includes (zerofrequency) frustrated translation modes.  Fig. 2 shows the plots the harmonic potential energies V harm vs. anharmonic energies (V DFT ) following from the MD simulations for both phases. The onset of anharmonicity around room temperature (corresponding to ≈ 2.5 eV) can be visually appreciated. Does it significantly contribute to relative phase stability? In the classical limit (i.e. all → 0) there would be no effect of anharmonicity and coupling in the potential on relative phase stability, while in the limit of large (specifically, when ≫ ), changes in would directly contribute to changes in G. The spectrum of spans the range up to about 175 cm -1 for α and up to about 130 cm -1 for β Sn i.e. all components are smaller than or on the order of at room temperature (see Fig. 3 for phonon DOS). To estimate the effect of anharmonicity on G, we express Eq. (2) as where the inner summation uses a sufficiently large n i , max for convergence, and , are the vibrational energies with n i quanta of excitation in mode i. We perform a non-linear regression of the dependence of V harm on V DFT , ℎ = ( ) as shown in Fig. 2 and apply the function f to estimate , : As changes in vibrational levels are typically only a fraction of changes in the potential energy values, 39 Eq. (6) can be considered as an upper limit of the effect due to anharmonicity. It results in a change in the relative energies of the two phases on the order of 0.001 eV/atom at room temperature. This would only change the transition temperature by about 10 K and is insignificant compared to the difference of about 0.02 eV/atom between the computed ∆E coh and Δ( − ) 290 . We therefore conclude that anharmonicity and coupling of tin vibrations are not the cause of the overestimation of T αβ reported in the previous subsection and that ∆E coh must be overestimated by DFT.

Analysis of electronic structure
To understand why there is persistent overestimation of ∆E coh with GGA DFT (PBE functional) regardless of the type of the basis set or pseudopotential, we turn to the differences in electronic structure between alpha and beta Sn and analyse how well they are modelled by DFT. From the Mulliken population analysis (performed with atom-centered basis set calculations in SIESTA), we find that α Sn has a higher degree of transfer of electrons from s to p population than β Sn when forming bonds in bulk tin: the s population in α Sn is about 1.5e and in β Sn about 1.7e (to compare to 2.0e in a free atom, as per the valence shell configuration 5s 2 5p 2 ). The differences in s vs p population directly contribute to the total energy, and thereby to ∆E coh , via the bandstructure component of the total energy. The population of orbitals with d character upon bulk formation is found to be insignificant. Fig. 4 shows projected density of states for both phases of tin.
In Table 2, we compare the average energies of s and p bands (computed as population-weighted energies of states with s and p character) and the difference between them (p-s "gap") in a Sn atom and in Sn 2 dimers computed with DFT to those computed with CCSD(T) as a reference. There is a similar pattern of differences (errors) in s and p -like band energies computed with DFT in all codes we used vs. CCSD(T); namely, both s and p -like levels computed with DFT are higher in energy (as expected) than the corresponding CCSD(T) values. More importantly, the s -like level is destabilized by a larger amount, and therefore the p-s gap is lower than with CCSD(T). When bonds are formed between tin atoms, in CCSD(T), the s band energy is lowered by 0.83 eV while the p band energy is slightly raised by 0.1 eV. In the DFT calculations, both s and p levels are stabilized by bonding, and p levels are stabilized more than s levels.
This qualitative difference is observed in the DFT calculations in all Gaussian 09, VASP, and SIESTA and is therefore not an artefact caused by a specific computational setup. In the periodic codes, we have confirmed that the s and p level energies are converged to within 0.01 eV in SIESTA and 0.025 eV in VASP with respect to the size of the cubic supercell (total energies and s-p gaps are converged much faster with respect to cell size than eigenvalues). It is therefore expected that in materials where Sn-Sn bonding is important, s and p band energies will be in error, and by a different amount. If the stronger destabilization of s vs p band with DFT using the PBE functional which is observed in a Sn atom and Sn dimer carries over into tin bulk, then one expects that a phase with a higher p population will be additionally (and erroneously) stabilized vs. a phase with a lower p population. This is exactly what is observed in alpha vs. beta Sn.
Conversely, it must be possible to reverse this effect by stabilizing the s-like band. This can be achieved by applying a positive Hubbard correction or U value to s-like states. The effect of lowering the s band energy and of increasing the p-s gap was indeed confirmed in Sn atom and dimer calculations. For example, applying U(s) = 1.5 eV in SIESTA reduces the destabilization of the s band of a Sn atom vs. CCSD(T) ( Table 2) from 3.3 eV to 2.8 eV while the destabilization of the p states remains at 3.1 eV; the s band destabilization in the dimer is reduced from 4.1 eV to 3.6 eV and the degree of underestimation of the p-s gap from 1.6 eV to 1.1 eV, while the p band destabilization remained at 2.5 eV. A similar effect is observed in VASP (see Appendix). We will see below that a similar order of magnitude of U is required to correct for the over-stabilization of the alpha vs. beta phase.  We have then applied different +U values to s -like states of both alpha and beta tin to test if this conjecture holds for bulk Sn. With the increase of the U value, ∆E coh indeed decreased. When using a localized basis set in SIESTA (which was tuned to E coh of α Sn), the application of +U led to a modest weakening of E coh (within 0.1 eV/atoms for U values of up to +2.0 eV) and only a small change of the lattice constant (on the order of 0.01 Å), see Table 1. This can easily be compensated by tuning the PAO.EnergyShift parameter. For example, with a DZP basis obtained with PAO.EnergyShift=0.002 Ry and giving E coh (α)=3.13 eV/atom (with +U), a U value of +1.5 eV results in ∆E coh = 0.025 eV/atom. To confirm that the effect of +U is not an artefact of the choice of the atom-centered basis set, we performed a calculation with a TZP basis (using the same PAO.EnergyShift value) and obtained E coh (α)=3.15 eV/atom, ∆E coh = 0.018 eV/atom with U=+1.25 eV and E coh (α)=3.17 eV/atom, ∆E coh = 0.022 eV/atom with U=+1.0 eV. With U=0, ∆E coh was 0.039 eV/atom with the same TZP basis.
The effect of a +U value applied to s states on the band structure can be seen from the PDOS (partial density of states) plot in Fig. 4, where we plot the density of s-like and p-like state. The PDOS are plotted up to the Fermi level, and E=0 refers to the average electrostatic potential. As expected, the +U value lowers the energy of the s band and affects little the p band. The population averaged energy of the s band (i.e. ) changes from -8.85 eV to -9.38 eV in α Sn and from -10.08 eV to -10.66 eV in β Sn when going from U=0 to U(s)=+1.75 eV. Simultaneously, the p-s gap increases from 0.26 to 1.06 eV in α Sn and from 2.49 to 2.99 eV in β Sn. The effect on the Fermi level is minimal: within 0.05 eV for both phases when U is increased from 0 to 1.75 eV; the changes in the average electrostatic potential are negligible when applying a U value of this magnitude.
We have therefore shown that (i) the application of +U affects similarly the band structure of tin atom and dimers and of bulk Sn; (ii) in tin atoms and dimers, the s band energy is destabilized more than the p band energy by DFT vs CCSD(T); (iii) alpha tin has a smaller fraction of s electrons than beta tin. Together (i)-(iii) allow us to conclude that (i) DFT with the PBE functional over-stabilizes alpha over beta phase of tin by the mechanism of destabilizing s band more than the p band; (ii) this can be fixed by applying a positive U value to s-like states which permits obtaining accurate cohesive energies of both phases with the same computational setup.
We have confirmed this result by performing calculations in VASP using a different type of basis set and of pseudopotential (i.e. plane waves and PAW, respectively). The results are also summarized in Table 1. Similarly to SIESTA calculations, the s population is higher in the beta phase by about 0.2 e/atom. The application of U=+1 eV on s levels does lead to the shift down in energy of the s-like band while inducing no significant shift in the p-like band (see Fig. A1 in the Appendix) and results in the change in ∆E coh from 0.039 eV/atom to 0.023 eV/atom. This simultaneously improves the cohesive energy of α Sn from 3.20 eV/atom to 3.16 eV/atom. That is, also in VASP we are able to choose a +U(s) value that results in correct cohesive energies of both α and β Sn. Similar values of +U(s) achieve ∆E coh of about 0.02 eV/atom in both VASP and SIESTA. That optimal values of +U a little different in VASP and SIESTA is not surprising given the somewhat differing definitions of s and p character with plane waves and with atom-centred bases and otherwise significant differences in the computational approaches and implementation. What is important is the common to these setups mechanism of destabilization of s levels vs p levels which leads to an overestimation of ∆E coh and which can be fixed by applying a +U value on s-like levels which allows reproducing accurate cohesive energies of both alpha and beta tin with one and the same computational scheme.

Conclusions
In the ab initio modelling of tin based materials, specifically doped tin and alloys, it is necessary to be able to correctly model the energetics of the two low temperature and low pressure phases, namely, alpha and beta Sn, with the same computational setup. The small difference in the cohesive energies of these phases (about ∆E coh =0.02 eV/atom) and the resulting proximity of phase transition temperature to the room temperature make such modelling difficult. Specifically, the most practical and by far the most widely used approximation -DFT with the PBE functional -overestimates the difference in cohesive energies. We have shown that values of ∆E coh ≈0.04 eV/atom are obtained with different types of basis sets and with different models of core electrons (all-electron or pseudopotentials of different types), which also leads to a significant overestimation of the α−β transition temperature.
To identify the reason for this overestimation, we have analysed for the first time possible contributions from anharmonicity to the vibrational part of the free energy as well as errors in the electronic structure. We have shown that anharmonicity does not have a significant effect on T αβ . By comparing DFT and CCSD(T) calculations on small model systems, we have shown that DFT with the PBE functional overestimates the energy of the s band to a larger extent than that of the p band. As s population is higher in beta Sn than in alpha Sn, this leads to overstabilization of α vs β Sn in DFT calculations. We have shown that it is possible to correct for this by applying a positive U value on s levels, whereby correct cohesive energies of both α and β Sn can be obtained with the same computational scheme.
The cohesive energies of α-Sn and β-Sn is therefore a highly practically relevant example of the need for developing of computationally efficient exchange correlation functionals that perform well for different phases of the same material with qualitatively different electronic structure (such as metalling and non-metallic phases).