A fast anharmonic free energy method with an application to vacancies in ZrC

TL;DR

This paper introduces a rapid and accurate method for calculating anharmonic free energies in crystals, demonstrated on ZrC with predictions on vacancy behavior and thermodynamics.

Contribution

A novel computational approach that significantly accelerates anharmonic free energy calculations with high accuracy, applicable to complex materials.

Findings

Achieves 10x speed-up over existing methods

Predicts vacancy concentrations in ZrC_x

Provides accurate heat capacity estimates

Abstract

We propose an approach to calculate the anharmonic part of the volumetric-strain and temperature dependent free energy of a crystal. The method strikes an effective balance between accuracy and computational efficiency, showing a $\times10$ speed-up on comparable free energy approaches at the level of density functional theory, with average errors less than 1 meV/atom. As a demonstration we make new predictions on the thermodynamics of substoichiometric ZrC$_x$, including vacancy concentration and heat capacity.

Tables

Table 1. Table 1: Computer resources consumption for Zr 32 C 32 test case calculation of F ah ​ ( V , T ) subscript 𝐹 ah 𝑉 𝑇 F_{\text{ah}}(V,T) on a 5 × 5 5 5 5\times 5 { V i , T i } subscript 𝑉 𝑖 subscript 𝑇 𝑖 \{V_{i},\,T_{i}\} mesh. Timings listed in CPU core-hours and quasiharmonic-free-energy-job units ( 1 / t ​ ( F qh ) 1 𝑡 subscript 𝐹 qh 1/t(F_{\text{qh}}) ) using the reference value t ​ ( F qh ) = 4800 𝑡 subscript 𝐹 qh 4800 t(F_{\text{qh}})=4800 core-hours.

Contributions	$t\left(F_{\text{ah}}\right)$ (core-hrs)			$t\left(F_{\text{ah}}\right)/t\left(F_{\text{qh}}\right)$
	this work	TU-TILD		this work	TU-TILD
Fit set DFT MD	${10}^4$	${10}^3$		10	0.6
MEAM fitting	${10}^3$	${10}^2$		0.4	0.03
MEAM TI	${10}^3$	${10}^2$		0.4	0.1
DFT TI	-	${10}^5$		-	110
DFT up-sampling	-	${10}^4$		-	6
Total	${10}^5$	${10}^6$		11	117

Table 2. Table 2: Concentration of vacancies ( c vac subscript 𝑐 vac c_{\text{vac}} ) in ZrC x in carbon atom %.

$T$ (K)	$c_{\text{vac}}$ (C at. %)
	this work	CALPHADFernández Guillermet (1995)
300	$9\times {10}^{-13}$	$1\times {10}^{-15}$
500	$6\times {10}^{-7}$	$1\times {10}^{-9}$
1000	0.01	0.003
1500	0.33	0.06
2000	1.87	0.33
2200	3.03	0.51

Equations

$$F=E_0+F_{\text{el}}+F_{\text{qh}}+F_{\text{ah}}+F_{\text{el-vib}}+F_{\text{config}},$$

$$F_{\text{ah}}={\int }_0^{1}d\lambda {⟨{\partial }_{\lambda }{E}_{\text{mix}}(\lambda )⟩}_{\lambda }.$$

$$F_{\text{ah}}^{\text{TUTILD}}$$

$$={\int }_0^1\text{d}\lambda {⟨E^{\text{MEAM}}(\mathbf{R},V)-E_{\text{qh}}^{\text{DFT}}(\mathbf{R},V)⟩}_{\lambda }$$

$$+{\int }_0^1\text{d}\lambda {⟨E^{\text{DFT}}(\mathbf{R},V)-E^{\text{MEAM}}(\mathbf{R},V)⟩}_{\lambda }.$$

$$F_{\text{ah}}^{}={\int }_0^{1}d\lambda {⟨E^{\text{MEAM}}(\mathbf{R},V)-E_{\text{qh}}^{\text{DFT}}(\mathbf{R},V)⟩}_{\lambda }.$$

$${⟨\Delta E⟩}_{\lambda }(\lambda )=\sum _{i=0}^{i=5}{a}_i{\lambda }^i,\lambda \in [0,1],$$

$$\delta =\sqrt{{\delta }_{\text{sys}}^2+{\delta }_{\text{rand}}^2+{\delta }_{\text{ref}}^2}.$$

$${\delta }_{\text{rand}}\sim \sqrt{\frac{{\sigma }_{\lambda }^2}{N_{\lambda }}\frac{\tau }{t}},$$

$${\delta }_{\text{sys}}={\int }_0^{1}d\lambda {⟨E^{\text{DFT}}(\mathbf{R},V)-E^{\text{MEAM}}(\mathbf{R},V)⟩}_{\lambda }.$$

$$F_{\text{ah}}^{\text{excess}}=F_{\text{ah}}\left({\text{Zr}}_{32}{\text{C}}_{31}\right)-F_{\text{ah}}\left({\text{Zr}}_{32}{\text{C}}_{32}\right),$$

$$G_{\text{form}}^{}=G\left({\text{Zr}}_{32}{\text{C}}_{31}\right)+\mu (\text{C})-G\left({\text{Zr}}_{32}{\text{C}}_{32}\right).$$

$$\frac{c_{\text{vac}}}{1-c_{\text{vac}}}=\text{exp}\left(-\frac{G_{\text{form}}}{k_{\text{B}}T}\right),$$

$$G_{\text{ah}}=\underset{V}{\min }\left[F_{\text{ah}}+F_{\text{qh}}+E_0+pV\right]-\underset{V}{\min }\left[F_{\text{qh}}+E_0+pV\right],$$

$$G_{\text{el-vib}}=\underset{V}{\min }\left[F_{\text{el-vib}}+F_{\text{qh}}+E_0+pV\right]-\underset{V}{\min }\left[F_{\text{qh}}+E_0+pV\right],$$

$$G_{\text{el}}=\underset{V}{\min }\left[F_{\text{el}}+F_{\text{qh}}+E_0+pV\right]-\underset{V}{\min }\left[F_{\text{qh}}+E_0+pV\right].$$

$$G_{\text{vac}}$$

$$=-\left(c_{\text{vac}}-c_{\text{vac}}^{\text{qh}}\right)k_{\text{B}}T,$$

$$\mu (\text{C})=\{\begin{array}{ll}{G}_{\text{diamond}}(T)+\left(H_{\text{graphite}}-H_{\text{diamond}}\right)\text{par}& 0\text{K}<T\le T_{\text{stn}}\\ {\sum }_{-3\le i\le 2}{a}_i(T^i-T_{\text{stn}}^i)+a_3[T\text{ln}(T)-T_{\text{stn}}\text{ln}(T_{\text{stn}})]+G_{\text{diamond}}(T_{\text{stn}})+\left(H_{\text{graphite}}-H_{\text{diamond}}\right)\text{par}& T_{\text{stn}}<T\le T_m.\text{par}\end{array}$$

$$a_{[-3,3]}=\{1.2\times 10^{10},-2.643\times 10^8,2562600,-17369,170.73,-4.723\times 10^{-4},-24.3\}.$$

$$\tau (T,V,\lambda )\equiv {\tau }_{\Delta U,\text{int}}.$$

$${\tau }_{\Delta U,\text{int}}=\left(\frac{1}{2}+\sum _{k=1}^{N}A(k)\right)\Delta t,$$

$$A(k)=\frac{⟨\Delta U_i\Delta U_{i+k}⟩-⟨\Delta U_i⟩⟨\Delta U_i⟩}{⟨\Delta U_i^2⟩-⟨\Delta U_i⟩⟨\Delta U_i⟩}.$$

Full text

anharmonicity, free energy, thermodynamic integration, DFT, MEAM, vacancy concentration, ultra-high temperature, substoichiometric zirconium carbide

I Introduction

Thanks to recent advances in computational thermodynamics, the thermal properties of metals such as aluminium and gold have been investigated up to the melting point, using thermodynamic integration (TI) with Langevin dynamicsGrabowski et al. (2009, 2015); Glensk et al. (2014) based on density functional theory (DFT). A two-step TI approach increases computational efficiency further, making predictions possible for more complex materials at the DFT level of theory. Examples so far include the ultra-high-temperature ceramics ZrC ($T_{m}=3700\,\text{K}$) and HfC ($T_{m}=4160\,\text{K}$),Duff et al. (2015a, 2018) and recent attempts at tackling the emerging class of multicomponent systems.Ikeda et al. (2018); Grabowski et al. Such calculations are not yet routine, but the course of our research and the recent methodological developments of others(Wu and Wentzcovitch, 2009; Moustafa et al., 2017a; Purohit et al., 2018) in this field is in that direction. Here we present some new developments that are a step towards the goal of routinely computing accurate free energies for hard matter systems, including binaries, ternaries, and high-entropy alloys, across the range of temperatures, pressures and chemical potentials, up to and eventually beyond the melting point.

In this work we compute the concentration of vacancies in ZrCx and associated ambient pressure thermodynamics for small deviations from stoichiometry. The ZrCx free energy and derivatives are analyzed in terms of the contributions(Zhang et al., 2018a)

[TABLE]

in which $E_{0}$ is the DFT energy of a static lattice at $T=0\,$K, $F_{\text{el}}$ is the Helmholtz free energy contribution from the thermal excitations of electrons, $F_{\text{qh}}$ is the quasiharmonic vibrational contribution, $F_{\text{ah}}$ is the anharmonic vibrational contribution, $F_{\text{el-vib}}$ is the electron-vibration contribution, and $F_{\text{config}}$ is the contribution of configurational entropy due to the number of distinct point-defect distributions. For each of the five temperature-dependent terms we have calculated the dependence on the independent variables volume and temperature up to the melting point. Considerable attention in this paper is given to the method we use to compute the challenging anharmonic term, $F_{\text{ah}}$. The method described achieves effective DFT accuracy in $F_{\text{ah}}$ (1 meV per bulk atom) with only an order of magnitude greater computational cost than ordinary quasiharmonic free energy calculations.

This paper is set out as follows. Sec. II.1 gives the context of our approach, theoretical details are in Sec. II.2, a description of the modified embedded atom method (MEAM) potential fitting, which reduces overall the number of expensive DFT calculations, in Sec. II.3, thermodynamic integration in Sec. II.4, and DFT technical details in Sec. II.5. Benchmarking is described in terms of accuracy and precision in Sec. III.1 and computational efficiency in Sec. III.2. Application to ZrCx provides insight into the nature of anharmonicity in substoichiometric binary crystals in Sec. IV.1, prediction of vacancy concentration in Sec. IV.2, and analysis of ZrCx heat capacity in Sec. IV.3.

II Methods

II.1 Background

There are a number of approaches in the literature to calculate the anharmonic vibrational properties of crystals,Alfè et al. (2001, 2002a, 2002b); Ackland (2002); Duff et al. (2015a); Errea et al. (2011); Glensk et al. (2014); Grabowski et al. (2009, 2015); Hellman et al. (2013); Monserrat et al. (2013); Klein and Horton (1972); Zwanzig (1954); Frenkel and Ladd (1984); Moustafa et al. (2017b) including thermodynamic integrationKirkwood (1935) (TI), which is the method used in this work. In TI the anharmonic part of the full Hamiltonian, $E-E_{\text{qh}}$, is switched on with the parameter $\lambda\in[0,\,1]$, in this instance linearly as $E_{\text{mix}}(\lambda)=E_{\text{qh}}+\lambda(E^{\text{}}-E_{\text{qh}}^{\text{}})$. Classical averages of $\partial_{\lambda}E_{\text{mix}}(\lambda)$ are obtained stochastically from molecular dynamics (MD), and numerically integrated along the coupling path to give the free energy due to $E_{\text{ah}}$:

[TABLE]

Note, $\partial_{\lambda}E_{\text{mix}}(\lambda)=E(\mathbf{R},\,V)-E_{\text{qh}}(\mathbf{R},\,V)$, where $E(\mathbf{R},\,V)$ is the full potential energy surface and $E_{\text{qh}}(\mathbf{R},\,V)$ the volume-dependent harmonic potential energy surface in Born-Oppenheimer nuclear coordinates $\mathbf{R}$.

The partitioning of the vibrational free energy $F_{\text{vib}}=F_{\text{qh}}+F_{\text{ah}}$ divides the problem conveniently into a simple quantum mechanical part, in which the vibrations are quantised as phonons, and an anharmonic part in which the vibrations are treated classically. Thus $F_{\text{vib}}$ has the appropriate low-temperature quantum statistics. The anharmonicity is treated classically but also non-perturbatively, which is important at high temperatures as the melting point is approached. In order to evaluate the anharmonic term, the expectation values $\left\langle\partial_{\lambda}E_{\text{mix}}(\lambda)\right\rangle_{\lambda}$ require between $10^{3}\ldots 10^{7}$ configurations for a typical $\lambda$-ensemble at a typical supercell size. To produce a free energy surface $F_{\text{ah}}(V,\,T)$ one must sample ensembles across dimensions of strain (here volume), temperature and coupling parameter, $N_{\lambda}\times N_{V}\times N_{T}\approx 10^{2}\ldots 10^{3}$. Thus the ball-park number of total energy calculations, between $10^{5}\ldots 10^{10}$ configurations, is prohibitive at the highly-converged DFT level of accuracy required.

In one approach to reduce computational complexity, $F_{\text{ah}}$ is obtained by cumulating a sequence of thermodynamic integrations. In an implementation of this approach referred to as TU-TILD,Duff et al. (2015a) which is expressed by Eqn. (3), much of $F_{\text{ah}}$ is captured using an inexpensive MEAM potential. This results in faster convergence of the expensive TI from MEAM to DFT, expressed in the last term of Eqn. (3).

[TABLE]

In practice to save computation time, the DFT MD calculations in a TU-TILD procedure were usually performed with a low-converged expansion of the wavefunctions, using a reduced number of plane-waves, and fewer k-points than required for maximum accuracy. The maximum accuracy was then obtained by up-sampling, as in the original UP-TILD methodGrabowski et al. (2009). The methodology we introduce below, inspired by these approaches, was devised in order to make significant further savings in computation time without sacrificing accuracy.

II.2 MEAM thermodynamic integration approach

The approach we propose in this work can be summarized by

[TABLE]

Our approach calculates the anharmonic free energy of a MEAM crystal referenced to a harmonic DFT crystal, which is formally the first stage in Eqn. (3). In the present method the quasiharmonic Helmholtz free energy at each volume is still explicitly represented by the volume-dependent dynamical matrix calculated with DFT, which captures much of the thermal expansion, but the anharmonic terms are now entirely described by the MEAM potential. The success of the method depends on being able to generate a MEAM potential of sufficient accuracy to replace the anharmonic DFT contribution. It is by no means obvious a priori that this is possible, or if it is, that the process of generating the potential is not too expensive to warrant the effort.

From the potential terms in Eqn.(4) it is clear that $F_{\text{ah}}^{\text{}}$ can be evaluated by this method to a high level of precision using modest computational resources, but the MEAM TI approximation introduces systematic potential errors with respect to DFT TI. Accuracy must be carefully controlled by generating custom MEAM potentials from high-quality DFT MD. Generating the training and validation data is time consuming, so fitting the potentials becomes the primary computational cost in predicting $F_{\text{ah}}$ in our TI approach. These costs incurred before doing the MEAM TI will be shown in Sec. III.1 to be comfortably small enough. Details of potential fitting and error control are presented in the following section.

II.3 Potential fitting

Interatomic potentials have been fitted using the reference-free modified embedded atom method (RF-MEAM).Duff et al. (2015b) The MEAM potentials fitted in this work lack transferability and are specialized to perform for the intended application. For instance, to model $F_{\text{ah}}(V,\,T)$ we fit a separate potential at each volume considered, which minimizes the possibility of strain-dependent errors. Obtaining the correct implicit volume dependence of a potential is important, as explicit anharmonic effects depend sensitively on the degree of lattice expansion, as demonstrated in Sec. IV.1.

Potentials for Zr32C32 and Zr32C31 have been fitted at the cell lattice parameters $a=\{4.685,\,4.730,\,4.759,\,4.801,\,4.850\}$ Å. In crystals of lower than cubic symmetry, thermal expansion may involve other modes of strain, but in our case only volumetric strain need be considered. At each volume the potential is fitted to configurations sampled from DFT MD runs between $T=200$ K and $T=3800$ K. The fitting set for each volume comprises $10^{3}$ Zr32C31 configurations, which supply energies and forces to the objective function, for minimizing force residuals and energy-residual variances. Fitting to the forces on each atom allows the potential to be specified by fewer distinct configurations than fitting to energies, as $N_{\text{forces}}=3N_{\text{atom}}N_{\text{energies}}$. To generate the fitting data points efficiently, low-quality DFT can be up-sampled to high-quality DFT to produce high-quality target forces and total energies.

The interatomic potentials are generated using a genetic-algorithm conjugate-gradient fitting procedure implemented and publicly available in the MEAMFIT2 code.Duff et al. (2015b); Duff (2016); MEA The code fits an RF-MEAM potential that permits locally positive and negative density terms in order to increase variational freedom, subject to a net positive background. The fitted potential has 3 embedding and 3 pairwise terms, within a radial cutoff of $4.8$ Å, which includes interactions up to third nearest-neighbor. These potential parameters provide a satisfactory compromise between accuracy and complexity, in terms of minimizing residual variances on hold-out data using the fewest degrees of freedom (78 parameters for the 3-3 potential). The quality of fits for energy and forces is presented in Fig. 1 for Zr32C31.

In this RF-MEAM application a different potential is fitted at each volume but we require a potential to be transferable across composition, i.e. we want the same potential to describe both Zr32C31 and Zr32C32, for a given volume at any temperature. This transferability ensures a systematic error cancellation in $F_{\text{ah}}$ for Zr32C31 and Zr32C32 that helps in calculating accurate fully anharmonic vacancy formation energies.

II.4 Thermodynamic integration

$F_{\text{ah}}$ is estimated by computing $\left\langle\Delta E\right\rangle_{\lambda}$, with $\Delta E=E^{\text{MEAM}}-E^{\text{DFT}}_{\text{qh}}$, for a series of $11$ ensembles at equal increments of $\lambda$, $\lambda_{i}=i/10$. In Fig. 2 we show the dependence of $\left\langle\Delta E\right\rangle_{\lambda}$ on $\lambda$ across a series of temperatures. We see that at each temperature, $\partial_{\lambda}\langle\Delta E\rangle_{\lambda}\leq 0$, a necessary condition that is easy to prove, as in a derivation of the Gibbs-Bogoliubov inequality. In Fig. 2 inset the convergence of $\left\langle\Delta E\right\rangle_{\lambda}$ is shown for the first 60,000 time steps of a simulation. Expectation values are generated using Langevin MD, with a one femtosecond time-step,Duff et al. (2015a) and a friction parameter of $\gamma=0.05$ fs-1 for Zr32C31 and $\gamma=0.01$ fs-1 for Zr32C32. At each pair, $\{V_{i}\,T_{i}\}$, the expectation value $\left\langle\Delta E\right\rangle_{\lambda}$ is fitted in $\lambda$ by least squares to the truncated power series

[TABLE]

for which the coefficients are alternating in sign and converging. If intrinsic defects form and migrate on the atomic vibration time-scale this series is expected to poorly converge. In this work we exclude any system configurations in which Frenkel defects have spontaneously formed, for example at the melting point, in order to ensure well-converged thermodynamic integrations of the anharmonic free energy of defect-free ZrCx.

Errors in predicting $F_{\text{ah}}$ are considered from two primary sources, namely statistical convergence and a systematic potential error. The DFT benchmark also has a small convergence error which is accounted for, but other sources of error, such as from DFT exchange-correlation, electron-phonon scattering, and other quantum effects beyond the harmonic approximation, are beyond the scope of this paper. The three countable contributions are shown schematically in Fig. 3, and give the total expected error of

[TABLE]

Precision error $\delta_{\text{rand}}$ arises from evaluating an observable from a finite number of samples in the MEAM MD, and the systematic error $\delta_{\text{sys}}$ is due to the energy difference between a MEAM potential and DFT. The convergence error in the benchmark $F_{\text{ah}}$ value from TU-TILD is $\delta_{\text{ref}}=0.1$ meV/at.Duff et al. (2015a)

The statistical convergence $\delta_{\text{rand}}$ is computed using stratified systematic sampling.Haile (1997) In the simulation the precision error scales asJanke (2002)

[TABLE]

where $N_{\lambda}$ is the number of integration path points sampled (with $N_{\lambda}=11$ in our case), $\sigma_{\lambda}$ is the norm of $\Delta E$ standard deviations, $\tau$ is the $\Delta E$ autocorrelation time (ca. $11$ fs, see Appendix), and $t$ is the simulation time (ca. $0.1$ ns). We emphasise that the TI method described differs from other approaches in that statistical convergence is not accuracy-limiting, for example, nanosecond simulations can comfortably be performed in a day on a low-performance computing platform.

The systematic potential error in the anharmonic free energy $\delta_{\text{sys}}$ can be computed by thermodynamic integration

[TABLE]

$\delta_{\text{sys}}$ is the primary error source in the method we describe to compute $F_{\text{ah}}$. Accuracy benchmarks in Sec. III.1 show $\delta_{\text{sys}}$ can be sufficiently controlled to satisfy 1 meV/at bulk convergence across the $F_{\text{ah}}(V,\,T)$ surface.

II.5 Technical details

Periodic plane-wave DFT calculations were performed using the VASP software,Kresse and Furthmüller (1996a, b) with the local density approximation (LDA) exchange-correlation function.Perdew and Zunger (1981) The projector-augmented wave (PAW) method was used,Kresse and Joubert (1999) with $4s$- and $4p$-Zr electrons included as valence states.

$E_{0}$ was computed on a mesh of 11 volumes, and at each volume the internal coordinates have been relaxed to give residual forces under $10^{-6}$ eV/Å. Self-consistent field (SCF) total energies and energy eigenvalues have been resolved to $10^{-9}$ eV. Methfessel-Paxton smearing has been used with a width of 0.1 eV.Methfessel and Paxton (1989) The kinetic energy cutoff has been set to 700 eV and k-point mesh $12\times 12\times 12$ for the $2\times 2\times 2$ supercell. The $E_{0}$ vacancy formation energy contribution is extrapolated to the dilute limit, using data points from Zr32C31, Zr108C107 and Zr256C255.

For the quasiharmonic Helmholtz free energy $F_{\text{qh}}$, the kinetic energy cutoff has been set to 700 eV and k-point mesh to $6\times 6\times 6$ for the $2\times 2\times 2$ supercell. Phonons were calculated using the small displacement supercell method with the PHONOPY code.Togo and Tanaka (2015) At each of the 11 $2\times 2\times 2$ supercells that span the range of lattice parameters $[4.575,\,4.875]$ Å, sets of 18 displacements were made for Zr32C31 and sets of four displacements for Zr32C32. The phonon q-points were sampled by a mesh of $25\times 25\times 25$ points for the $2\times 2\times 2$ supercells.

The electronic Helmholtz free energy $F_{\text{el}}(V,\,T)$ has been calculated using the Mermin finite-temperature formulation of DFT,Mermin (1965) on a mesh of 10 temperatures and 8 volumes sampled between $V_{\text{eq}}(T=0\,\text{K})$ and $V_{\text{eq}}(T_{m})$. Electron states are self-consistent to at least $10^{-7}$ eV/atom. We used $384$ bands, which was sufficient to span all states with partial occupation up to the melting point $T_{m}$. A kinetic cutoff energy of 700 eV was used, with k-point sampling at $12\times 12\times 12$ for the $2\times 2\times 2$ supercells.

The electron-vibration Helmholtz free energy $F_{\text{el-vib}}(V,\,T)$, has been calculated from low-converged MD configurations that are subsequently up-sampled, as in the proceedure recently performed for a number of transition metals.Zhang et al. (2017) The electronic free energy is calculated for each MD configuration, using the Mermin formulation at electronic temperature corresponding to the MD ensemble temperature. At each volume-temperature mesh point, the electronic free energies are averaged over the ensemble configurations, and referenced to the perfect crystal, in order to find the electron-vibration coupling contribution to the Helmholtz free energy.

The anharmonic Helmholtz free energy $F_{\text{ah}}(V,\,T)$ was determined using a mesh of six temperatures and five volumes. Temperatures span $0\,$K to $T_{m}$, and volumes $V_{\text{eq}}(T=0\,\text{K})$ to $1.15V_{\text{eq}}(T_{m})$. Potentials were fitted to MD configurations from DFT that used a 700 eV cutoff and k-point sampling mesh of $6\times 6\times 6$ for the $2\times 2\times 2$ supercell.

III Benchmarks

III.1 Accuracy and precision

Our MEAM TI approach predicts $F_{\text{ah}}(V,T)$ to within a target accuracy of 1 meV/at compared to DFT endpoint TI. This is demonstrated in Fig. 4a. Perfect Zr32C32 is shown with error bars (TU-TILD reference) for 25 volumes and temperatures up to the melting point. $F_{\text{ah}}(V,T)$ energies are converged to sufficient precision that the error bars are in effect systematic potential error bars. The mean absolute error (MAE) is 0.5 meV/at, with a mean signed deviation of 0.05 meV/at. The MAE values at the lattice parameters $\{4.685,\,4.730,\,4.759,\,4.801,\,4.850\}$ Å are $\{0.72,\,0.38,\,0.35,\,0.46,\,0.64\}$ meV/at. Errors resolved at the temperatures $\{760,\,1900,\,2500,\,3200,\,3805\}$ K have the MAE values $\{0.39,\,0.48,\,0.43,\,0.59,\,0.66\}$ meV/at.

On the basis of adequately small errors for bulk ZrC, we propose using the MEAM thermodynamic integration approach for more complex systems. In this regard Zr32C31 is a useful test case for two reasons. The carbon-vacancy introduces complexity in terms of physical interactions. It removes inversion symmetry at sites around the vacancy, so there are terms in the energy of odd-order in atomic displacements, previously excluded by symmetry in perfect ZrC. Secondly, making free energy predictions per vacancy increases computational complexity considerably due to the nature of statistical error scaling for TI predictions on a per-vacancy basis.

For the vacancy system Zr32C31, $F_{\text{ah}}(V,T)$ is shown in Fig. 4b. Obtaining comparable DFT TI values for systems with vacancies like Zr32C31 is prohibitively expensive in general but we have computed a DFT benchmark for Zr32C31 at $a=4.801$ Å and $T=3200$ K. The MEAM thermodynamic integration is found to overestimate the DFT TI reference $F_{\text{ah}}$ value by $0.4$ meV/at, which is comparable to the MAE in the perfect bulk ZrC.

III.2 Computational cost

In Table III.2 timings are reported for the MEAM-based TI in this work and TU-TILD (DFT) calculations. Both methods compute the $F_{\text{ah}}(V,\,T)$ surface across 25 mesh points for a Zr32C32 test case. The MEAM approach does not have DFT TI and DFT TI up-sampling steps, which account for the majority of the TU-TILD $F_{\text{ah}}$ cost. In the MEAM approach, the main CPU-time overhead is high-quality DFT calculations on selected MD configurations, which are used create data to fit the MEAM potentials. Furthermore the time required to optimizeDuff (2016) the MEAM potential with a large fitting set is substantially longer (ca. $\times 10$), compared to a MEAM potential used in the intermediate TI steps in TU-TILD. Despite this the former scheme still gains a factor of at least $\times 10$ in efficiency overall due to having no DFT TI or DFT TI up-sample. For anharmonic predictions where 1 meV/atom convergence is sufficient, the TI method described in this work is likely to be a cost effective choice for metals and alloys. It would be of interest to compare the efficiency of less specialised machine learned potentials(Cheng et al., 2019; Behler, 2017; Grabowski et al., ) to the MEAM type as applied here, in terms of parameter fitting time, required training DFT data and potential compute time.