Structural and vibrational properties of lithium under ambient conditions within density functional theory

TL;DR

This study uses advanced first-principles calculations to explore the phase diagram and vibrational properties of lithium at ambient conditions, highlighting the importance of anharmonic effects and entropy in phase transitions.

Contribution

It introduces a quasi-anharmonic approach to map lithium's potential energy surface and assess phase stability, incorporating quantum nuclear effects and anharmonic vibrations.

Findings

Anharmonic effects significantly influence lithium's phase free energies.

The phase transition sequence aligns with previous studies, confirming quasi-harmonic sufficiency.

Entropy from phonons drives the martensitic transition from close-packed to BCC lithium.

Abstract

We apply a general first-principles approach to derive the phase diagram of metallic Lithium at ambient pressure between 0 and 350 K, including identification of candidate phases. We use ab initio random structure searching (AIRSS) to identify competing phases and supplement the results with calculations of vibrational properties and relevant derived neutron diffraction patterns. Strong quantum nuclear effects are present, prompting a careful treatment of vibrations. We directly map the Born-Oppenheimer surface of Li for the first time, allowing the extension of the normal quasi-harmonic treatment of vibrations to a "quasi-anharmonic" approach, where the effects of anharmonicity are included. The Gibbs free energies of the FCC, BCC, HCP and 9R phases are derived using a variety of equations of state. We find that the anharmonic contribution to the Gibbs free energies of Li phases is of the same order as the differences between phases. Anharmonicity also makes a noticeable difference to the 0 K phonon dispersion of the BCC phase, with the largest difference at the N-point phonon. The ordering of phase transitions that we find agrees with the calculations of Ackland et al. [Science 356, 1254-1259, (2017)], even when anharmonic effects are included, suggesting that a quasi-harmonic treatment is sufficient for correct phase behaviour. We show explicitly that the Martensitic phase transition from close-packed to BCC lithium upon heating is driven by entropic contributions to the phonon free energy.

Tables

Table 1. Table 1: The lowest-energy results of a simple Li AIRSS search. For each space group the number of times it was found ( N 𝑁 N ) is shown. We also report, for each space group, the lowest Helmholtz free energy found ( F m ​ i ​ n subscript 𝐹 𝑚 𝑖 𝑛 F_{min} , relative to the FCC phase) and the corresponding volume ( V 𝑉 V ). The space groups are reported in order of increasing F m ​ i ​ n subscript 𝐹 𝑚 𝑖 𝑛 F_{min} . Only electronic energies are calculated, vibrational energies are neglected in these calculations. The conventional names for the four phases which are investigated more deeply in this work are shown in brackets next to the space group. In total 469 structures were generated sup . It is interesting to note the absence of the DHCP structure (see supplement sup ).

Space group	$F_{min}$ (meV/atom)	$V$ (Å3/atom)	$N$
$Fm\overline{3}m$ (FCC)	0.0	18.96827	39
$P{6}_3/mmc$ (HCP)	0.8115	18.97343	17
$Im\overline{3}m$ (BCC)	1.3684	18.94680	26
$R\overline{3}m$ (9R)	2.0593	18.94582	3
$I4/mmm$	2.1320	18.94480	30
$C2/m$	2.5273	18.97458	10
$Cmcm$	2.7300	18.93219	5
$Immm$	5.2318	18.97877	2
$P{2}_1/m$	11.781	18.99312	3
$P6/mmm$	14.167	19.05003	3

Equations

$$H^{(2)}=\sum _i-\frac{1}{2m_i}{\nabla }_i^2+\frac{1}{2}\sum _{i,j}{\delta }_i{\delta }_j\frac{{\partial }^{2}V(x)}{\partial x_i\partial x_j}.$$

$$H^{(2)}=\sum _{\mathbf{q},\sigma }-\frac{1}{2}\frac{{\partial }^2}{\partial p_{\mathbf{q},\sigma }^2}+\frac{1}{2}{\omega }_{\mathbf{q},\sigma }^2{p}_{\mathbf{q},\sigma }^2.$$

$$H=\left(\sum _{\mathbf{q},\sigma }-\frac{1}{2}\frac{{\partial }^2}{\partial p_{\mathbf{q},\sigma }^2}\right)+V(\{p_{\mathbf{q},\sigma }\},\beta ).$$

$$H_{\sigma ,\nu }(\mathbf{q})=⟨\mathbf{q},\sigma |H|\mathbf{q},\nu ⟩$$

$$G(T,P)=\underset{V}{\min }[F(T,V)+PV]$$

$$I(2\theta )=a{I}_{\text{BCC}}(2\theta )+b\sum _i{I}_i(2\theta )\frac{w_i}{\exp \left(\frac{E_i-\mu }{k_{B}T}\right)-1}$$

$$\begin{array}{rl}{F}_{ph}(T,V)=\frac{1}{2}& \int \omega D(\omega ,V)d\omega \\ +& \int \frac{\omega }{\exp (\omega /T)-1}D(\omega )d\omega -T{S}_{ph}(T,V)\end{array}$$

$$\begin{array}{rl}{S}_{ph}(T,V)=& \int \frac{\omega /T}{\exp (\omega /T)-1}D(\omega ,V)d\omega \\ -& \int \ln [1-\exp (-\omega /T)]D(\omega ,V)d\omega \end{array}$$

Full text

I Introduction

Lithium is the lightest simple metal in the periodic table. It has two stable isotopes, Li-7 (which we focus on in this work) and the much rarer Li-6, both of which exhibit a Martensitic (diffusionless) phase transition from BCC to a close-packed mixed phase upon cooling Ackland et al. (2017); Schaeffer et al. (2015); Deemyad and Zhang (2018); Pichl et al. (2003). Such phase transitions are of significant practical importance, for example in the martensite-austenite transition in steel, upon which many thermomechanical treatments and alloying methods are based Dmitrieva et al. (2011). Martensitic transformations are also crucial to phenomena such as shape-memory of materials James and Hane (2000). Li is the simplest metal that undergoes such a transformation and, because it is such a light element, quantum-nuclear effects make a signifigant contribution to the free energies Ashcroft (1989); Staikov et al. (1997); Deemyad and Zhang (2018). In particular, at low temperatures Li exhibits close-packed polytypism, whereby the Martensitic transition results in a mixture of metastable close-packed phases that can be described by hexagonal layers. The formation of such a mixed close-packed Martensitic phase means that Li is a useful model system for understanding transformation kinetics and characterization of mixed phases, especially where vibrational effects are important Ackland et al. (2017); Caspersen and Carter (2005); Ackland and Jones (2006). Here we demonstrate a general framework for deriving the phase behaviour of such systems.

One of the most challenging aspects of predicting the structures and stability of Li phases is the small energy differences between candidate phases, which are of the order of a few meV per atom Ackland et al. (2017); Staikov et al. (1997). This is because the electronic structure is free-electron like and therefore relatively insensitive to the nuclear positions Naumov et al. (2015). The computational results obtained may therefore be extremely sensitive to numerical accuracy and convergence parameters. We have calculated well-converged energies for the most important structures of Li, such as the BCC, FCC, HCP and 9R phases. These energies include vibrational effects, and we report both quasi-harmonic and fully anharmonic results.

DFT calculations were performed using the castep Clark et al. (2005) plane-wave pseudopotential code and ultrasoft pseudopotentials. We use the ab initio random structure searching approach (AIRSS Pickard and Needs (2011); Needs and Pickard (2016); Pickard and Needs (2006)) to derive candidate phases and to study a wider range of relevant structures and the resulting energy scales. The AIRSS approach consists of the construction of pseudo-random initial structures with a set of optional constraints on the symmetry, inter-atomic distances, coordination numbers, structural units (based on the bonding that is indicated by chemical considerations) and the dimensionality of the system, etc. These structures are then relaxed using DFT geometry optimization in order to sample candidate low-energy structures.

II Vibrational energies

Vibrational contributions are of crucial importance to the thermodynamics of Li phases Staikov et al. (1997); Ackland et al. (2017); Gorelli et al. (2012). A hierarchy of approximations allow a systematic, fully quantum-mechanical treatment of such systems. The Hamiltonian within the harmonic approximation is

[TABLE]

where $x_{i}$ are the nuclear coordinates, $\delta_{i}$ are the displacements of the nuclei from their equilibrium positions and $V(x)$ is the static lattice energy. Using a suitable change of variables Monserrat et al. (2013), Eq. 1 can be rewritten as

[TABLE]

where $p_{\mathbf{q},\sigma}$ are the phonon coordinates. This Hamiltonian corresponds to a series of non-interacting simple harmonic oscillators with frequencies $\omega_{\mathbf{q},\sigma}$. In order to include anharmonic effects, the harmonic potential in Eq. 2 is replaced by the static lattice energy of a given phonon perturbation, $V(\{p_{\mathbf{q},\sigma}\},\beta)$ Monserrat et al. (2013), leading to

[TABLE]

where $\beta$ is the inverse temperature, which can affect the static lattice energy through entropic effects (Mermin entropy) and smearing of the electronic Fermi surface. We found these effects to be negligible in our calculations, helped by the use of a dense electronic k-point sampling which provides an accurate resolution of the Fermi surface sup .

We calculate an anharmonic correction to the phonon dispersion at each phonon q-point independently by diagonalizing $H$ in the basis of single phonon states at that q-point, $\{\ket{\mathbf{q},\sigma}\}$,

[TABLE]

Clearly there is an infinite number of excited states ($\sigma$’s) at each q-point; in order to evaluate the matrix elements we truncate this set to the first 20 exited states, which is more than sufficient for our purposes. The anharmonic potential $E_{el}(\{p_{\mathbf{q},\sigma}\},\beta)$, is mapped using DFT calculations for a discrete set of 31 amplitudes for each mode. The potential is then interpolated to a large number of points (5000) using a quadratic spline. The integral in Eq. 4 is then carried out over these points to obtain the matrix elements of the Hamiltonian, which is then diagonalized. The resulting eigenvalues are the anharmonic excitation energies at a particular q-point, giving an anharmonic correction to the phonon dispersion. These excitation energies can then be used to construct an anharmonic phonon free energy at any given temperature (including entropic and zero-point effects) Fultz (2010); Baroni et al. (2001). We find that this method is around 50 times more expensive than traditional harmonic phonon calculations.

Application of a Legendre transform to the free energies provides the Gibbs free energy at a given temperature and pressure via the following minimization:

[TABLE]

The thermodynamically most stable phase has the lowest Gibbs free energy at a given temperature and pressure. We perform this minimization by fitting $F(T,V)$ data to a suitable equation of state (see results section III.2). Within this method the effects of thermal expansion and any anharmonic contributions to the equilibrium volume of the system are included.

III Results and discussion

III.1 Structure searching

The FCC, BCC, HCP and 9R structures have all been proposed in the past to explain the experimental data for Li Ackland et al. (2017); Overhauser (1984); Barrett (1956); Barrett and Trautz (1948). We find that these structures can be recovered very quickly from first principles calculations using ab initio random structure searching (AIRSS Pickard and Needs (2011); Needs and Pickard (2016); Pickard and Needs (2006)) using only very simple constraints on the initial structures. In particular we constrain the volume per atom to be within 50% of the known value (based on the density of solid BCC Li) and require that no two atoms are within 1Å of one another (approximately 1/3 of the nearest neighbour distance in BCC Li). The unit cells are generated to contain between 1 and 6 atoms. Once a cell has been generated according to these rules we perform a DFT geometry optimization. This minimizes the electronic free energy and locates a locally stable structure (neglecting expensive-to-calculate vibrational effects). The results of an AIRSS search with these constraints are shown in Table 1. We see that it is possible to obtain a good heuristic understanding of the energy landscape, even when the effects of vibrations are neglected. The electronic energies of these phases are extremely similar to one another, with differences on the order of meV/atom. Several other searches were also carried out with additional constraints on the symmetry of the initial structures sup , but this was not found to be useful in this case.

III.1.1 Close packed structures

In the previous section we saw that the AIRSS searches quickly recovered the FCC, BCC, HCP and 9R structures. Most of the other structures found in the searches consist of various close-packed polytypes. Any repeating sequence of hexagonal layers, where no two adjacent layers are of the same type (the layer types are labeled A, B and C; see Fig. 1), gives a close-packed structure. As a result, there is an infinite number of different close packed polytypes. Many of these are realized in the Martensitic phase of Li that is observed experimentally at ambient pressure and low temperatures; this phase consists of a mixture of close-packed phases and remnants of the high-temperature BCC phase Ackland et al. (2017). This is an example of one of the most common forms of Martensitic transformation (BCC $\rightarrow$ close-packed) Katsnelson et al. (1994), in which anharmonic effects may play an important role Gooding and Krumhansl (1988). Close packed crystal structures with a short repeating sequence of hexagonal layers are more likely to be realized in such a mixed phase, as several repeats are more likely to fit between defects. Layer-layer interactions also diminish rapidly with inter-layer separation Loach and Ackland (2017). This justifies the truncation of the infinite set to smaller sequences that are below a certain repeat length. In particular we investigate the energetics of all of the possible close-packed Li structures with a sequence of 15 or fewer layers. There are 10,922 such structures111This number is $\sim\frac{2}{3}2^{14}$. In general the number of close-packed structures with $\leq n$ hexagonal layers is $\sim\frac{2}{3}2^{n}$. The factor of $2^{n}$ arises from the number of leaves in a tree such as Fig. 1 and the $\frac{2}{3}$ arises from the $\frac{1}{3}$ of such structures which have the same initial and final layer in the sequence (meaning they do not correspond to a close packed structure)., 489 of which are unique, labeling equivalent structures following Ref. Loach and Ackland (2017)222We note a minor error in Loach and Ackland (2017), where it is stated that the number of unique close-packed sequences with up to 10 atomic layers is 43. There are in fact only 38.. For a given unique structure $i$, we call the number of equivalent structures $w_{i}$.

The powder neutron diffraction pattern is simulated for each of these unique structures and the resulting overall mixed-phase pattern is obtained by thermodynamic weighting of the individual patterns. Neglecting defect energies the resulting neutron pattern at temperature $T$ is given by sup :

[TABLE]

where $I_{i}(2\theta)$ is the pattern for the $i^{\text{th}}$ structure, $E_{i}$ is the energy per atom and $w_{i}$ is the multiplicity of the structure (as defined in the previous paragraph). $a$, $b$, $\mu$ and $T$ are fitted to the experimental pattern. The resulting combined patterns are shown in Fig. 2.

In agreement with Ackland et al. the experimental data are well described by a Martensite consisting of a mixture of close-packed phases and BCC remnants Ackland et al. (2017). We give the results for such weighted close-packed patterns both at a fixed volume per atom, and after optimizing the volume of each phase. The latter case is equivalent to allowing density variations across the sample based on the local crystal structure. We see little difference in the resulting patterns, likely because the optimized volumes are all very similar.

Simply weighting the close packed phases thermodynamically in this way reproduces many of the non-BCC characteristics of the experimental pattern, but not as accurately as the large-scale molecular dynamics calculations performed in Ref. Ackland et al. (2017). This is because transformation kinetics are not fully included by a simple thermodynamic weighting of phases, as energy barriers along transition pathways between phases also play an important role Caspersen and Carter (2005). Such kinetic effects can be captured by molecular dynamics simulations using suitable interatomic potentials Ackland et al. (2017); Ko and Jeon (2017). We are investigating the exploration of transition states from first principles in order to better understand these effects.

For completeness, the pattern that results from combining the structures found in our AIRSS search using Eq. 6 is also given in Fig. 2.

III.2 Free energy calculations

Quasi-harmonic calculations of the Gibbs free energy are performed for the FCC and BCC phases of Lithium, using a well-converged parameter set. We use plane-wave DFT with a PBE functional, a plane-wave cutoff of 3 keV and an electronic k-point grid with a spacing of 0.02 Å*-1*. The phonon Brillouin zone is sampled using the highly-efficient non-diagonal supercells method Lloyd-Williams and Monserrat (2015) which is significantly faster than a normal phonon supercell calculation, without any loss of accuracy. The phonon dispersion is calculated on a 10x10x10 grid in reciprocal space and then interpolated to a 40x40x40 grid using Fourier interpolation. The Helmholtz free energy is calculated for a range of volumes around the equilibrium volume and fitted to the Birch-Murnaghan equation of state Birch (1947) to extract the Gibbs free energy. An example of such a fit at 300 K is shown in Fig. 5. The resulting Gibbs free energy is shown in Fig. 3 for a range of temperatures between 0 and 350 K. The Murnaghan, Rose-Vinet and Poirier-Tarantola equations of state Murnaghan (1944); Vinet et al. (1987); Poirier and Tarantola (1998) give essentially identical results. A transition from the FCC to a BCC phase occurs upon heating at 217$\pm$13 K. This phase transition is observed on isobaric heating of the FCC phase experimentally; however, the calculated transition temperature is somewhat above the experimental range of 110-200 K Ackland et al. (2017); Pichl et al. (2003); Smith (1987); McCarthy et al. (1980); Schwarz and Blaschko (1990); Crisp (1991). The reverse BCC $\rightarrow$ FCC transition is not seen experimentally upon isobaric cooling at ambient pressure. The FCC phase is instead prepared via a high-pressure route Ackland et al. (2017) to avoid formation of the Martensitic phase investigated in section III.1.1. Isothermal compressibilities and thermal expansion coefficients closely match experimental results sup ; Kittel (2005); Cohen et al. (2003).

Using the method outlined in section II, we calculate an anharmonic correction to the Helmholtz free energy, effectively moving from the quasi-harmonic regime into the “quasi-anharmonic”. The resulting correction to the Gibbs free energy is shown in Fig. 6. The effect on the FCC $\rightarrow$ BCC transition is included in Fig. 3. Anharmonic effects are much stronger in the BCC phase, with the FCC $\rightarrow$ BCC transition temperature increasing by 15 K (to 232 K) as a result.

It is also interesting to note the difference between the lowest-lying eigenvalues of the harmonic and anharmonic Hamiltonians (Eq. 2 and 4 respectively), shown in Figs. 10 and 11. These corrections are present even at zero temperature as they result directly from the anharmonicity of the 0 K Born-Oppenheimer surface. The most significant effect is a stiffening of the N-point $(0,0.5,0)$ phonon in BCC Li by $\sim$ 1 meV (around 30 % stiffer than the harmonic case). It has been suggested that stiffening of this mode could be important in stabilizing BCC phases Chen et al. (1985), however we find that this mode is dynamically stable to begin with (contrary to BCC Zr in Pinsook and Ackland (1999); Chen et al. (1985)).

To complement these results we perform a second set of similar calculations which include the 9R and HCP structures. These calculations are carried out with an LDA functional and with less intensive convergence parameters (1 keV plane-wave cutoff, 0.035 Å*-1* electronic k-point grid spacing and an 8x8x8 phonon q-point grid). This is done in order to investigate the sensitivity of our results with respect to convergence parameters and the functional used. The resulting Gibbs free energies are shown in Fig. 7. The sequence of phase transitions with increasing temperature is the same as in the PBE calculations in Fig. 3. However, the FCC $\rightarrow$ BCC transition temperature in the LDA calculations is 47 K higher at 264$\pm$26 K. This is due to the fact that the transition temperature is extremely sensitive to the convergence of the differences in the Gibbs free energy; we estimate an increase of 117 K in transition temperature per meV increase in the BCC-FCC energy difference. The appearance of a metastable BCC $\rightarrow$ 9R transition upon cooling through 120 K is compatible with the Martensitic transition, however the predicted transition temperature is once again somewhat above the experimental range Smith (1987); Schwarz and Blaschko (1990); Smith et al. (1990); McCarthy et al. (1980); Pichl et al. (2003); Crisp (1991).

We go on to show that the FCC $\rightarrow$ BCC transition is driven by entropic contributions to the phonon free energy. The phonon free energy is given by:

[TABLE]

where $D(\omega)$ is the phonon density of states. The first term is the zero-point energy, the second arises from the thermal occupation of phonon states and the entropic term, $S_{ph}$, is given by Fultz (2010); Baroni et al. (2001):

[TABLE]

The importance of the entropic term becomes apparent when its contribution is explicitly neglected, resulting in the Gibbs free energy landscape shown in Fig. 8 in which phase transitions are no longer present. The individual contributions to the phonon free energy for BCC Li are shown in Fig. 9, in which clearly that the entropic contribution dominates the zero-point and occupational terms. The FCC, HCP and 9R phases show similar behaviour. Because the low-lying phonon modes in BCC Li are softer than those of FCC (see Figs. 10 and 11), the entropic effects are stronger, which leads to the FCC $\rightarrow$ BCC phase transition upon heating.

IV Conclusions

We have investigated the application of a general framework for deriving the phase behaviour of materials with strong vibrations from first principles to metallic Li. The AIRSS approach Pickard and Needs (2011); Needs and Pickard (2016) is found to quickly pick out relevant low-energy phases with little effort, including both the BCC phase and numerous close-packed phases. We have found that DFT predicts a zero-temperature FCC phase and a room temperature BCC phase; a metastable transition from BCC to the 9R phase is also predicted in the region of FCC stability. These results are consistent with experimental results, which show a (partial) Martensitic transition from BCC to various close-packed forms on cooling Ackland et al. (2017); Smith (1987); McCarthy et al. (1980); Schwarz and Blaschko (1990). Experiments also show that the FCC phase is stable at low temperatures, but may be formed via a high-pressure pathway Ackland et al. (2017). The predicted FCC $\rightarrow$ BCC transition temperature is found to be extremely sensitive to changes in the Gibbs free energy landscape ($\sim$ 117 K/meV). As a result, effects which would normally be considered negligible become potentially important, including the effects of anharmonic vibrations. We calculate an anharmonic correction of up to 0.4 meV/atom, which increases the predicted FCC $\rightarrow$ BCC transition temperature by 15 K to 232 K, but find that it does not qualitatively change the phase diagram. Our calculations show that the N-point phonon in BCC Li is dynamically stable, contrary to the case for BCC Zr in Ref. Pinsook and Ackland (1999); Chen et al. (1985) and is around 30% stiffer when anharmonic effects are included. The anharmonic free energy calculations are around 50 times more expensive than the quasi-harmonic ones, but we have found that they can be accelerated without loss of accuracy by employing symmetry considerations and sampling anharmonic potentials at fewer amplitudes. We calculate the various contributions to the Gibbs free energy and find that entropic contributions to the phonon free energy are of crucial importance, without which such phase transitions disappear entirely.

Acknowledgements.

M.H. would like to acknowledge the EPSRC Centre for Doctoral Training in Computational Methods for Materials Science for funding under grant number EP/L015552/1. We thank Bartomeu Monserrat, Joseph Nelson and Nick Woods for helpful discussions. R. J. Needs acknowledges financial support from EPSRC under Critical Mass Grant EP/P034616/1 and the Support for the UKCP consortium grant EP/P022596/1. This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (http://www.csd3.cam.ac.uk/), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk). Supporting research data may be freely accessed at https://doi.org/10.17863/CAM.31655, in compliance with the applicable Open Data policies.