Influence of the exchange-correlation functional on the quasi-harmonic lattice dynamics of II-VI semiconductors

First-principles computation is well established as a powerful tool in materials science. The most widely used theoretical “workhorse” at present is the Kohn-Sham density-functional theory (DFT) formalism,^1,2 which recasts the problem of an interacting many-electron system as a system of non-interacting particles experiencing an effective potential. The ability to perform atomistic quantum-mechanical calculations at a manageable computational cost not only facilitates the use of theory to interpret experimental measurements but also the characterization of new materials in silico, a capability which has enabled ambitious undertakings such as the Materials Project.³ 

The key technical parameter in DFT calculations, leaving aside implementation details such as the basis set used to expand the Kohn-Sham wavefunctions, is the functional form used to calculate the exchange-correlation (XC) energy for a given spatial electron density, n(r). Following Perdew and co-workers,⁴ XC functionals are typically classified according to a “Jacob’s ladder” of approximations, with each rung in the hierarchy representing a better description of the true many-electron physics.

The simplest functionals are based around the local (spin)-density approximation (L(S)DA), in which the exchange and correlation energies of the electron density at a point in space are approximated by that of a homogenous electron gas with the same density. For some systems, this approximation benefits from a fortuitous cancellation of errors, but for most, it does not yield satisfactory results.

The LDA can be improved upon by including the density gradient, ∇n(r), which is the basis of the semi-local generalised-gradient approximation (GGA) functionals. GGA functionals tend to yield improved energetics and to lengthen bonds, thereby increasing cell volumes and lattice constants and softening phonon frequencies. GGA functionals thus correct for the phenomenological over binding exhibited by the LDA, but in many cases overcompensate and under bind, overestimating lattice constants and underestimating phonon frequencies.⁵ Perhaps the most widely used GGA at present is the ubiquitous Perdew-Burke-Enzerhof (PBE) functional,⁶ which was devised as an improvement to the earlier PW91 GGA, with a simpler functional form and derivation. Other well-known GGA functionals include the two revised PBE functionals of Hammer et al. (revPBE and RPBE),⁷ which can yield improved energetics compared to PBE, AM05,^8,9 which was designed to better treat surface effects, and PBEsol, a variant of PBE optimised for solids.⁵ It has been shown^5,6 that GGAs suffer from a fundamental limitation, in that a larger dependence on the density gradient yields better atomisation energies, but poorer lattice parameters and surface energies, and vice versa, through differences in the description of the exchange energy. As a result, GGAs optimised specifically for solids can be constructed with a reduced density dependence, which is the approach adopted in PBEsol.⁵ 

Building on the GGA, meta-GGA functionals include the second derivative of the density, typically via the kinetic-energy density τ(r). This family of functionals include the original Perdew-Kurth-Zupan-Blaha (PKZB) functional¹⁰ along with the well-known Tao-Perdew-Staroverov-Scuseria (TPSS)¹¹ model and its subsequent revision,¹² and the fitted M06-L functional developed by Zhao and Truhlar.¹³ In principle, the greater flexibility of the meta-GGA formalism allows TPSS to describe both solids and molecules with good accuracy, overcoming the trade-off inherent in the GGA,¹¹ although applying the same principle used to construct PBEsol to TPSS, which is the spirit of revTPSS, was found to improve its description of solids.¹² 

Hybrid functionals such as PBE0¹⁴ and the HSE family,^15–17 which improve on (meta-)GGAs by including a fraction of the non-local Hartree-Fock exact exchange energy, are becoming increasingly affordable with modern computational resources. Hybrids are routinely used for accurate electronic-structure calculations, but at present are impractically expensive for tasks such as structural optimisation, at least in the general case.

One well-documented limitation of (semi-)local DFT functionals is their poor description of the non-local electron correlation which gives rise to dispersion (van der Waals) forces. Since a theoretical treatment of these weak interactions can be very challenging, several approximate approaches to introduce dispersion corrections into (meta-)GGAs have been developed.^18–25 The vdW-DF method of Lundqvist et al.²³ and its derivatives^19,20,22 uses a parameterisation of the local electron density and its gradient to correct the LDA correlation energy, with the exchange energy being calculated using an underlying GGA functional. However, while these functionals perform well for systems where dispersion interactions are significant, when used in calculations on bulk solids, they tend to exhibit the same issues as the GGAs they are based on.¹⁹ 

The Grimme DFT-D functionals adopt a semi-empirical approach, adding damped dispersion corrections to the energies and forces obtained in DFT calculations with a scaling factor adapted to the XC functional being used. In the popular DFT-D2 scheme,²⁴ the corrections are based on pairwise London interactions parameterised by atomic ionisation potentials and static dipole polarizabilities. The newer DFT-D3 scheme²¹ improves on this in several ways, including accounting for the change in atomic polarizability as a result of bond formation, allowing the dispersion correction to be varied based on the atomic coordination environment, and including three-body as well as pairwise interactions. An alternative scheme proposed by Tkatchenko and Scheffler²⁵ uses a similar principle to DFT-D2, but with charge-density dependent dispersion coefficients and damping functions, and a revised version of this DFT-TS scheme accounts self-consistently for the screening effects of neighbouring ions.¹⁸ 

Most DFT-based material studies have focussed on investigating the “static” properties of athermal structures, such as the equilibrium lattice parameters and bulk moduli, the electronic structure and energy gaps, and the optical properties. However, an increasing number of studies aim to take into account the effect of temperature, for example, through the use of ab initio molecular-dynamics (MD)^26,27 and lattice-dynamics calculations.^28–30 The latter provide an easily interpretable view of the microscopic motion of atoms in solids³¹ and can also provide thermodynamic information such as free energies³² and heat capacities,^28,29 which can in turn be used to predict finite-temperature properties.^30,32 Lattice-dynamics calculations can also be used to model and interpret data from a range of routine experimental measurements with near-quantitative accuracy,^33,34 and thus can provide references to, for example, aid the use of quantitative spectroscopic techniques as a tool for characterising structural defects.³⁴ 

Whereas the performance of different DFT functionals in reproducing static properties has been well characterised,^35–37 the ability of common functionals to model lattice dynamics and finite-temperature structure and properties is poorly studied. In principle, lattice-dynamics calculations are a stringent test of functional performance, since they probe not only the predicted equilibrium structure but also the response to atomic displacements. It has also been pointed out that lattice-dynamical effects such as the vibrational zero-point energy and thermal expansion should be taken into account when comparing experimental results to theoretical predictions.^38,39 Building on these previous studies investigating the effect of the zero-point energy on the predicted 0 K structure and properties of bulk materials obtained with different functionals, we have carried out a systematic comparison against experimental data of the finite-temperature properties of four II-VI semiconductors, viz., PbS, PbTe, ZnS, and ZnTe, obtained from quasi-harmonic lattice-dynamics calculations with eight popular exchange-correlation functionals. Our tests include the LDA, the PW91,^40,41 PBE⁶ and PBEsol⁵ GGAs, the TPSS¹¹ and revTPSS¹² meta-GGAs, and the semi-empirical PBE-D2 and D3 functionals,^21,24 all of which are well established in the theoretical literature.

The rocksalt and zinc-blende structures of the lead and zinc chalcogenides, respectively, represent common high-symmetry geometries adopted by a large number of binary systems, and experimental data for comparison to the calculations is readily available in data handbooks.^42–44 Despite their simple structures, however, the four materials display some interesting dynamical properties. PbTe has been widely studied as a candidate thermoelectric material,^45,46 and its thermoelectric efficiency has been attributed in part to its low thermal conductivity, which arises from unusual features in its lattice dynamics.^47–49 Among these is a possible low-temperature phase transition from the rocksalt to a lower-symmetry structure,⁴⁷ which is proposed also to occur in PbS. However, the literature is currently divided on this phenomenon, with different measurement techniques giving conflicting results,^47,50,51 and the conclusions of the original experiments recently being disputed.⁵² Both zinc chalcogenides exhibit a small negative thermal expansion at low temperature,^53,54 a subtle feature that could potentially be difficult to model.

Our work is organised as follows: Section II presents an overview of the quasi-harmonic lattice-dynamics methodology, and Section III gives the technical details of the calculations performed in this study. In Section IV, we consider the equilibrium structure and harmonic lattice dynamics of the four materials, and in Section V we compare the temperature-dependent material properties obtained from the quasi-harmonic calculations. The results are discussed in Section VI, with particular focus on the ability of the functionals to quantitatively reproduce structural and vibrational properties at 300 K. We then close with some concluding remarks in Section VII.

The theory of lattice dynamics provides a framework for modelling the phonons in periodic solids.⁵⁵ A crystal of N atoms is modelled by sets of 3N independent quantum-harmonic oscillators with associated reciprocal-space wavevectors q, which define the wavelengths and propagation directions of the atomic-displacement waves.

The central quantities required for a lattice-dynamics calculation are the interatomic force-constant (IFC) matrices, Φ_αβ, which capture the changes in force that arise in response to the displacement of atoms from their equilibrium positions,

The subscripts α and β denote the Cartesian directions, while l and l′ are the unit cells in which atoms i and j, respectively, reside relative to one another. The IFCs can be computed from numerical differentiation using the “direct” finite-displacement approach, or from a perturbation-theory technique such as density-functional perturbation theory (DFPT). In the former approach, the Parlinski-Li-Kawazoe supercell method⁵⁶ is commonly used to include IFCs whose range extends beyond a single unit cell, i.e., q-vectors away from the centre of the Brillouin zone.

The phonon frequencies and eigenvectors are then obtained from diagonalization of a 3N × 3N dynamical matrix, D_αβ, which can be constructed from the force-constant matrices, for a given q-vector, according to

where m_i is the mass of atom i and r(il) is the position of atom i in the lth unit cell.

In addition to providing access to the phonon density of states (DOS) and band dispersions in q-space, the phonon frequencies integrated over the Brillouin zone can also be used to calculate the vibrational contributions to the constant-volume (Helmholtz) thermodynamic free energy, A, according to

where U_L and U_V are the lattice and vibrational internal energies, respectively, and S_V is the vibrational entropy. We note that U_L is the only quantity available from standard total-energy calculations, and thus calculating the temperature dependence of the free energy requires the contribution from the lattice dynamics to be taken into account.

In practice, A is typically obtained via the bridge relation from statistical mechanics, $A ( T ) =− k B TlnZ T$⁠, where the partition function, $Z T$⁠, is defined in this model as follows:

In this expression, the phonon frequencies ω are indexed by a q-vector and a phonon (band) index, v.

The quasi-harmonic approximation (QHA) extends the harmonic approximation to account for anharmonic effects due the variation in lattice volume with temperature by allowing the phonon frequencies, and hence the free energy, to become volume dependent. In the QHA, harmonic-phonon calculations are performed at a range of expansions and compressions about the (athermal) equilibrium volume, and the Gibbs free energy is obtained at a target temperature T and pressure p from the expression

where the notation $min V [ expr ]$ indicates that, for each value of T and p, expr is minimised with respect to the volume, V, most commonly by fitting the free energy/volume curve to an equation of state such as the Murnaghan or Vinte-Rose expressions.^57,58 The implicit assumption is that the phonons remain harmonic at each temperature, an approximation which, as a “rule of thumb,” is valid up to around 2/3 of the melting temperature, T_m, before higher-order anharmonic contributions dominate.^59,60

From the temperature-dependent Gibbs energies, volumes and bulk moduli obtained from the equation-of-state fits, the temperature dependence of a number of derived properties can hence be calculated, e.g., expansion coefficients, constant-pressure heat capacities, and the (average) Grüneisen parameter. Under the assumption that the main effect of temperature is captured by the volume of the lattice, the QHA in principle allows the temperature dependence of a number of other important material properties to be computed from first principles, making it a powerful tool for the characterisation of bulk materials.^30,32

Quasi-harmonic lattice-dynamics calculations were performed according to the procedure used in our previous work.^30,32 The Phonopy package⁶¹ was used for the setup and post-processing of supercell finite-displacement phonon calculations, with the Vienna ab initio Simulation Package (VASP) code⁶² used as the force calculator.

We carried out calculations with eight well-known DFT functionals implemented in VASP. As a baseline, we used the Perdew-Zunger parameterisation of the local-density approximation,⁶³ which we compared against the PW91,^40,41 PBE,⁶ and PBEsol⁵ GGA functionals. We also performed calculations with the TPSS¹¹ and revTPSS¹² meta-GGA functionals, and the dispersion-corrected PBE-D2 and D3 functionals.^21,24 The elemental C₆ and R₀ PBE-D2 parameters for S, Zn, and Te used were drawn from the VASP internal database, and the value for Pb was taken from Ref. 64. The values of the other parameters for PBE-D2 and D3 were left at the VASP defaults for use with the PBE functional; we note that the PBE-D3 calculations were performed with the original formulation of the method, rather than the variant which uses the Becke-Jonson damping scheme.⁶⁵ 

For consistency, we kept other key technical parameters as similar as possible. A kinetic-energy cutoff of 550 eV for the plane-wave basis was adopted in all calculations. For the calculations on the two-atom primitive cells, Γ-centred Monkhorst-Pack k-point meshes with 8 × 8 × 8 subdivisions were used to sample the Brillouin zone, with the grid being reduced proportionally for the supercell finite-displacement calculations. Convergence testing (see Section III of the supplementary material¹⁰²) showed that these settings were generally sufficient to converge the absolute total energies to within 1 meV per atom, and the stress tensor to within 1 kbar (0.1 GPa). All simulations were performed with the LDA projector augmented-wave (PAW)^66,67 pseudopotential set distributed with VASP, treating the outermost s and p electrons in Pb, S, and Te, and the outermost s, p, and d electrons in Zn, as valence states. The PAW projection operators were applied in reciprocal space, and the precision of the charge-density grid was chosen automatically to be sufficient to avoid aliasing errors. Electronic wavefunctions were optimised to a tolerance of 10⁻⁸ eV, and structural optimisations were set to be converged when the magnitude of all forces was less than 10⁻² eV Å⁻¹, although in practice, the relaxations were invariably constrained by the high symmetry of the structures.

For the LDA and GGA functionals, the Born effective-charge tensors and the electronic-polarisation components of the macroscopic dielectric constants, needed to compute the non-analytical corrections to the phonon frequencies to account for long-range electrostatic interactions (longitudinal-optic (LO)/transverse optic (TO) splitting), were obtained using the DFPT routines in VASP.⁶⁸ This functionality is currently not available for the meta-GGA and PBE-D2/D3 functionals, so in these calculations the quantities were computed by application of a finite electric field (10⁻³ eV Å⁻¹) according to the PEAD approach of Nunes and Gonze.⁶⁹ This technique requires that the systems are insulating, and the small gap of PbS, when combined with the usual issue of the Kohn-Sham bandgap obtained from (semi-)local functionals being underestimated with respect to the electronic gap,⁷⁰ meant that it was not possible to compute the corrections for these systems with these functionals. It has previously been found that convergence of the Born effective-charge tensors and dielectric constants requires a denser k-point sampling than is needed for electronic and force convergence,^30,71 and so we employed Γ-centred meshes with 16 × 16 × 16 and 12 × 12 × 12 subdivisions in the DFPT and finite-field calculations, respectively, with the reduced sampling used in the latter due to the higher computational demand of these calculations. When testing the two procedures with ZnS and PBE (Table S3 in the supplementary material¹⁰²), we found that they gave appreciably different results, perhaps due to the different convergence criteria employed. However, as LO/TO splitting has only a very minor effect on the phonon density of states, and hence on the finite-temperature properties computed from the quasi-harmonic calculations, we did not investigate this in detail. We also avoid using the Γ-point LO frequencies, which depend strongly on the splitting correction, to compare functionals in our analysis.

Force calculations were performed on 4 × 4 × 4 expansions of the primitive cells. For the quasi-harmonic calculations, harmonic-phonon calculations on a set of eleven expansions and contractions about the equilibrium volume (∼ ± 5%) were used in the free-energy equation-of-state fits, which were performed using the Vinet-Rose functional form.⁵⁸ For calculating the phonon DOS curves and partition functions, the Brillouin zone was integrated by sampling with a 48 × 48 × 48 Γ-centred q-point mesh.

Finally, as discussed in Section IV, and also Section I of the supplementary material,¹⁰² for the Born-charge and supercell-force calculations on the zinc chalcogenides with the meta-GGA functionals, we found it necessary to increase the number of points in the charge-density grid by a factor of 1.5× along each dimension from the “High” default in VASP.

The starting point for the discussion is the athermal equilibrium primitive-cell volumes and bulk moduli obtained from fits of the energy/volume curves computed with the eight functionals to the Murnaghan equation of state⁵⁷ (Fig. 1). In keeping with the phenomenological over/under-binding characteristics widely reported in the literature,^{35–37,72–74} the LDA consistently predicts the smallest equilibrium volume for all four materials, whereas PW91 and PBE, which yield very similar results, invariably predict the largest; the expansion predicted by the latter two generally amounts to a 5%-10% increase in volume over the LDA. The volumes obtained from PBEsol are intermediate between the two extremes. Both meta-GGAs yield smaller volumes than PBE, with revTPSS giving the lower of the two, but both predict equilibrium volumes closer to the PBE values than to the LDA results. Adding the Grimme D2 correction to PBE appears consistently to give similar results to PBEsol, whereas the PBE-D3 results tend to be more similar to PBE/PW91 than to PBE-D2, at least in the case of the lead chalcogenides.

The smaller equilibrium volumes lead to larger bulk moduli, which is mainly due to the factor of 1/V in the definition—for a given material, the equations of state calculated with different functionals have very similar curvature when shifted to the corresponding equilibrium volumes. The variation in the moduli obtained from the eight functionals is more pronounced than that in the equilibrium volume, with PW91/PBE underestimating the LDA values by 10%-20%.

These comparisons provide a quantitative measure of the variation in equilibrium structure that can be expected from studies performed with different functionals, at least for the simple bulk materials studied in the present work, and reflects the variation in values frequently observed between different modelling studies (e.g., for PbTe^{28–30,78–80}).

Comparing the predicted equilibrium volumes and lattice constants to low-temperature measurements, the best-performing functionals are broadly those optimised for solids (i.e., the PBE GGA and the revTPSS meta-GGA), with the former yielding better predicted bulk moduli. For some systems, the dispersion-corrected functionals also yield results close to the experimental values, although not as consistently. However, it is important to note that these comparisons do not take into account dynamical effects (i.e., zero-point energy corrections);^38,39 this point is considered further in Sec. V.

Within the QHA, the key quantity is the volume dependence of the phonon frequencies.^30,32 It is therefore of interest to compare the phonon frequencies predicted by the eight functionals at the calculated equilibrium volumes, again to see what sort of variation can be expected. Fig. 2 shows the calculated equilibrium phonon dispersions and DOS curves for PbS and PbTe; corresponding plots for the zinc chalcogenides may be found in Fig. S12 in the supplementary material.¹⁰² 

The shapes of the dispersions and DOS curves are broadly similar between functionals, but with some bands varying in frequency by up to ∼0.5 THz (∼17 cm⁻¹). In keeping with the variation observed in Fig. 1, the phonon frequencies predicted by the LDA and PBE/PW91 generally form upper and lower bounds, respectively, with the other functionals yielding intermediate values. The discrepancy in the phonon dispersions between functionals is most pronounced at the zone centre (Γ); this is due to the inclusion of LO/TO splitting, which effectively incorporates the variation in the electronic-polarisation component of the dielectric tensor and the Born effective charges predicted by the different functionals into the dispersions. We note that, as mentioned in Section III, due to the small bandgap of PbS we were not able to compute the LO/TO splitting for this material with TPSS, revTPSS, PBE-D2 and PBE-D3.

During our calculations on the zinc chalcogenides, we encountered problems with numerical noise in the force constants calculated with the GGA, meta-GGA and dispersion-corrected functionals, which led to visible artefacts in the phonon dispersions and DOS curves, in particular unphysical negative rigid-translation modes at the Γ point. After carefully verifying the technical parameters, we linked the problem to the precision of the charge-density grid, which is used to calculate the density derivatives required by the (meta-)GGA functionals. This is discussed in detail in Section I of the supplementary material.¹⁰² 

We were able to correct the artificial imaginary modes in the GGA dispersions straightforwardly by symmetrising the force constants within Phonopy. This procedure enforces the acoustic-sum rule, and is performed by default in many software packages. However, we were only able to obtain satisfactory results with the meta-GGAs by increasing the size of the charge-density grid by a factor of 1.5× (from the “High” default in the VASP code). Increasing the plane-wave cutoff and using a denser k-point sampling both made no visible difference to the phonon dispersions, and while a supercell expansion of the conventional cell appeared to yield considerably better results, visible artefacts remained, and the subsequent QHA calculations displayed clear anomalies when compared to the results from the other functionals. It is worth noting that increasing the precision of the charge-density grid led only to sub-meV differences in the total energy (Table S1 in the supplementary material¹⁰²) and, given that the geometry of the zinc blende structure is completely constrained by symmetry, would therefore not be expected to influence the equilibrium volumes and bulk moduli in Fig. 1.

A second issue we encountered, detailed in Section II of the supplementary material,¹⁰² is that the PBE-D2 calculations on PbTe yielded obviously unphysical phonon band structures, despite the same elemental C₆ and R₀ parameters for Pb and Te working well for PbS and ZnTe, respectively. Our interpretation of this is that the atomic-polarizability correction used in the D2 method is a poor approximation for systems such as PbTe with heteropolar bonding between highly polarisable ions. We therefore opted also to exclude the PBE-D2 results on PbTe from this study.

We found that PBE-D3, which attempts to correct for interactions between atoms, gave better results, although the calculated phonon density of states exhibited a discontinuous lowering of the optic-mode frequencies under some of the larger lattice expansions (Fig. S11 in the supplementary material¹⁰²), and these volumes had to be excluded from the quasi-harmonic calculations in order to obtain physical results (see Sec. V). This suggests that PbTe apparently represents a problematic system for both the D2 and D3 corrections.

Finally, we noted that VASP does not enforce symmetry in the Born effective-charge tensors. In the four two-atom primitive cells studied here, the effective charge on the anion should be equal and opposite to the charge on the cation. In the present calculations, the Born charges computed with DFPT were invariably close to being symmetric, whereas those computed using the finite-field approach deviated substantially. To test this further, we recalculated the ZnS and ZnTe PBE Born charges with the finite-field approach and found that these also required symmetrisation. We therefore symmetrised the Born charges obtained with the finite-field method, i.e., those calculated with the meta-GGA and DFT-D functionals, using the approach implemented in Phonopy.⁸¹ 

As a basic check on the phonon frequencies and the general computational setup, we verified that the constant-volume heat capacity (C_V) curves calculated with the various material/functional combinations all tend to the high-temperature Dulong-Petit limit of 6R per mole per formula unit (3R per mole per atom, Fig. S13 in the supplementary material¹⁰²).

A feature of note in the equilibrium phonon dispersions is the presence of an imaginary mode at the X point in the PbS/LDA curve (Fig. 2(a)). We found that this anomaly was also present at additional lattice volumes close to the predicted equilibrium, but was absent under larger expansions and compressions (Fig. 3). After some investigation, we found that this was due to the small bandgap of PbS (0.29 eV at 4 K experimentally⁴⁴). Since the bandgap in this material decreases with volume,⁸² and, as noted previously, the Kohn-Sham bandgap obtained from DFT is typically lower than the electronic gap,⁷⁰ the LDA predicts the system to be on the border of a metal-insulator transition at its equilibrium volume. The negative frequency then arises due to deviation from the expected linear-response behaviour during the finite-displacement force calculations, i.e., small ionic displacement leads to a disproportionately large change to the underlying electronic structure and hence the calculated forces.

To explore this further, we plotted the volume dependence of the Kohn-Sham bandgap of PbS with the eight functionals (Fig. 4). All predict similar bandgap deformation, with a shift of around 0.2 eV between the LDA, which appears to predict the smallest gap, and PBE and TPSS, which predict the largest. This qualitative similarity between the functionals is also evident in the calculated volume dependences of the bandgaps in PbTe and ZnS/ZnTe (Fig. S14 in the supplementary material¹⁰²), although the range observed for a given volume depends on the size of the gap. It is also worth noting that the PBE and PBE-D2/D3 bandgaps lie on overlapping curves, reflecting the fact that the dispersion corrections adjust the PBE forces, and hence equilibrium volume, but not the underlying electronic structure. While only the LDA predicts a zero gap at its equilibrium volume, PBEsol and PBE-D2 both predict a crossover from a semiconducting to a metallic state under compression. We would therefore expect the dispersions computed with PBEsol and PBE-D2 at these volumes to show the same phenomenon as the equilibrium LDA band structure, which we indeed observed to be the case (Figs. S15 and S16¹⁰²).

In addition to this artefact, we also observed imaginary modes in the phonon dispersions of the lead chalcogenides calculated with PW91 and PBE (Figs. S9 and S17-S19 in the supplementary material¹⁰²) at large (∼3%-5%) expansions, which we ascribed to the large equilibrium volume predicted by these functionals. Similar phenomena are also observed for both compounds with PBE-D3 (Figs. S11 and S20¹⁰²). While this may be responsible for the optic-mode softening observed in the PbTe phonon density of states noted above, the same problem did not occur in the PbS/PBE-D3 calculations. In all cases, the imaginary modes do not make a significant contribution to the overall phonon DOS (Figs. 3, S9, S11, and S15-S20¹⁰²), and therefore, at least to first approximation, they are not expected to significantly influence the results of the QHA calculations discussed in Sec. V.

As a final point, it is of interest to compare the computational cost of performing supercell-force calculations with the eight different functionals. We thus performed the phonon calculations on PbS presented in Fig. 2 with an identical set of resources (96 cores, made up of six dual-CPU Intel Xeon E5-2650v2 nodes with 128 Gb RAM and an Intel TrueScale QDR Infiniband interconnect), and recorded the total time taken for the calculation, plus the average time required for an electronic SCF cycle (Table S4 in the supplementary material¹⁰²). In these calculations, the GGAs and PBE-D2/D3 took 1.2-1.5× longer than LDA for a complete set of force calculations, while TPSS and revTPSS took 5.17× and 4.6× longer, respectively. It is apparent from comparing the total and SCF-cycle times that the differences in performance are to a large extent due to the calculations taking different numbers of SCF cycles to converge, in particular those with the meta-GGAs; individual SCF cycles with the GGAs and PBE-D2/D3 took a comparable time to those run with LDA, while SCF steps with the meta-GGAs took ∼2.5× longer. It is perhaps also worth noting that the PBE-D2 and PBEsol calculations were slightly faster than those run with PBE and PW91, which is perhaps strange given their (notionally) comparable complexity, but may be due to their being better optimised for use with bulk systems.

In this section, we discuss the functional dependence of the finite-temperature properties obtained from the QHA calculations, and compare the calculations against available experimental data.

As a starting point, it is of interest to quantify the differences in equilibrium structure obtained when including the vibrational zero-point energy in the equation of state fits, i.e., to compare the 0 K cell volume, lattice constants, and bulk moduli obtained from the QHA calculation to the athermal parameters presented in Fig. 1. Across the 31 sets of QHA calculations performed in the present study, we observed an average 0.4% variation in the cell volume, with the minimum and maximum being 0.29% and 0.64%, respectively. This translates to an average variation of 0.13% in the lattice constant. Given its greater sensitivity to the volume, we observed a larger variation in the bulk modulus, with an average difference of −0.95% and a maximum absolute deviation of 13.48%. The latter was observed in the PbTe calculations using the PBE-D3 functional, and a similarly large 7.38% difference was observed in the PbS/PBE-D3 calculation, suggesting that phonon calculations on both lead chalcogenides may represent problematic cases for the DFT-D3 functional. Excluding these outliers, the absolute deviation observed between the athermal and 0 K bulk moduli across the rest of the datasets was a considerably smaller 2.9%. Except for the PBE-D3 calculations on PbTe, we found that including thermal effects invariably increased the 0 K volume and hence the lattice parameter; in most cases, this led to a corresponding decrease in the bulk moduli, although not always. The raw data from these comparisons is presented in Tables S5-S8 in the supplementary material.¹⁰² 

Considering the temperature dependence of the lattice constant (Fig. 5), the same trend as was observed for the equilibrium volume is clearly evident, i.e., that the LDA underestimates the lattice constant compared to experiment, whereas PBE and PW91 both overestimate it. For the lead chalcogenides, the overestimation with the two GGAs is further exaggerated, such that for PbS the ±5% range of volumes used in the QHA calculation was only sufficient for structure predictions up to 390 K. PBE-D3 exhibits a similar issue, and the overestimation of the PbS lattice constant again suggests that this functional has problems describing the dynamical properties of both lead chalcogenides The meta-GGAs improve on PBE/PW91, with revTPSS matching the experimental data better than TPSS, but neither comes as close to experiment as PBEsol, which performs consistently well across all four of the materials studied. In contrast, for the zinc chalcogenides, the temperature dependence of the lattice constants is best captured by revTPSS and PBE-D3.

PBEsol and PBE-D2 also perform quite well for these systems, whereas TPSS appears to overestimate the lattice constants in a similar manner to PW91 and PBE. Whereas PBE-D2 yields similar temperature dependences to PBEsol for the zinc chalcogenides, for PbS, it predicts a result closer to the LDA, i.e., tending to over bind with respect to experiment. Taken together with the issues we had during the calculations on PbTe using this functional, this suggests PBE-D2 may also have problems with describing both lead chalcogenides.

Interestingly, while PBEsol best reproduces the absolute value of the lattice constant, the linear-expansion coefficients $α L = 1 a d a d T$ are consistently reproduced extremely well by the LDA (Fig. S21 in the supplementary material¹⁰²), with PBE-D3 also performing well for ZnS. All the functionals tested successfully capture the subtle low-temperature negative thermal expansion of ZnS and ZnTe, although again the LDA does so most successfully.

This raises the interesting question of whether the differences between the functionals in Fig. 5 are primarily due to differences in the equilibrium volume rather than the thermal expansion. To test this, we replotted Fig. 5 with the curves for the eight functionals shifted to their respective equilibrium volumes (Fig. S22¹⁰²). For the zinc chalcogenides, the main contributor to differences in the temperature dependence of the lattice constant is indeed the athermal starting point, with the difference between the thermal expansion predicted by the functionals being on the order of 10⁻² Å at 500 K. On the other hand, for the lead chalcogenides, which display a considerably larger thermal expansion, the difference between functionals is more of a balance between both the athermal starting point and the predicted expansion, although the former still plays a significant role.

In principle, one could use a more sophisticated functional to evaluate the equation of state and then a cheaper functional (e.g., LDA) to perform the phonon calculations to obtain the temperature dependence of the lattice constant using the QHA. However, the present results suggest that both components can be obtained fairly accurately using a GGA functional optimised for solids, which does not require much more computing power than the LDA. Moreover, for more complex systems with a larger number of degrees of freedom (e.g., internal ion positions or the cell shape), simpler functionals may not capture the volume dependence of the structure correctly.

The temperature dependence of the bulk moduli of the four materials is shown in Fig. 6. As remarked in previous work,³⁰ the temperature dependence of this property can be particularly strong, which implies that, if an athermal calculation gives a good reproduction of a bulk modulus measured at 300 K, and the material undergoes moderate thermal expansion, the functional may not be describing the structure well. For example, for PbTe, the 300 K bulk modulus of 39.8/40 GPa^85,86 is well matched by the equilibrium moduli predicted by PBE and PW91, whereas the 300 K values predicted from the QHA calculations with these functionals are too soft.

For PbTe, the temperature dependence is best reproduced by PBEsol. For PbS, the spread among different 300 K experimental values is a significant fraction of that among the eight functionals, but if these data points are taken to be reliable then the LDA, PBEsol, and PBE-D2 yield the best reproduction of the finite-temperature moduli for this system. For ZnS, PBEsol and the two DFT-D functionals perform very well. For ZnTe, none of the functionals yield values as close to the experimental measurements as for PbTe and ZnS, but the finite-temperature moduli, and the shallow temperature dependence, appear to be best reproduced by PBE-D2/D3 and the LDA, with PBEsol also performing well up to ∼200 K.

Finally, we also calculated the constant-pressure heat capacity (C_p) curves obtained with the eight functionals (Fig. S23 in the supplementary material¹⁰²), as well as the temperature dependence of the Gibbs free energy (Fig. S24¹⁰²). Aside from the divergence in the heat capacities of PbS and PbTe predicted by PW91, PBE, and PBE-D3 at moderate temperatures, most of the functionals are able to produce the experimentally-measured values of C_p with good accuracy. The Gibbs energy curves likewise display only subtle differences in curvature; while this may be important if free energies for different phases calculated with different functionals were compared, it is difficult to see when such a situation would arise.

We have shown that the differences in athermal properties predicted by the different exchange-correlation functionals are in effect magnified at finite temperature within the quasi-harmonic approximation. We now consider the ability of the functionals to predict quantitatively the room-temperature (300 K) structure and properties of the four chalcogenides; we note that, for all four, this temperature is well below the validity limit of the QHA (∼2/3 T_m).^59,60

The 300 K lattice constants, linear-expansion coefficients, bulk moduli, and constant-pressure heat capacities obtained from the QHA calculations are collected together with experimental data in Table I. The lattice constants of the lead chalcogenides are best reproduced by PBEsol, whereas for the zinc chalcogenides, revTPSS and the two dispersion-corrected functionals also perform well. The reproduction of the linear-expansion coefficients is more variable, but the experimental values are consistently well reproduced by the LDA. This functional also yields a reasonable reproduction of the finite-temperature bulk moduli, while PBEsol performs consistently well across all four materials. The dispersion-corrected PBE-D2 and D3 functionals give good results for the zinc chalcogenides, while revTPSS, which gives a reasonable reproduction of the 300 K lattice constants of ZnS/Te, performs less well. For the zinc chalcogenides, most functionals predict constant-pressure heat capacities in good agreement with experiment. We were unable to obtain room-temperature measurements for the lead chalcogenides, but comparing the functionals shows a similar consistency, except in the case of PBE, PW91, and PBE-D3, where the larger thermal expansion leads to higher values.

Since phonon frequencies can be strongly volume dependent,³⁰ it is a common practice to perform phonon calculations with the unit-cell volume fixed at the experimentally-determined value (e.g., Refs. 72 and 91–95), or to explicitly investigate the volume dependence (e.g., Refs. 74 and 96). For a given functional, the QHA calculations yield the 300 K lattice volume on the DFT free energy surface, which would, in principle, allow finite-temperature phonon calculations to be performed without a priori knowledge of the corresponding lattice parameters. It is therefore an interesting exercise to compare these two approaches, to see which gives the better agreement with experimental data.

Fig. 7 compares the PbS and PbTe phonon dispersions between the symmetry points X, Γ, and L, calculated with the different functionals at the experimental and QHA 300 K lattice constants, to experimental measurements.

As might be expected given the variation in the predicted 300 K lattice constants in Table I, the calculations performed at the experimental volume yield more consistent phonon frequencies. This is in line with the findings of van de Walle and Ceder (Ref. 72). The largest variation in the fixed-volume calculations is seen in the frequency of the LO modes at the Γ point, which, as noted in Section IV, is mainly due to the fact that the size of the LO/TO splitting depends on the underlying electronic structure predicted by the different functionals. In keeping with this, the variation between functionals is far more pronounced for PbS than for PbTe, which is consistent with the small gap in the former causing small perturbations to lead to disproportionately large electronic responses with some functionals (i.e., the same phenomenon that led to the appearance of an imaginary mode at the X point in the equilibrium-structure LDA phonon dispersion).

Aside from the PbS LO modes, all the functionals tested appear to perform qualitatively well at predicting the phonon frequencies. For a more quantitative comparison, we calculated the average and standard deviation of the differences between the various sets of calculated frequencies and the experimental measurements (Table S9 in the supplementary material¹⁰²), excluding the Γ-point LO frequencies to avoid the comparison being skewed by the large differences in the size of the LO/TO splitting. From this analysis, we found that, at the fixed 300 K volumes, PBEsol gave the most quantitative results for both materials, while the LDA performed well for PbTe but not for PbS. Given the near-quantitative prediction of the 300 K lattice constants obtained with PBEsol, the frequencies obtained with this functional at its QHA volume were similarly good, while in these calculations, the LDA, revTPSS, and PBE-D2 functionals performed very well for PbS, and both meta-GGAs gave good results for PbTe. PBE-D3 gave better results than PBE-D2 for PbS at the experimental lattice constant, but the large thermal expansion predicted by the former leads to a higher deviation at the predicted QHA volume.

The agreement in the phonon frequencies predicted in the fixed-volume calculations implies that the phonon frequencies may be more strongly volume dependent than they are functional dependent. This is evident in a plot of the calculated Γ-point TO frequencies of the four materials against volume (Fig. 8, a corresponding plot showing the volume dependence of the Γ-point LO-mode frequencies is presented in Fig. S25¹⁰²). The calculated frequencies for the zinc chalcogenides are generally more consistent between the functionals than for the lead systems, which may account for the smaller variation in the predicted temperature dependence of the lattice constant when shifted to the athermal equilibrium volume (see Fig. S22 in the supplementary material¹⁰²). Hence, while the different functionals predict different equilibrium volumes, for a given volume/structure, they appear to predict a similar restoring force in response to displacement of the ions from their equilibrium positions, at least in the high-symmetry systems modelled here. This being the case, and as evident in Fig. 7, the practice of calculating the phonon frequencies at (fixed) experimental volumes, where available, will generally give better agreement with experiment for most DFT functionals, provided the internal strain is minimized.

Finally, another noteworthy observation is that at the QHA-predicted 300 K volumes, both PW91 and PBE predict imaginary optic modes at the Γ point, which would correspond to a permanent displacement of the ions from their ideal positions in the rocksalt structure. In the absence of LO/TO splitting, the imaginary modes are triply degenerate, while including this correction causes one mode to become real, leaving two degenerate modes with the same (imaginary) frequency as that calculated without the splitting.

To confirm the presence of these instabilities, we mapped out the total energy of the primitive cell as a function of displacement along the two imaginary-mode coordinates. Within the harmonic approximation, the displacement of the jth atom in the lth unit cell at time t is given by the following:

where N is the number of unit cells, Q(q, v) is the normal-mode coordinate, which absorbs the time dependence, and e(j, q, v) is the component of the phonon eigenvector on atom j. Considering only Γ-point modes (q = (0, 0, 0)) within a single primitive cell, Eq. (6) simplifies to

The potential surface obtained by mapping the imaginary modes in the PbS/PBE calculation is shown in Fig. 9; a similar plot for PbTe with PBE is shown in Fig. S26 in the supplementary material.¹⁰² 

The eigenvectors of the imaginary modes correspond to a mixture of offsite displacements along the three Cartesian directions, which is superficially similar to the polar rhombohedral distortion seen in the binary chalcogenides GeTe and SnTe; however, in this case, following both soft modes simultaneously would not lead to a “clean” displacement of an ion at (0.5, 0.5, 0.5) in the undistorted structure to (0.5 + x, 0.5 + x, 0.5 + x), as happens in the germanium and tin tellurides. The maps indicate a symmetric double-well potential, with the rocksalt structure lying <1 meV above the distorted structure in all cases, which is substantially less than the value of k_BT at 300 K of approximately 25 meV.

It is worth noting that as there is no macroscopic polarization exactly at the zone center, and the magnitude of the LO/TO splitting depends on the direction of approach, it would technically be more correct to map out the three Γ-point modes obtained in the absence of the splitting. Furthermore, with degenerate phonons, as here, the choice of the eigenvectors defining the basis set used to map the potential surface is somewhat arbitrary, since any linear combination of the eigenvectors would also be valid, and the maps in Figs. 9 and S26 in the supplementary material¹⁰² could thus be rotated into different bases. The present analysis is, however, sufficient to verify the presence of structural instabilities implied by the dispersions in Fig. 7.

The possible existence of instabilities is of particular interest because it would evidence the widely disputed low-temperature phase transition initially postulated by Bozin et al.^47,50–52 In the present case, however, when set against the results obtained from the other functionals, it appears to be an artefact of the anomalously-large thermal expansion predicted by PW91 and PBE. In support of this, the imaginary modes are not observed with any of the functionals at the experimental 300 K lattice volume.

We note that this result also has potential implications for constant-pressure MD simulations. If such modelling were to be performed using functionals that predict similarly anomalous thermal expansion, this could manifest in unrealistic average structures, or unexpected results such as the onset of melting at a low temperature. On the other hand, since performing variable-cell MD calculations with plane-wave codes frequently involves additional complications related to the changes in basis-set quality as the cell shape and/or volume change, MD calculations are often performed with the volumes constrained to values obtained from experimental measurements (e.g., in the theoretical literature on the dynamics of chalcogenide glasses^97–100), which the present findings suggest may largely mitigate this issue.

In summary, we have performed a quantitative comparison of the finite-temperature material properties of PbS, PbTe, ZnS, and ZnTe obtained from quasi-harmonic lattice-dynamics calculations with eight exchange-correlation functionals to experimental data.

The overarching trend is that the GGAs parameterised to favour energetics (PW91 and PBE) yield too large equilibrium volumes and too soft phonon frequencies, issues which are magnified in the QHA calculations. PBEsol, for which the GGA parameterisation is biased in favour of structural properties, yields a much better correction to the LDA and is the best general-purpose functional among those tested for these calculations. TPSS partially compensates for the PBE under binding, although revTPSS, which is optimised for bulk materials, generally performs better. However, neither meta-GGA functional displays better general-purpose performance than PBEsol in these calculations, and both are significantly more computationally expensive than GGAs. It is debatable whether dispersive interactions are important in the four chalcogenides studied here, but we do find that the D2 and D3 corrections generally offer an improvement over bare PBE. However, the improvement is inconsistent, and so these functionals may not be as generally applicable to bulk materials as standard GGAs or meta-GGAs.

In our QHA calculations, we found PBEsol to be the most capable at predicting finite-temperature lattice constants, with revTPSS and PBE-D2/D3 also performing well for the zinc chalcogenides, whereas the LDA best reproduced the experimentally-measured expansion coefficients. Accurate heat capacities could be obtained using all the functionals tested, except in the calculations on PbS and PbTe with PW91, PBE, and PBE-D3, where the anomalous thermal expansion predicted by these functionals led to a divergence of the heat capacity at moderate temperature. The variation between functionals is due in part to differences in the predicted athermal equilibrium geometries, although the relative importance of this and the predicted thermal expansion is system dependent.

All eight functionals predicted reasonable 300 K phonon frequencies, with the best consistency being observed when the calculations were performed at the cell volume given by the experimentally-measured lattice constants. Again, however, PBEsol showed the best general performance, both in the fixed-volume calculations and in those performed at the QHA volumes. The imaginary modes obtained for PbS and PbTe with PBE and PW91 at the QHA volumes were found to correspond largely to offsite displacements of the Pb cation; although this is interesting in light of the disputed low-temperature phase transition in these materials, when set against the results of the other calculations it appears to be an artefact of the anomalous thermal expansion predicted by these two functionals.

The results presented here reinforce previous studies showing that quasi-harmonic lattice-dynamics calculations represent a powerful and economical means of modelling the finite-temperature properties of bulk solids. This possibility may be of considerable value to general materials modelling and contemporary materials-discovery projects alike and could also act as a means of benchmarking the performance of electronic-structure methods. In light of its balance between computational cost and performance, we recommend PBEsol as a good starting point for lattice-dynamics calculations on bulk materials. We note, however, that the other functionals tested, e.g., the two meta-GGAs, the dispersion-corrected GGAs, or indeed other techniques such as PBE(sol) + U, may of course be more suitable for specific cases, for example, in systems with substantial dispersive interactions or strongly-correlated systems such as transition-metal oxides.

Aside from our main findings, the present study has also yielded a handful of other noteworthy observations. A reduction of the lattice constant in narrow-gap semiconductors can cause the system to undergo a metal-to-insulator transition, which for PbS appears to lead to the appearance of an imaginary phonon mode at the X point. For the LDA, which predicts the smallest lattice constants among the seven functionals tested, this occurs at the predicted equilibrium volume. For electronic-structure codes which evaluate the exchange-correlation potential from an electron-density grid, accurate calculations with meta-GGAs may require a finer grid than lower-level functionals, which would further increase the computational cost associated with using them. Finally, the PBE-D2 functional gives spurious phonon frequencies for PbTe, at least with the parameter combination used in these calculations, which we ascribe to the fact that a dispersion correction based on atomic polarizabilities may be insufficient to describe covalent bonding between highly-polarisable ions. The newer PBE-D3 correction rectifies this issue, but the lead chalcogenides nonetheless appear to represent challenging cases for both PBE-D2 and PBE-D3.

A final remark is that, in the interests of future benchmarking studies, it would be highly desirable to build a systematic database of the properties of a representative set of solids, and their temperature dependence, obtained from high-quality experimental measurements and/or state-of-the-art theoretical techniques.¹⁰¹ This would provide a reference against which new and existing functionals, and also technical methods such as the quasi-harmonic approximation, could be quantitatively evaluated.

J.M.S. gratefully acknowledges funding support from a UK Engineering and Physical Sciences Research Council Programme Grant (Grant No. EP/K004956/1). This work was carried out using the ARCHER supercomputer through membership of the UK HPC Materials Chemistry Consortium, which is funded by EPSRC Grant No. EP/L000202. We also made use of the Blue Joule facility provided by the SFTC Hartree Centre. The Hartree Centre is a research collaboratory in association with IBM, providing High Performance Computing platforms funded by the UK’s investment in e-Infrastructure. The centre aims to develop and demonstrate next generation software, optimised to take advantage of the move towards exa-scale computing. Finally, we are also grateful to Bath University Computing Services, who maintain the Aquila and Balena HPC systems.

The raw data discussed in this work is available in a repository accessible at https://github.com/WMD-Bath/Phonons. This data includes all of the structures used in the calculations, the force sets used as input to the phonon calculations, the data obtained from post processing, including the calculated phonon dispersions, densities of states, and thermal properties for each structure, and the temperature-dependent properties from the quasi-harmonic calculations.