Ab initio up to the melting point: Anharmonicity and vacancies in aluminum

We propose a fully ab initio based integrated approach to determine the volume and temperature dependent free-energy surface of nonmagnetic crystalline solids up to the melting point. The approach is based on density-functional theory calculations with a controlled numerical accuracy of better than 1 meV/atom. It accounts for all relevant excitation mechanisms entering the free energy including electronic, quasiharmonic, anharmonic, and structural excitations such as vacancies. To achieve the desired accuracy of $<1\text{meV}/\text{atom}$ for the anharmonic free-energy contribution without losing the ability to perform these calculations on standard present-day computer platforms, we develop a numerically highly efficient technique: we propose a hierarchical scheme—called upsampled thermodynamic integration using Langevin dynamics—which allows for a significant reduction in the number of computationally expensive ab initio configurations compared to a standard molecular dynamics scheme. As for the vacancy contribution, concentration-dependent pressure effects had to be included to achieve the desired accuracy. Applying the integrated approach gives us direct access to the free-energy surface $F\left(V,T\right)$ for aluminum and derived quantities such as the thermal expansion coefficient or the isobaric heat capacity and allows a direct comparison with experiment. A detailed analysis enables us to tackle the long-standing debate over which excitation mechanism (anharmonicity vs vacancies) is dominant close to the melting point.

Figure 1

(Color online) $\lambda$ dependence of the ensemble average quantity ${⟨\partial U/\partial \lambda ⟩}_{\lambda }$ for the perfect (perf) bulk and vacancy (vac) calculations for a small and a large volume each. The results were obtained for a $2^3$ supercell, the LDA functional, and a temperature of 900 K. The standard deviations (Ref. 44) $\left(\approx 0.5\text{meV}/\text{atom}\right)$ are of a similar size as the symbols indicating the calculated values. The solid and dashed lines are guides for the eyes.

Figure 2

(Color online) Explicitly anharmonic free energy $F^{\text{ah}}$ at 900 K for the LDA functional. Results for the $2^3$ and $3^3$ supercell (sc) and for three different $\lambda$ samplings, $\lambda =0.5$, $\lambda$ $\text{set}1=\left\{0.1,0.3,0.5,0.7,0.9\right\}$, and $\lambda$ $\text{set}2=\text{set}1+\left\{0.02,0.98\right\}$, used for the integration in Eq. (5) are shown. The symbols represent the calculated values, the corresponding vertical lines show the standard deviation (Ref. 44), and the dashed and dotted lines in between are guides for the eyes. The vertical black dashed lines indicate the 0 K (including zero-point vibrations) and the melting temperature equilibrium volumes of LDA, $V_{\text{T}0}$ and $V_{\text{Tm}}$, respectively.

Figure 3

(Color online) Illustration of the UP-TILD method. (a) Convergence of the average quantity ${⟨\partial U/\partial \lambda ⟩}_{\lambda }^{\text{low}}$ (at $V=15.7{\text{Å}}^3$, $T=900\text{K}$, and $\lambda =0.5$) as a function of the number of MD steps for a “usual” thermodynamic integration run as discussed in Sec. 2C1. The superscript indicates that the run is performed with low-converged parameters. (b) Convergence of the average energy difference ${⟨\Delta E⟩}_{\lambda }$ for a set of uncorrelated MD structures taken from the run in (a). The difference is obtained from the low and high converged energies according to Eq. (8). (c) The $\lambda$ dependence of the converged quantities from (a) and (b) and their difference ${⟨\partial U/\partial \lambda ⟩}_{\lambda }^{\text{low}}-{⟨\Delta E⟩}_{\lambda }={⟨\partial U/\partial \lambda ⟩}_{\lambda }^{\text{high}}$, which corresponds to the upsampled quantity ${⟨\partial U/\partial \lambda ⟩}_{\lambda }^{\text{high}}$, are shown.

Figure 4

(Color online) Explicitly anharmonic free energy $F^{\text{ah}}$ at 900 K for the LDA functional and the $2^3$ supercell. Results for low and high converged parameters are shown: low: 8 Ry, 256 kp⋅atom and high: 14 Ry, 2048 kp⋅atom. Free energies obtained directly from a usual thermodynamic integration calculation are marked as “direct,” while “UP-TILD” indicates that the UP-TILD method was employed.

Figure 5

Schematic illustration of the concept of calculating the free energy $F\left(\Omega ,T;N,n,{\Omega }^v,N^v\right)$ of a periodic crystal with vacancies. The larger box represents a supercell of volume $\Omega$ at temperature $T$ containing $N$ atoms and $n$ vacancies. A light-gray shaded box with a white circle represents a cell of volume ${\Omega }^v$, containing $N^v$ atoms, exactly one vacancy, and having the free energy $F^v\left({\Omega }^v,T;N^v\right)$. The dark-gray shaded region represents the perfect bulk without vacancies, with volume $\Omega -n{\Omega }^v$, $N-n{N}^v$ atoms, and free energy $F^b\left(\Omega -n{\Omega }^v,T;N-n{N}^v\right)$. The dashed lines indicate periodic boundary conditions.

Figure 6

(Color online) The contribution of the electronic (el), explicitly anharmonic (ah), and vacancy (vac) excitations to the free energy $F_P$ at constant pressure, the expansion coefficient $\alpha$, and the heat capacity $C_P$ at constant pressure (all quantities at zero pressure). The black solid lines show the sum of all excitations. Thick lines show GGA results. LDA results are either shown as thin lines or coincide with the GGA result. The right axes are scaled with respect to the “full” GGA values at the melting temperature $T^m$ (indicated by the vertical dashed line) given in Table . The inset shows an enlargement of the vacancy contribution to $F_P$ at high temperatures.

Figure 7

(Color online) Quasiharmonic expansion coefficient and heat capacity at zero pressure. The solid lines correspond to fully converged quantities. The dotted-dashed/dotted/dashed lines represent quantities based on approximation (a)/(b)/(c) [see Eqs. (26, 27, 28) for details]. The inset shows the corresponding free energies $\Delta F^{\text{qh}}$ (per atom) at 900 K referenced to the fully converged free energy. The curves correspond to GGA results. LDA results show the same behavior. The vertical thin dashed line marks the melting temperature $T^m$.

Figure 8

(Color online) Explicitly anharmonic free energy $F^{\text{ah}}$ at 900 K for the two investigated XC functionals. The dots/triangles represent the calculated values, while the corresponding vertical lines indicate the standard deviation (Ref. 44). The solid lines are fits using the analytical model (see text). The vertical dashed (dotted-dashed) lines indicate the 0 K (including zero-point vibrations) and melting temperature equilibrium volumes of LDA (GGA), $V_{\text{LDA}0}$ $\left(V_{\text{GGA}0}\right)$ and $V_{\text{LDAm}}$ $\left(V_{\text{GGAm}}\right)$, respectively.

Figure 9

(Color online) Explicitly anharmonic free energy $F^{\text{ah}}$ at three different volumes for the LDA functional. The calculated values are represented by the dots/diamonds/triangles. At each point the statistical error is represented by the vertical solid lines. The dashed lines in between are a guide for the eyes and the solid lines are fits using the analytical model (see text). GGA results show the same qualitative dependence. The melting temperature $T^m$ is indicated by the vertical dashed line.

Figure 10

(Color online) Illustration of the analytical model (see text) describing the origin of the anharmonic free energy. (a) Volume independent anharmonic frequency ${\overline{\omega }}^{\text{ah}}$. (b) Effective frequency resulting from the sum of the quasiharmonic frequency ${\overline{\omega }}^{\text{qh}}$ and ${\overline{\omega }}^{\text{ah}}$. (c) Model anharmonic free energy $F^{\text{ah}}$ at the melting temperature as obtained from Eq. (29) and the frequencies from (b).

Figure 11

(Color online) Equilibrium vacancy concentration as a function of the inverse temperature multiplied by the melting temperature $T^m$. Results for the two investigated exchange-correlation functionals, LDA, and GGA, are shown. The solid lines correspond to a volume-optimized (volOpt) calculation, i.e., one which is based on solving the coupled system of Eqs. (16, 17, 18, 19), with all free-energy terms included in Eq. (24). The dashed lines correspond to a volume-optimized calculation excluding the anharmonic free-energy terms, $F^{b,\text{ah}}$ and $F^{v,\text{ah}}$, in Eq. (24). The dashed lines also correspond to a volume-optimized calculation excluding the anharmonic, $F^{b,\text{ah}}$ and $F^{v,\text{ah}}$, and electronic, ${\overline{F}}^{b,\text{el}}$ and ${\overline{F}}^{v,\text{el}}$ (indicated by the parentheses enclosing the electronic contribution in the legend) terms in Eq. (24). The dotted lines correspond to a calculation based on the approximation (approx) Eq. (21) and including all free-energy terms. The squares indicate experimental values from Ref. 62 (differential dilatometry). The diamonds/circles indicate experimental values from Ref. 63 (differential dilatometry/positron annihilation).

Figure 12

(Color online) Explicitly anharmonic free energy $F^{\text{ah}}$ of the vacancy calculation at three different volumes for the LDA functional. The dots/diamonds/triangles represent the calculated values. At each point the statistical error is represented by the vertical solid lines. The dashed lines in between are guides for the eyes and the solid lines are fits using the analytical model (see Sec 3B). GGA results show the same qualitative dependence. The melting temperature $T^m$ is given by the vertical dashed line.

Figure 13

(Color online) Free energy $F_P$ at constant pressure for the two investigated exchange-correlation functionals: LDA and GGA, in comparison with CALPHAD results (obtained using THERMOCALC version $Q$ and the “PURE4-SGTE pure elements database”). Both the ab initio and the CALPHAD data were shifted to zero at $T=200\text{K}$. For the difference free energy $\Delta F_P$ the CALPHAD data are taken as a reference at each temperature. The melting temperature $T^m$ is given by the vertical dashed line.

Figure 14

(Color online) Heat capacity at constant pressure for the two investigated exchange-correlation functionals, LDA and GGA, in comparison with experiment. The melting temperature $T^m$ is given by the vertical dashed line. Experimental data older than 1950 are indicated by the open circles (Refs. 7, 8, 9). The remaining experimental data are indicated by the filled squares (Refs. 2, 5, 6, 10, 11, 12, 13, 14, 15). The inset shows the low-temperature region with experimental data from Ref. 6.

Figure 15

(Color online) Thermal expansion coefficient for the two investigated exchange-correlation functionals, LDA and GGA, in comparison with experiment. The melting temperature $T^m$ is given by the vertical dashed line. Experimental data older than 1950 are indicated by the open circles (Refs. 67, 68, 69, 70). The remaining experimental data are indicated by the filled squares (Refs. 71, 72, 73, 74, 75, 76). The inset shows the low-temperature region with experimental data from Ref. 77.