# A density functional study of lithium bulk and surfaces

## I Introduction

Lithium has been the subject of considerable interest over many years. Although its electronic structure is relatively simple, its structural properties still pose a significant challenge to both experiment  and simulation . Lithium is very soft; the determination of its elastic constants and surface energies requires experiments of high accuracy and simulations of high numerical stability. The calculation of surface formation energies is particularly delicate as has been discussed recently .

The aim of this article is to present a comprehensive and systematic study of the band structure, cohesive energy, elastic constants, phase stability, and surface energies of lithium. Very few systematic studies of the dependence of results on the computational parameters are available. This however is especially important when energy differences are required as for example in calculation of surface energies or phase stabilities. We present the results of fully converged, full potential, all-electron calculations based on both density functional theory and the Hartree-Fock approximation. We examine the use of finite temperature density functional theory as a technique for accelerating convergence with respect to reciprocal space sampling. We expand the crystalline orbitals as a linear combination of Gaussian type orbitals (LCGTO). This approach is very well established for insulators . We find, in accord with a recent study of magnesium , that this approach is also well suited to the simulation of a free-electron metal.

The paper is organized as follows: in section II, we discuss the computational parameters. In sections III and IV, we discuss results for lithium bulk and surfaces, respectively, and summarize our conclusions in section V.

## II Basis set and Method

All the calculations were performed with the program package crystal . The main numerical approximation in our approach is the choice of the Gaussian basis set. The difficulties of selecting basis sets for metallic systems have been explored in previous studies . In principle the quality of a calculation can be systematically improved by adding additional functions to the basis set and optimizing their exponents in a suitable reference state - usually the bulk crystal. In practice one must balance the overcompleteness of the basis set, which leads to linear dependence, with the need for additional variation freedom. For molecular systems and insulating solids these problems have largely been overcome and Gaussian basis sets are very widely used. For metallic systems and their surfaces in particular there have been very few systematic studies.

In a solid the tails of the atom centered functions overlap strongly and so diffuse basis functions optimized for the description of atomic or molecular systems are not useful and indeed may give rise to linear dependence. We are thus unable to simply use basis sets from the many libraries developed for the description of molecular systems. We have therefore developed an hierarchy of basis sets of increasing quality in order to examine the convergence of computed properties.

The smallest basis set used has 3 $`s`$-symmetry functions and 2 $`p`$-symmetry functions and is denoted as $`[3s2p]`$ . The $`[1s]`$ radial function was taken from Ref. . The exponents of the two additional $`sp`$-shells were optimized in local density approximation (LDA) calculations (with Dirac-Slater exchange  and Perdew-Zunger’s correlation functional ) for the solid at the experimental lattice constant. The lowest energy was obtained with exponents of 0.50 and 0.08. Similar results were obtained when using the Perdew Wang gradient corrected approximation (PWGGA). However, as an exponent of 0.08 gives rise to a very diffuse basis function close to numerical instability, instead exponents of 0.50 and 0.10 were chosen. This $`[3s2p]`$ basis set is very robust and computationally efficient — it does not give rise to linear dependence even when the bulk is strongly distorted (for example, to determine the elastic constants).

A $`[4s3p]`$ basis set was obtained by using three exponents (0.50, 0.20, and 0.08) - which were chosen to be ”even tempered”, i.e. the ratio between the exponents is kept fixed (2.5 in this case). This ratio is close to the lowest which can be tolerated before on-site (atomic) linear dependence is seen. It is however also known to converge the atomic energy to within less than $`10^4E_h`$ ($`E_h`$ = 27.2114 eV) of the exact Hartree-Fock ground state energy (see the analysis in Ref. ). Finally, an additional polarization function of $`d`$-symmetry was added and the exponent optimized within a PWGGA calculation to be 0.15. However, the $`d`$-function leads only to a minor change in the total energy. The energy varies only by $`5\times 10^5E_h`$ when, e.g., changing the exponent to 0.5. The basis sets developed in this manner are displayed in table I.

Both at the Hartree-Fock (HF) and B3LYP  (invoving a hybrid of Fock exchange and a modification of the Becke gradient corrected exchange functional, and the Vosko-Wilk-Nusair local correlation functional V  and the gradient corrected correlation potential by Lee, Yang and Parr ) levels, an optimization of basis set exponents was not possible. Instead, the outermost exponent became more and more diffuse until finally the solution became unstable. This is a well known pathology of the use of Fock exchange in metallic systems (see also the discussion in Ref. ). When features of the HF solution are discussed in this article, they were obtained with the $`[3s2p]`$ basis set (outermost exponents 0.50 and 0.10).

In order to compute binding energies the free atom is calculated within a spin-polarized formalism with the same $`[1s]`$ function but with additional $`s`$-exponents 0.60, 0.24, 0.096, 0.04, and 0.016 to describe the long range behaviour of the atomic wavefunction.

For the LDA and PWGGA calculations we also expand the exchange and correlation potentials in an auxiliary Gaussian basis set which consists of 13 even tempered $`s`$-functions with exponents from 0.1 to 2000, 3 even tempered $`p`$-functions with exponents from 0.1 to 0.8, and 2 $`d`$-functions with exponents of 0.12 and 0.3. This is sufficient to integrate the charge density to an accuracy of 10<sup>-7</sup> $`|e|`$. For the free atom, we use an auxiliary basis set with 18 even tempered $`s`$-functions with exponents from 0.0037 to 4565.

Reciprocal space sampling is a delicate problem especially in metals. The sampling is performed on a Pack-Monkhorst net  where the density of points is determined by a shrinking factor. The Fermi energy and shape of the Fermi surface are determined by interpolation onto a ”Gilat” net. This net is simply related to the Pack-Monkhorst net by an additional subdivision factor. To further improve convergence, the finite temperature generalization of density functional theory can be used to apply Fermi surface smearing . In table II, the dependence of the total energy on the density of points in the Pack-Monkhorst net is displayed. For this purpose, PWGGA calculations on a body centered cubic (bcc) lattice at the equilibrium lattice constant of 3.44 Å  were performed.

At zero temperature, even with the largest net used, the energy is still slightly decreasing as more points are used. For a smearing of 0.001 $`E_h`$ (which corresponds to a temperature of $`T=\frac{0.001E_h}{k_B}=316K`$ with Boltzmann’s constant $`k_B=3.1667\times 10^6\frac{E_h}{K}`$), the energy is stable up to a few $`\mu E_h`$. A higher number of sampling points in the Gilat net leads to a systematic improvement at zero temperature. At finite temperature, the number of sampling points in the Gilat net does not influence the results when a sufficiently high number in the Pack-Monkhorst net is chosen. As shown in table II, the difference in energy for the different number of sampling points in the Gilat net is of the order of only a few $`\mu E_h`$ for a fixed shrinking factor of 24 in the Pack-Monkhorst net. Note that for the most dense net at $`k_BT=0.001E_h`$, the difference in energy between smeared and unsmeared results is less that $`10^4E_h`$. At an even higher temperature of 0.02 $`E_h`$, the energy converges to a value which is 6 $`mE_h`$ higher than at 0.001 $`E_h`$. An estimate of the zero temperature energy is possible by using the approximation $`E(0K)=\frac{1}{2}(E(T)+F(T))`$ (with $`F=ETS`$ being the free energy and exploiting the fact that the energy increases quadratically with temperature for low temperature) as suggested in Ref. . The electronic entropy $`S`$ is defined as   
$`S=k_B_i^{N_{states}}f_i\mathrm{ln}f_i+(1f_i)\mathrm{ln}(1f_i)`$   
with $`f_i`$ being the Fermi function. This leads to a value of $`E(0)`$ extrapolated from $`k_BT=0.02E_h`$ which deviates by less than $`10^3E_h`$ from the value obtained at a temperature of 0.001 $`E_h`$. The functions $`E(T)`$, $`F(T)`$ and $`\frac{1}{2}(E(T)+F(T))`$ for a fixed value of 145 sampling points (corresponding to a shrinking factor of 16 in the Pack-Monkhorst net) are also displayed in Figure 1. Indeed, even for relatively high temperature, $`E(0)`$ is well approximated by $`\frac{1}{2}(E(T)+F(T))`$.

In conclusion, as we need a high accuracy for the energy difference between the different phases of Li, we chose for the calculations on lithium bulk a shrinking factor of 16 for the Pack-Monkhorst net which gives 145 sampling points in the irreducible Brillouin zone of the non-distorted bulk and a temperature of 0.001 $`E_h`$. This ensures that convergence of the energy to at least $`10^4E_h`$ with respect to reciprocal space sampling is achieved.

## III Results for bulk Lithium

### A Band structure

In figure 2 the LDA band structure for the bcc structure is displayed. The occupied bands and the lower unoccupied bands are in excellent agreement with results earlier obtained with Gaussian basis sets and $`X_\alpha `$ exchange , by the Kohn-Korringa-Rostoker (KKR) method  and Slater exchange, augmented plane wave (APW)  and modified APW (MAPW)  calculations both using LDA, or linear muffin tin orbital (LMTO) calculations  with a combination of exchange using the Langreth and Mehl functional  and LDA correlation. When using PWGGA instead of LDA, the band structure does not exhibit major differences. The experimentally known conduction band width ($``$ 4 eV)  is slightly lower than the result calculated here (4.6 eV). The slow decay of the density matrix in metals leads to difficulties when Fock exchange is involved: the summation of the exchange series in direct space is very long ranged and is truncated at a large but finite distance. This cutoff for large distances in direct space results in numerical instabilities for small $`\stackrel{}{k}`$. Thus, both in the Hartree-Fock and B3LYP band structures, artificial oscillations can be found around the Gamma point. As usual for Hartree-Fock calculations, we find that the bandwidth is roughly twice as large as the experimental bandwidth.

### B Cohesive properties

Table III gives results for ground state properties of bcc Li (lattice constant, bulk modulus, cohesive energy and elastic constants C<sub>11</sub>, $`C_{11}C_{12}`$, and C<sub>44</sub>). The elastic constants were obtained by applying a rhombohedral distortion for $`C_{44}`$, a tetragonal distortion for $`C_{11}`$, and an orthorhombic distortion for $`C_{11}C_{12}`$ to the solid at the equilibrium lattice constant. Cohesive energy and lattice constant agree well with experiment and depend only weakly on the basis set. The bulk modulus is obtained from the energy as a function of volume and agrees within the accuracy of the fit with that obtained using the relation $`B=\frac{1}{3}(C_{11}+2C_{12})`$. Elastic constants have a strong dependence on the basis set and the deviation from experiment improves especially when going from $`[3s2p]`$ to $`[4s3p]`$; the $`d`$-function has only minor impact on the results. We find that PWGGA is closer to experiment with results similar to Ref. .

### C Relative stabilities

The experimental crystal structure of lithium at zero temperature is still unclear (at room temperature a bcc structure is favoured). Both face-centered cubic (fcc) and hexagonal (hex) structures have been suggested as well as more sophisticated structures such as 9R (a nine layer sequence of close-packed planes ABCBCACAB)  or mixtures of these phases. Therefore, we also investigated the relative stability of the different phases. In Figure 3, total energies of bcc, fcc and hex phase, obtained with the PWGGA functional and the best basis set ($`[4s3p1d]`$), are displayed.

The c/a ratio of the hexagonal phase remains close to the ideal close packed value of 1.633 (it varies between 1.631 and 1.635 which is within the accuracy of the fit). We find that the closed-packed structures are slightly lower in energy than the bcc structure. Although it is reasonable to conclude that the close packed structures are lower in energy than bcc it is not possible to resolve the difference in energy between the fcc and hex phases. The variation of the energy with basis set is given in table IV. As the basis set is systematically improved the energy difference between bcc and fcc increases from 1 $`\times 10^5E_h`$ to 4 $`\times 10^5E_h`$, and between bcc and hex decreases from 7 $`\times 10^5E_h`$ to 4 $`\times 10^5E_h`$. We note that even the $`d`$-function influences this energy splitting. These results, which do not include the zero point energy, indicate a preference for closed-packed structures which is in agreement with most of the previous calculations  except for one  (see table IV).

HF calculations were only possible with the smallest $`[3s2p]`$ basis set where the order of phases is different with hcp being the lowest in energy, followed by bcc and fcc being highest. The same is found in LDA with the smallest $`[3s2p]`$ basis set, but changes when the basis set is increased to $`[4s3p]`$ where both fcc and hcp are 2 $`\times 10^4E_h`$ lower in energy than bcc (again, calculations with the best $`[4s3p1d]`$ basis set were not possible because of linear dependence) — table IV.

As the energy splitting between the close packed and bcc phases is rather small the zero-point energy difference cannot be neglected. The first published calculation of this, based on the harmonic approximation,  gave an additional stabilization of $`9\times 10^5E_h`$ of the fcc phase compared to the bcc phase. In these calculations however the hcp phase was found to be higher in energy that the bcc phase (table IV). Very recently, in a calculation also including anharmonic effects, the stabilization was calculated to be about $`1.5\times 10^5E_h`$ both for fcc and 9R phases relative to the bcc phase . In addition, the authors computed the variation of the vibrational free energy as a function of temperature and found a phase transition from a closed packed to a bcc phase at a temperature of $`T200K`$.

## IV Results for surfaces

A further test of this approach is the calculation of surface energies. We model lithium surfaces by using a slab and varying the numbers of layers (a bcc structure is assumed). Surface energies can be calculated in two ways, either by deriving a bulk energy by subtracting energies of two slabs with $`n`$ and $`m`$ layers:

$$E_{surface}=\frac{1}{2}(E_{slab}(n)(E_{slab}(n)E_{slab}(nm))\frac{n}{m})$$
(1)   
which has the advantage of a systematic error cancellation (in particular the reciprocal space sampling is consistent between the bulk and slab energies) or by using an independent bulk energy

$$E_{surface}=\frac{1}{2}(E_{slab}(n)E_{bulk}\times n)$$
(2)   
All the quantities $`E_{surface}`$, $`E_{slab}(n)`$, and $`E_{bulk}`$ are energies per atom.

In Figure 4, results for surface energies using both equation 1 (with $`m`$=1) and 2 are displayed. Equation 1 leads to relatively strong oscillations (dotted line with pluses) and gets more stable the larger $`m`$. Numerical noise in the expression $`(E_{slab}(n)E_{slab}(nm))\frac{n}{m}`$ is reduced for larger values of $`m`$. Equation 2, however, leads at zero temperature to a slight linear decreasing of the surface energy as a function of the number of layers (thin line without additional symbol). The reason for the nonvanishing slope is that the energy difference $`E(n)E(n1)`$ is not identical to the energy of the bulk due to the systematic errors in the convergence of the total energy with respect to reciprocal lattice sampling (this was also emphasized in Refs. ). As shown in table II, at zero temperature the bulk energy is still changing on the order of $`10^4E_h`$ with increasing sampling point number. Similarly the bulk energy varies when extracted from the slab. This slight discrepancy gives rise to a variation of the surface energy with the number of layers with a slope of $`10^4E_h/atom/layer`$.

The origin of the poor convergence with respect to reciprocal lattice sampling is due to the sharp cutoff imposed by the Fermi energy. One possibility to obtain the surface energy would be to use the intercept from a linear fit, but a better and simpler way to alleviate the difficulties associated with reciprocal lattice sampling is smearing the Fermi surface with a finite temperature. Already at a temperature of 0.001 $`E_h`$, the slope is virtually zero (thick line without additional symbol). This is consistent with table II as the bulk energy converges much faster at finite temperature. At a higher temperature of $`k_BT=0.02E_h`$, $`E(T)`$ clearly deviates from $`E(0)`$ (thin line with stars) and the approximation $`E(0)=\frac{1}{2}(E(T)+F(T))`$ should be applied. This works very well when comparing the corresponding results (thin line with crosses) with results calculated at $`k_BT=0.001`$ $`E_h`$ (thick line without additional symbol). The surface energy obtained this way is only slightly higher than that from a calculation at $`k_BT=0.001`$ $`E_h`$ which is consistent with figure 1.

A higher smearing temperature also reduces the oscillations in the curves both when using equation 1 or 2. It even allows to substantially reduce the number of sampling points as shown in Figure 5. At low temperature, a high number of sampling points is necessary to obtain the correct result whereas at high temperature already a shrinking factor of 4 (resulting in 6 sampling points in the irreducible zone) is sufficient. It should be noted that in this case we used Equation 1 so that all the data is consistently extracted from calculations on slabs.

In conclusion, calculations at zero temperature are very cumbersome when results from calculations on slabs and on the bulk have to be combined. The error cancellation can be maximized by extracting the surface energy from calculations on slabs only. The problem of error cancellation between bulk and slab is already improved at very low temperature when it is possible to fully converge the bulk energy with respect to reciprocal lattice sampling. Higher temperature also leads to a smoother behaviour of the surface energy as a function of number of layers. Finally, calculations at high temperature can be performed with a strong reduction of the number of sampling points and using $`E(0)=\frac{1}{2}(E(T)+F(T))`$ as an approximation for zero temperature results as suggested in Ref. .

In table V, LDA and PWGGA results for the unrelaxed (100), (110) and (111) surface are summarized. The lattice constant was chosen as the bulk equilibrium lattice constant, a temperature of 0.001 $`E_h`$ and a higher shrinking factor of 24 resulting in 91 sampling points for the slabs and 413 sampling points for the bulk was used.

As found in Ref.  for the jellium model without additional long-range corrections, surface energies in PWGGA are lower than in LDA (the different lattice constants for PWGGA and LDA are not the reason, an evaluation at the LDA equilibrium lattice constant leads to a change of the PWGGA surface energy which is negligible compared to the difference between LDA and PWGGA surface energies). The energies are reduced when going from the smaller $`[3s2p]`$ to the $`[4s3p]`$ basis set, introducing a $`d`$ function leads only to minor changes. Our results agree well with the literature .

## V Conclusion

In this study of lithium metal, we presented accurate results using a full potential, all electron density functional scheme. Results for cohesive properties, elastic constants, band structure, and surface energies are in full agreement with experiment and calculated values from literature. The results in best agreement with experiment were obtained with the gradient corrected functional of Perdew and Wang. Hartree-Fock calculations for lithium are very difficult as it is impossible to optimize the exponents because of the very long range of the exchange interaction; the same problems appear in functionals involving an admixture of Fock exchange. We demonstrated the convergence of the different properties with respect to the computational parameters by using an hierarchy of basis sets, different reciprocal lattice samplings and different smearing temperatures. We showed that finite temperature calculations can be used to improve convergence and still an extrapolation to zero temperature is possible and accurate. Quantities like cohesive energy and lattice constant are already stable with the smallest basis set; elastic constants and surface energies are more sensitive. The most difficult quantity to calculate is the energy splitting between the different phases where we have reached the limit of numerical accuracy. We can not make a prediction about the preferred crystal structure; the energy difference is so small that, from the computational point of view, even subtle changes such as introducing a $`d`$-function are important and zero point energies must be included . We confirm the finding of Ref.  that an approach based on Gaussian type functions provides a reliable and very efficient description of metallic systems.