Calculation of the Anisotropic Coefficients of Thermal Expansion: A First-Principles Approach

TL;DR

This paper introduces a first-principles method that incorporates anharmonic effects to accurately predict anisotropic thermal expansion coefficients and other thermal properties of materials across wide temperature ranges.

Contribution

The authors develop a novel computational approach that extends beyond quasiharmonic approximation to include anharmonic effects for thermal property prediction.

Findings

Accurately predicts thermal expansion coefficients for Cu and AlN.

Shows good agreement with experimental data over wide temperature ranges.

Provides a comprehensive method for thermal property calculations in anisotropic materials.

Abstract

Predictions of the anisotropic coefficients of thermal expansion are needed to not only compare to experimental measurement, but also as input for macroscopic modeling of devices which operate over a large temperature range. While most current methods are limited to isotropic systems within the quasiharmonic approximation, our method uses first-principles calculations and includes anharmonic effects to determine the temperature-dependent properties of materials. These include the lattice parameters, anisotropic coefficients of thermal expansion, isothermal bulk modulus, and specific heat at constant pressure. Our method has been tested on two compounds (Cu and AlN) and predicts thermal properties which compare favorably to experimental measurement over a wide temperature range.

Tables

Table 1. Table 1: Our ground state (0 K) calculated Debye temperature, isentropic bulk modulus, and lattice parameters compared to experimental values. The experimental lattice parameters of Cu and AlN come from experimental data extrapolated to zero temperature.

	Cu			AlN
	Calc.	Exp.	Ref.	Calc.	Exp.	Ref.
${\Theta }_D$ (K)	$348$	$347$	65	$1048$	$971$	66
$B$ (GPa)	$143.8$	$137.6$	67	$196.5$	$201$	68
$a$ (Å)	$3.633$	$3.615$	69	$3.129$	$3.107$	70
$c$ (Å)				$5.012$	$4.973$	70

Equations

$$B$$

$$G$$

$$\displaystyle+3(C_{44}+C_{55}+C_{66})\bigg{)}.$$

$$v_l=\sqrt{\frac{B+4/3G}{\rho }}$$

$$v_s=\sqrt{\frac{G}{\rho }},$$

$${\Theta }_D=\frac{h}{k_B}{\left[\frac{9m}{4\pi }\right]}^{1/3}{\left[\left(\frac{1}{v_l^3}+\frac{2}{v_s^3}\right)\right]}^{-1/3},$$

$$F(\mathbf{x},T)$$

$$F_{\text{vib}}(\mathbf{x},T)$$

$$\displaystyle+k_{B}Tg(\omega_{\lambda})\ln\left(1-\exp\left(\frac{\hbar\omega_{\lambda}}{k_{B}T}\right)\right)d\omega_{\lambda}\Bigg{]},$$

$$F_{\text{fit}}(\mathbf{x})=\sum _{i,j}{f}_{i,j,k}{a}^i{b}^j{c}^k.$$

$${\alpha }_{x_i}=\frac{1}{x_i^0}\frac{d{x}_i}{dT}=\frac{d\ln (x_i)}{dT}$$

$$C_p=C_v+VT\frac{\mathrm{Tr}(\alpha {)}^2}{{\beta }_T},$$

Full text

I Introduction

While experimental measurements of the coefficient of thermal expansion (CTE) can be done using a number of experimental techniques Belousov and Filatov (2007); Cliffe and Goodwin (2012); Jones et al. (2013); Samovilov et al. (2006); Dudescu et al. (2006), most theoretical tools available for predicting the CTE are limited to the quasiharmonic approximation (QHA) or the Debye-Gruneisen approximation. There are a few programs Togo and Tanaka (2015); Togo et al. (2010) and several models Jin and Wu (2002); Quong and Liu (1997); Tohei et al. (2016); Liu, Zhou, and Wang (2017); Wang (2006); Wang and Reeber (2011); Ding and Xiao (2015); Lamowski et al. (2011); Liu and Zhou (2018); Gilev (2018) available for scientists to calculate the CTE, and its associated thermal properties, from first-principles calculations. These currently available methods principally rely on the QHA. They allow the user to determine the temperature dependence of the isotropic lattice parameters, the coefficient of volume expansion, isothermal bulk modulus and specific heat at constant pressure. Since these programs determine the coefficient of volume expansion, they are limited to calculations of systems with a considerable amount of structural symmetry. A limited amount of work has been published on the thermal properties of anisotropic systems employing anharmonic effects Shulumba et al. (2016, 2015); Schmerler and Kortus (2014), but this work is not available as a comprehensive package of tools that can automatically predict such properties from the crystal structure of a material.

In this work, an algorithm is developed to determine the thermal properties of both isotropic and anisotropic systems using density functional theory (DFT) and density functional perturbation theory (DFPT) calculations. This algorithm utilizes the Vienna Ab initio Simulation Package (vasp) Kresse and Hafner (1993); Kresse and Furthmuller (1996a, b) for DFT and DFPT calculations and the Temperature Dependent Effective Potential (tdep) Hellman and Abriksov (2013); Hellman et al. (2013); Hellman, Abriksov, and Simak (2011) software package to calculate the ground state and interatomic force constants using a method designed to work for both harmonic and strongly anharmonic systems Hellman, Abriksov, and Simak (2011).

To demonstrate the abilities of our methods, two systems with distinctive symmetries are compared to a variety of experimental results. The first system, pure copper (Cu) Hahn (1970); White (1973), possesses cubic symmetry and therefore the CTE is a single temperature-dependent value. The temperature dependent properties of copper are thoroughly investigated experimentally Touloukian et al. (1975); Hahn (1970); Adenstedt (1936); Carr, McCammon, and White (1964); Carr and Swenson (1964); Esser and Eusterbrock (1941); Fraser and Hallett (1965); Lifanov and Sherstyukov (1968); Rubin, Altman, and Johnston (1954); Arblaster (2015) and Cu is consequently used to determine the quality of our calculation methods. The second material, aluminum nitride (AlN) Figge et al. (2009), has a hexagonal crystal structure and therefore two unique non-zero CTE’s and is used to demonstrate our ability to calculate the expansion coefficients of an anisotropic material with a known temperature dependent anharmonic interaction Shulumba et al. (2016).

In section II, an outline of the calculation methods is provided and in section III, we analyze the results of our methods via comparison to available experimental data for Cu and AlN. Finally, in section IV, we provide a concluding discussion. An appendix is provided that outlines the algorithm used to generate the necessary files for the vasp and tdep calculations, gather the necessary results, and perform the needed temperature dependent calculations.

II Calculation Methods

Our software launches both DFT and DFPT calculations using vasp Kresse and Furthmuller (1996a, b). The present benchmark calculations were performed with the projector-augmented wave (PAW) method Blochl (1994) to describe the core electrons by utilizing PAW pseudopotentials Kresse and Joubert (1999). Exchange and correlations effects were described by the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) Perdew, Burke, and Ernzerhof (1996). To describe the electronic system, a plane-wave energy cutoff equal to 500 eV was used and a Monkhorst-Pack Monkhorst and Pack (1976) mesh of points was generated for each grid in reciprocal space assuming a $k$-point density of at least five points per Å*-1*. For the relaxation, ground-state, and configuration calculations the same $k$-point mesh was used. DFPT calculations were done on a $k$-point mesh twice as dense as the grid used in the ground state calculations.

For each ground state and configuration calculation, the iterations of the total energy were stopped once the differences in energy between successive iterations were less than 0.01 meV. To calculate the dielectric, elastic, and Born effective charge of our material, we used finite-differences as implemented within vasp where only symmetry-inequivalent atomic displacements were used to calculate the Hessian matrix. The elastic tensor was used to determine the Debye temperature and the isentropic bulk modulus, as outlined below. The dielectric tensor and Born effective charge tensor were used to account for the long-range electrostatic interactions in semiconducting compounds within tdep.

It is important to note that there is a fundamental relationship between the isoentropic bulk modulus (via the elastic tensor) and the isothermal bulk modulus (found by fitting the equation of state) directly related to the ratio of specific heats. While this ratio is nearly identical to one (see below) we will distinguish here due to the differences in the methodology used to determine the quantities.

Using the elastic tensor, one can determine both the Debye temperature and isentropic bulk modulus of the system using the following relationships. First, the Voigt average of the bulk ($B$) and shear ($G$) moduli (upper bound for polycrystalline materials) are found from the components of the elastic tensor, $C_{ij}$ ($i,j=1..6$), as Nye (1957); Pike et al. (2018)

[TABLE]

With the bulk and shear moduli one can then calculate the longitudinal and shear sound velocities ($v_{l}$ and $v_{s}$) as Jia, Chen, and Zhang (2017)

[TABLE]

and

[TABLE]

where $\rho$ is the density of the material. Then, the Debye temperature can be calculated as Jia, Chen, and Zhang (2017)

[TABLE]

where $h$ is Plank’s constant, $k_{B}$ is Boltzmann’s constant, and $m$ is the number of atoms per volume and, unlike in Ref. 46, we do not separate the acoustic bands and use all the phonon bands to calculate the Debye temperature.

Since the Debye temperature is calculated from the elastic tensor, which is calculated at zero temperature within vasp, it also represents the ground state (0 K) Debye temperature. While one could, in principle, calculate the elastic properties as a function of temperature Steneteg et al. (2013), and therefore the Debye temperature as a function of temperature, we have found that the Debye temperature at zero temperature provides an adequate starting point for our calculations since variations in the configuration temperature, outlined below, do not significantly modify our results.

The Debye temperature is used within tdep to generate an initial guess for the force constants of the system. They are in turn employed to generate configurations based on a canonical ensemble. Here, twelve configurations are generated as follows: with the Debye temperature and symmetry of the cell, one generates an initial guess for the interatomic force constants and solves the equations of motion for the system by finding the resulting eigenvalues and eigenvectors. Then, these eigenvalues and eigenvectors are used to determine the amplitude and velocity of each atom chosen such that these quantities are normally distributed over a canonical ensemble. The configurations consisted in this study of 208 atoms for Cu and AlN. The ensemble is generated at a finite temperature (here designated the configuration temperature) with Bose-Einstein statistics used to determine the mean normal mode amplitude  West and Estreicher (2006). Vasp is then used to calculate the energies, displacements, and forces for each of the generated configurations. These are combined within tdep to generate the finite-temperature force constants by fitting the Born-Oppenheimer energy surface. This procedure allows us to go beyond the quasiharmonic approximation by explicitly including anharmonic effects through the canonical ensembles.

The extracted force constants are used to find the phonon frequencies and phonon density of states on a $q$ mesh in reciprocal space (in this work $30\times 30\times 30$ points were used) using the tetrahedron integration approach Lehmann and Taut (1972). The density of states is used to calculate thermal properties, as outlined below. Using this grid and integration method ensured the convergence of the free energy to within 0.01 meV/atom.

The calculated lattice parameters ${\bf x}(T,p)$ depend on temperature and pressure. Determining the lattice parameters at ambient pressure ${\bf x}(T,0)$ requires minimizing the free energy $F({\bf x},T)$ at each temperature. $F$ can be expressed as a function of the total electronic energy $U(\bf{x})$ and the vibrational free energy ($F_{\textrm{vib}}({\bf x},T)$) as

[TABLE]

where $\lambda$ is the phonon mode index and both the vibrational density of states ($g(\omega)$) and phonon frequencies $\omega_{\lambda}$ depend on $\bf{x}$. $U$ is the volume dependent total electronic energy of the system. In order to minimize the free energy, a set of lattice parameters is generated automatically by applying strain to the system around the ground state equilibrium lattice parameters. This was in the present work done by finding a series of six lattice parameters in each symmetry unique direction such that the new lattice parameters range between one percent compressive and four percent tensile strain. This provided a smooth free energy surface at each temperature.

The optimized lattice parameters ${\bf x}_{0}(T)$ are found by fitting the free energy to a polynomial function Tohei et al. (2016); Wang (2006); Schmerler and Kortus (2014) of the lattice parameters $x_{i}=a,b$, and $c$:

[TABLE]

$f_{i,j,k}$ are the coefficients of the polynomial fit. Here, we used fourth-order polynomials, and varied $a$ ($a$ and $c$) in the case of Cu (AlN).

To determine the minimum of Eq. (6), we use the constrained BFGS minimization method as implemented in the Scipy optimize package Jones et al. (1); Byrd et al. (1995); Zhu et al. (1997) within Python. After storing these minimizing parameters they are fit to an eighth-order polynomial as a function of temperature to account for the numerical noise in the free energy calculations. The CTE is then calculated by computing the derivative of the smoothed data Slack and Bartram (1975) as

[TABLE]

for the different components of the CTE. Where $x^{0}_{i}$ is either the zero temperature lattice constant or the temperature-dependent lattice constant, as shown here. Using either of these produces practically identical results as assumed by Slack Slack and Bartram (1975).

With the integrated phonon density of states, the specific heat at constant volume can be calculated with tdep for the relaxed geometry. To determine the specific heat at constant pressure ($C_{p}$) the well-known thermodynamic relationship between the specific heat at constant volume, ($C_{v}$) and the trace of the thermal expansion coefficient ($\alpha$) calculated with Eq. (7), can be used:

[TABLE]

where $\beta_{T}$ is the isothermal compressibility which is the inverse of the isothermal bulk modulus $B$. To determine $B$, we fit the free energy at a fixed temperature versus volume of the cell using the Birch-Murnaghan equation of state Birch (1947, 1978); Murnaghan (1944). By fitting this equation of state, we can extract the cohesive energy, isothermal bulk modulus, and the pressure derivative of the bulk modulus as a function of temperature.

III Results

To test the accuracy of our vasp and tdep calculations and the accuracy of our algorithm, we will compare our calculated values to experimental results for well-known systems: bulk Cu and AlN. Bulk Cu is a face-centered cubic structure (space group $Fm\bar{3}m$) from 0 K to its melting point of approximately 1357 K Brand et al. (2006); Preston-Thomas (1990). Bulk AlN has a wurtzite crystal structure (space group $P6_{3}mc$) from 0 K to its melting point of approximately 3270 K MacChesney, Bridenbaugh, and O’Conner (1970). For each of these materials, our calculations can be compared to experimental measurement as shown in Table 1 and in the figures within the text. Table 1 contains a comparison of our ground state (0 K) calculated Debye temperature, isentropic bulk modulus, and DFT lattice parameters to experimental measurement. The agreement between our values and experiment is well within the typical errors seen in ground state GGA-DFT calculations.

The temperature dependent lattice parameters of Cu and AlN are shown in Figs. 1 and 2, respectively. The difference between our calculated values and experimental measurement comes from the inherent error in predicting the lattice parameters at the specific level of theory (i.e.: the PBE GGA of DFT). To compare to experiment over a broad temperature range for Cu, we have used the recommended values of the thermal lattice expansion, given in Ref. 29.

We performed a polynomial fit to the extracted lattice parameters, shown as the solid lines in Figs. 1 and 2. In addition we plot lattice parameters shifted by the ground-state error (difference between DFT calculated and experimentally extrapolated values at $T=0$ K) with black dashed lines for the sake of comparison. The un-shifted polynomial fit was used to calculate the CTE, specific heat, and bulk modulus. We obtain excellent agreement between the calculated CTE and experimental measurements for Cu and AlN between approximately 40 K and the Debye temperature, as shown in Figs. 3 and 4. However, there is some deviation above the Debye temperature, most notably for Cu where the calculated CTE is markedly higher than the experimental values as the temperature approaches the melting temperature. This is most likely due to higher-order anharmonic interactions not sufficiently accounted for at these temperatures. For AlN, it has been shown that the temperature dependent $c/a$ ratio is not correctly reproduced in anharmonic DFT/tdep calculations Shulumba et al. (2016), and our extracted anisotropic CTE of AlN are consistent with those results.

The calculated specific heats of these materials also agree very well with experimental measurement as shown in Figs. 5 and 6. We also calculated $C_{p}$ for both materials, but the difference between $C_{v}$ and $C_{p}$ was too small over the entire temperature range of data to be discernible in the plot—at any temperature the difference $1-(C_{p}/C_{v})$ was smaller than $10^{-6}$. This calculation, therefore, is further evidence that one can safely ignore the differences between these two quantities when comparing experimental and theoretical values Cardona et al. (2007).

Finally, during our fitting of the equation of state, we extract the isothermal bulk modulus and its pressure derivative at zero external pressure as a function of temperature. The isothermal bulk modulus and the CTE are used to calculate the specific heat at constant pressure for both Cu and AlN. In Figs. 7 and 8, we plot the isothermal bulk modulus and pressure derivative vs temperature for both Cu and AlN, respectively. Literature values of the isoentropic bulk modulus are shown for the sake of comparison and lie above our extracted values. A technical explanation of our method is outlined in the appendix.

IV Discussion and Conclusion

Our algorithm calculates the thermal properties of both isotropic and anisotropic materials using the DFT, DFPT and TDEP methods. Several related thermal properties are also calculated including the isothermal bulk modulus and specific heat at constant pressure. Agreement between our calculations and experimental results for the temperature dependent lattice parameters, specific heat at constant pressure, and isothermal bulk modulus give us confidence in our calculation methods. The correspondence between our calculations and experiment for the CTE is only moderate at temperatures above the Debye temperature. This may e.g. be due to unattributed numerical inaccuracies, the present level of theory (PBE-GGA), or due to an inadequate sampling of the anharmonic interactions at elevated temperatures due to the finite number of configurations used in this calculation. Nevertheless, most of our predictions are very well aligned with experiment, which makes us confident that this software can be used to determine the CTE for many more materials of experimental interest either with individual calculations or as part of a high-throughput framework.

We have undertaken several careful convergence studies of the number of configurations, super-cell size, and $k$-point mesh density to minimize the error between our calculations and experimental measurement over the entire temperature range. We have determined that using the $k$-point mesh density given above provides a reasonable accuracy between our calculated CTE and experiment. Additionally, we have determined that the configuration temperature used by tdep should be slightly smaller than the Debye temperature to account for a reasonable amount of the anharmonic interactions as the temperature approaches the melting temperature. Here, we show results with the configuration temperature equal to 80% of the Debye temperature and have tested both higher and lower configuration temperatures as well with similar results.

While not shown here, the present algorithm can be used to calculate the CTE from a completely anisotropic system. Additionally, this algorithm can easily be integrated into existing high-throughput software workflows that we hope will enable subsequent researchers in e.g. designing the next generation temperature-dependent electronics. While one must use caution when calculating temperature-dependent properties, and carefully consider the convergence of many of the DFT, DFPT, and tdep input parameters, this software thus provides an algorithm for automatically calculating the CTE and related temperature-dependent properties of any material with a periodic unit cell. The symmetry and size of this unit cell may be limiting factors of the choice of materials, since it may be prohibitively expensive to perform calculations on very large unit cells with low symmetry with currently available high performance computing resources.

Acknowledgments

We would like to thank the Research Council of Norway through the Frinatek program for funding and Martin Fleissner Sunding, Olle Hellman, and Nina Shulumba for helpful conversations and support. The computations were performed on resources provided by UNINETT Sigma2 - the National Infrastructure for High-Performance Computing and Data Storage in Norway through grant numbers nn2615k and nn9462k.

Appendix A Summary of the Algorithm

Our algorithm for calculating the CTE and associated thermal properties is diagrammed in Fig. 9 which involves both serial and parallel computations using vasp and tdep. To start a calculation, the user enters a single vasp formatted POSCAR file describing the material’s crystal structure into the directory where the entire calculation will take place. This POSCAR file should represent the conventional cell of the system in the desired symmetry and, preferably, should represent the DFT-relaxed ground state structure. To launch the calculations, the user should execute the script with the tag ”–relaxation” which generates the necessary input files and bash scripts, creates the directories for the vasp calculations, and checks the convergence of the calculations once completed. Following a relaxation calculation to ensure convergence, the script will automatically launch the finite-differences calculations for the needed tensor quantities. Should the algorithm detect an error in the data or convergence procedure, the code will abort. Once the error is fixed, the user can relaunch the calculation.

After the finite-differences calculation, our algorithm automatically starts the next stage of the calculation (which can be manually executed with –build_cells) that constructs a set of unit cells corresponding to different volumes. The algorithm generates the necessary lattice parameter perturbations, as described above, to determine the CTE for the compound. During this calculation, the script will analyze the symmetry of the relaxed unit cell, calculate the Debye temperature for the creation of the initial set of configurations used by tdep to determine the inter-atomic force constants, and generate the necessary files for the related vasp calculations. Once generated, the algorithm will rerun the relaxation of the newly perturbed unit cell by allowing relaxation of only the atomic positions. After this second relaxation, the script will generate the configurations and launch each configuration in parallel. Since this is the most time-consuming part of the calculation, the user can relaunch this calculation as needed and only the parts of the calculation that are incomplete will be executed.

The finishing stage of this calculation can be executed with the tag ”–thermal_expansion” in which the script, for each lattice perturbation, will launch tdep to post-process the results of the previous vasp calculations. This post-processing generates the phonon density of states and various thermal properties including the free energy of the system. After gathering the calculated free energies, the script will determine the set of lattice parameters that minimize the free energy as a function of temperature. This temperature dependent set of lattice parameters is used to determine the CTE, isothermal bulk modulus, and specific heat at constant pressure as outlined above.

During the concluding stage of the calculation, several output files are produced. The calculated free energies for each volume as a function of temperature are printed to a single data file ”out.free_energy_vs_temp” with temperatures given in Kelvin and free energies given in eV. The temperature dependent lattice parameters are printed to a file called ”out.thermal_expansion” with the lattice parameters given in Å. The thermal expansion coefficients, given in Eq. (7), are printed to the file ”out.expansion_coeffs” in units of K*-1*. Calculations of the isothermal bulk modulus, its pressure derivative, and cohesive energy are printed to the file ”out.isothermal_bulk” where the isothermal bulk modulus is given in GPa, the derivative of the bulk modulus with respect to pressure is unit-less, and the cohesive energy is in eV. Finally, the specific heat at constant volume and constant pressure are calculated using Eq. (8). These are printed to the data file ”out.cv_cp” where both specific heats are given in units of kJ kg*-1* K*-1*.

This software, and corresponding documentation, is available by direct request to the corresponding author.