PHYSICAL REVIEW B 108, 184107 (2023) Dynamically stabilized phases with full ab initio accuracy: Thermodynamics of Ti, Zr, Hf with a focus on the hcp-bcc transition Jong Hyun Jung , Axel Forslund , Prashanth Srinivasan ,* and Blazej Grabowski Institute for Materials Science, University of Stuttgart, Pfaffenwaldring 55, 70569 Stuttgart, Germany (Received 29 August 2023; accepted 30 October 2023; published 27 November 2023) Certain systems feature phases that are dynamically unstable at 0 K but are stabilized by vibrations at higher temperatures. Treatment of these phases by conventional 0-K methods is not feasible and effective harmonic models introduce approximations. Here, we significantly advance the direct upsampling methodology [npj Comput. Mater. 9, 3 (2023)] to obtain free energies including the anharmonic contribution to full ab initio accuracy also for such dynamically stabilized phases. The centerpiece behind the procedure is accurate machine- learning potentials (moment tensor potentials) which are used to efficiently scan the volume-temperature space to uncover the stability regime and to perform thermodynamic integration on a dense grid within the stable window. We apply the methodology to the prototype systems Ti, Zr, and Hf and calculate hcp-bcc transition properties and thermodynamic properties of both phases. We find a very good agreement for the heat capacities with existing experimental/CALPHAD data, and an overall best agreement for Ti. The transition properties agree well on a relative temperature axis, where the temperature is scaled with respect to the transition temperature. Anharmonic free energies increase the transition temperature by up to one thousand kelvin. Electronic effects are smaller and bring down the transition temperature by as much as 172 K. We establish a new definition of the 0-K energy-volume curve for the dynamically stabilized bcc phase. Instead of using static lattice ab initio values, an extrapolation of the high-temperature high-accuracy free-energy surface to 0 K provides a physically more meaningful description. With this effective 0-K definition, discrepancies existing in the literature between CALPHAD and ab initio values are addressed. DOI: 10.1103/PhysRevB.108.184107 I. INTRODUCTION Thermodynamic properties such as the Gibbs energy, heat capacity, expansion coefficient, and bulk modulus are crucial for materials design [1–3]. The Gibbs energy for example determines phase stability and facilitates the prediction of phase diagrams. To calculate thermodynamic properties with ab initio accuracy, a precise representation of the free-energy surface is needed. Two factors are crucial in this regard—a very dense volume-temperature grid on which free-energy calculations are performed, and the inclusion of the rele- vant physical excitation mechanisms. Recently, an efficient methodology—direct upsampling—was introduced to per- form such high-accuracy calculations [4,5]. However, the prototype systems in these previous works were single-phase elements with the considered phases being thermodynami- cally and dynamically stable across the entire temperature range from 0 K up to the melting temperature. Severe challenges arise for systems featuring a high- temperature phase that is dynamically unstable1 at 0 K [6–8]. *prashanth.srinivasan@imw.uni-stuttgart.de, prashanth.s.iitm@gmail.com 1Conventionally, dynamically unstable phases are those that corre- spond to a saddle point on the potential energy surface and thus show imaginary 0-K phonon frequencies. We consider here such phases with the additional requirement that they are stabilized by thermal While the phase is stabilized at higher temperatures through anharmonic vibrations, it inevitably transforms into a differ- ent, low-symmetry phase at lower temperatures. Prototype systems are the elements in group IV of the periodic table— titanium (Ti), zirconium (Zr) and hafnium (Hf)—which are dynamically stable in the body-centered cubic (bcc) phase at high temperatures and transform into the hexagonal ω or hexagonal close-packed (hcp) phase at lower tempera- tures depending on the boundary conditions of the simulation box [9,10]. Conventional 0-K approximations cannot be ap- plied to dynamically unstable phases, since it is impossible to obtain free energies from imaginary phonon frequencies. Methods exist in literature that enable the treatment of dynamically unstable phases. Large-displacement phonons were used to describe the bcc phase of Ti, Zr and Hf [11]. Quasiharmonic Debye-Grüneisen models were applied to bcc Ti [12,13]. The phase transition in Zr was modeled with a temperature-dependent effective potential (TDEP) [14,15], and using machine-learned force fields [16]. A piecewise polynomial potential partitioning was applied to describe the local minima near the equilibrium bcc configuration [17]. Table I lists a compilation of relevant works where thermo- dynamic properties of the prototype systems were calculated. vibrations before melting sets in. Since we treat these phases pri- marily in their metastable and thermodynamically stable temperature regime, we mostly refer to them as “dynamically stabilized” phases. 2469-9950/2023/108(18)/184107(22) 184107-1 ©2023 American Physical Society https://orcid.org/0000-0002-2409-975X https://orcid.org/0000-0002-0419-3546 https://orcid.org/0000-0002-9199-4340 https://orcid.org/0000-0003-4281-5665 http://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevB.108.184107&domain=pdf&date_stamp=2023-11-27 https://doi.org/10.1038/s41524-022-00956-8 https://doi.org/10.1103/PhysRevB.108.184107 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) TABLE I. Ab initio studies of the isobaric heat capacity (Cp), the thermal expansion coefficient (α) and the bulk modulus (B) for hcp and bcc Ti, Zr, and Hf. The employed exchange-correlation (XC) functionals and the potentials representing the core electrons are listed. The columns under “Contributions to F” indicate which terms were included to the total free energy [0-K static energy (E0K ), electronic (el), quasiharmonic (qh), and anharmonic (ah)]. For works that included the anharmonic contribution, the number of volume V (or pressure P) and temperature T data points is mentioned. Studies including the explicit anharmonic contribution up to all orders and to the accuracy of DFT are printed in boldface. All of the used abbreviations are expanded in Ref. [19]. DFT methodology Contributions to F Year Ref. Elements XC Potential E0K el qh ah Grid for ah Cah p αah Bah 2007 [20] hcp Ti LDA, PBE PAW x x x 2007 [21] hcp Ti, hcp Zr LDA NC x x x 2008 [22] hcp Zr GGA PAW x x x 2009 [12] hcp, bcc Ti PBE PAW x x x, bcc: qh Debye 2010 [13] hcp, bcc Ti PBE PAW x x x, bcc: qh Debye 2011 [23] hcp, bcc Zr PBE PAW x x x, bcc: qh Debye 2012 [11] bcc Ti, Zr, Hf - PAW x large displacement 2012 [24] hcp, bcc Zr PBE PAW x qh Debye 2013 [25] hcp Hf PBE PAW x x x 2017 [17] bcc Ti PBE PAW x x x piecewise polynomial 1P × 7T x potential partitioning 2019 [26] hcp, bcc Hf PBE PAW x x x, qh Debye 2020 [27] bcc Ti PBE PAW x x AIMD, 1P × 4T x x stress-strain curve 2022 [28] bcc Ti PBE PAW x large displacement 2023 this work hcp, bcc Ti, Zr, Hf PBE PAW x x x direct upsampling � 12V × 7T x x x While such effective models typically provide access to good renormalized phonons, they are approximate in nature with respect to the fully anharmonic ab initio description, especially at high temperatures [4,5,18]. The methods are hence inadequate for accurate prediction of thermodynamic and phase transformation properties, which are strongly af- fected even by small changes in the free energy. Thus, further development of computational methods that take full anhar- monicity into account is needed. In the present work, we address this challenge and sig- nificantly advance the direct upsampling methodology to be applicable also for dynamically stabilized phases. We demonstrate the performance of the method on Ti, Zr, and Hf. The key ingredient of the methodology is optimized machine learning potentials, specifically, moment tensor po- tentials (MTPs) [29,30]. For each element, we train accurate MTPs to high temperature density-functional-theory (DFT) data of the hcp and bcc phase. The MTPs are used to obtain stability maps, which are crucial in determining the volume- temperature (V, T ) grid for the free-energy calculations. For the bcc phase, the stability window is restricted and requires a rigorous adjustment of the (V, T ) grid. At lower temperatures, there is a constraint due to the transformation to the hexagonal ω phase [9], whereas at high temperatures, Frenkel defects oc- cur. Within the stability regime, careful treatment is required to fit an effective harmonic potential for the dynamically sta- bilized phase, and utilize it as a reference for thermodynamic integration. In addition to the accurate consideration of the full vibrational and electronic contributions, we also propose a novel approach to include the effect of vacancies for the dynamically stabilized phase. Utilizing the methodological advancements, we compute accurate thermodynamic properties of the hcp and the bcc phases and hcp-bcc transition properties for Ti, Zr, and Hf, and compare the DFT prediction to data from experiments and from the CALculation of PHAse Diagrams (CALPHAD) method. Based on the accuracy of the free-energy surface in the high-temperature regime, we establish a new understand- ing of the 0-K behavior for the dynamically stabilized phase. The introduced effective energy-volume curve corresponds to the extrapolation of the high-accuracy free-energy surface to 0 K and provides a physically consistent description of the 0-K energies. The resulting effective enthalpies resolve the long-standing mismatch in the comparison of T = 0 K DFT static lattice data and CALPHAD [6,31,32]. II. METHODOLOGY We apply the direct upsampling methodology separately to the hcp and bcc phases of Ti, Zr, and Hf, with refine- ments mainly for the dynamically stabilized bcc phase. Within the direct upsampling methodology, the free energy for each phase is adiabatically decomposed as [33] F (V, T ) = E0K(V ) + F el(V, T ) + F qh(V, T ) + F ah(V, T ). (1) Here, E0K(V ) represents the energy-volume (E -V ) curve at 0 K, F el(V, T ) the electronic free energy (including coupling 184107-2 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) to vibrations), F qh(V, T ) the quasiharmonic free energy, and F ah(V, T ) the explicitly anharmonic contribution. The specifics of the direct upsampling methodology and the computational details to obtain each of the free-energy contributions and, from F (V, T ), the thermodynamic prop- erties, were described in detail in Refs. [4,5] (cf. also the respective Supplementary Information). The main steps of the methodology comprise the training of the MTPs, ther- modynamic integration from the quasiharmonic reference to the MTP, and a two-stage upsampling (free-energy pertur- bation [34]) to DFT. The final DFT energy includes the electronic contribution calculated from the finite-temperature formalism of DFT [35]. In practice, the molecular dynamics (MD) runs during thermodynamic integration are performed with LAMMPS [36], and the DFT runs during upsampling using VASP [37] within the projector augmented-wave (PAW) method [38,39]. For the exchange-correlation functional, we use here the generalized gradient approximation (GGA) by Perdew, Burke, and Ernzerhof (PBE) [40]. For the hcp phase we use a fixed c/a ratio determined at 0 K. Further computa- tional details and parameters can be found in Ref. [19]. For systems exhibiting dynamically stabilized phases, re- finements to the previously published procedure [4,5] are required, as described in the following subsections. The MTPs (Sec. II A) carry out additional tasks in the refined method- ology, in particular, the efficient detection of the stability regime in which the free-energy calculations should be per- formed (Sec. II B). The fit of the effective quasiharmonic reference, which is used to calculate F qh(V, T ) and which serves as a reference for calculating F ah(V, T ), needs a careful adjustment in the stability regime (Sec. II C). The conven- tional 0-K energy-volume curve becomes meaningless for the dynamically stabilized phase, while the high-temperature free-energy surface leads to a better definition of an effective 0-K E -V curve (Sec. II D). Lastly, a new approach is nec- essary to analyze and incorporate the effect of vacancies on the bulk thermodynamics of the dynamically stabilized phase (Sec. II E). A. MTPs: The backbone of direct upsampling The main ingredient that accelerates the direct upsam- pling workflow is accurate interatomic potentials. In line with previous works on direct upsampling [4,5,41], we use machine-learning interatomic potentials—the moment tensor potentials (MTPs)—also in this work. For systems exhibit- ing dynamically stabilized phases, the MTPs serve multiple purposes in the overall methodology. Firstly, they are vital in identifying the stable temperature regime of the individ- ual phases, as will be discussed in Sec. II B. Secondly, the accuracy of the MTPs decides the efficiency of direct upsam- pling. In Ref. [4], the number of snapshots needed for the convergence of the perturbative upsampling was formalized to be quadratically dependent on the energy root-mean-square error (RMSE) of the MTP. Thirdly, the MTPs are pivotal in the efficient computation of the vacancy contribution as elaborated in Sec. II E. Hence, a well-trained MTP is crucial in the overall performance of the methodology. For each element, we train three different potentials—a separate MTP for each individual phase, and an MTP trained simultaneously on both phases. The training set contains 96 and 128-atom configurations in the hcp and bcc phase, respec- tively, for a set of eight volumes at the experimental melting point of the element which are denoted by the green dots in Fig. 1. For the phase-specific MTPs, 720 configurations of the target phase are chosen for the training set. For the combined MTP, we merge both individual training sets leading to 1440 configurations in total. In each case, we train MTPs with dif- ferent levels using the MLIP code [29,30]. The corresponding training and test errors as a function of the level of the MTP are plotted in Fig. S3 in Ref. [19]. For estimating the sta- bility regime and performing thermodynamic integration, we choose the level 24 MTPs with 913 fitting parameters (level 22 for bcc Hf), above which there is no significant improve- ment in the accuracy of the MTPs. Further details regarding the training set, fitting weights, and the fitting procedure are provided in Ref. [19]. Figure 2 highlights the accuracy of the MTPs. The MTP- RMSEs [cluster of symbols on the bottom left of Fig. 2(a)] are contrasted with RMSEs of conventional classical potentials [embedded atom method (EAM)/modified embedded atom method (MEAM)] from literature for Ti, Zr, and Hf in both of the relevant phases, hcp and bcc [42–44]. The RMSEs in energies and forces of the MTPs are less than 2.5 meV/atom and 0.11 eV/Å, respectively, almost an order of magnitude better than most EAM/MEAM potentials. In Fig. 2(b), we zoom-in on the RMSEs of the MTPs (including, additionally, previously fitted MTPs for refractory unaries [5]). One needs to be aware that the range of values on the x and y axes is much smaller than in the upper plot in Fig. 2(a). The slightly higher force RMSE for the bcc phase comes from the higher temperature of the snapshots chosen for the test data [the highest temperature in the (V, T ) grid]. We observe that, for all elements, the MTP trained on both phases is almost as accurate as the MTPs trained separately on each phase. For achieving the same target accuracy of < 1 meV/atom in the free energy at a given (V, T ), a dif- ference in the energy RMSE of 0.62 meV/atom (between the Zr-hcp and Zr MTPs for instance) would require ≈ 10 additional snapshots. The single MTP fitted to hcp and bcc (for Zr) also has a higher force error by 0.011 eV/Å than the Zr-hcp MTP, leading to a difference of ≈ 5% in the diffu- sion probability (red contours in Fig. 1). Likewise, the bcc-ω transformation temperature (blue lines in Fig. 1) is lower by 50-100 K using the MTP trained on both phases, compared to the pure bcc MTP. These differences should be consid- ered in the prediction of the stable temperature regime of the individual phases. However, the final free-energy surface prediction does not get affected within the stability regime of the individual phases. B. Stability regimes of the hcp and bcc phases A dense and optimized (V, T ) sampling grid, positioned appropriately based on the stability of a particular phase, is crucial to enable a reliable parametrization of the free-energy surface. For the systems studied earlier in Refs. [4,5], the grid was chosen with the temperature ranging from 0 K to above the melting point, since these systems (e.g., Al, Mg, V, Mo) exhibit a single stable phase across the entire temperature 184107-3 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) FIG. 1. Free-energy surfaces illustrating various quantities of relevance to the direct upsampling methodology: (1) the (V, T ) grids represented by the crossings of the gray lines; (2) the probabilities for Frenkel defect formation and diffusion by the red contours; (3) the transformation temperatures from bcc to ω by the blue lines including an uncertainty estimate as given by the blue numbers; (4) the thermal expansion at ambient pressure by the light green lines; (5) the PBE-predicted hcp-bcc transition temperatures T hcp→bcc PBE by the filled black circles; (6) the experimental hcp-bcc transition temperatures T hcp→bcc exp by the hollow black circles; (7) the volumes at T melt exp chosen for the MTP training by the green dots, and (8) the volumes and temperatures for the low- or intermediate-temperature effective quasiharmonic fitting by the green crosses. The probabilities for Frenkel defect formation correspond to counts of the formation of such defects in 10 ps long MTP MD runs (in an initially defect-free bulk cell) on a much finer (V, T ) grid than the one used for the free-energy calculations. The leftmost labels hcp and bcc indicate hexagonal or cubic cell shapes, respectively, and the labels (ω, bcc, hcp) on the surfaces indicate the stability regimes under the constraint of the specific cell shape. range. Determination of the volume-temperature regime for dynamically stabilized phases is more challenging. We ad- dress this challenge by performing MD runs with the MTPs to rapidly scan the relevant (V, T ) space and thereby locate the stable regime. We execute this procedure separately for the low-temperature hcp and the dynamically stabilized bcc phase. The resulting (V, T ) meshes optimized for the explicit free-energy calculations within direct upsampling are marked by the gray grids in Fig. 1. The hcp phase is dynamically stable at 0 K for the here considered elements (Ti, Zr, Hf) and the temperature regime for the free-energy calculations thus starts at 0 K. As the highest temperature for the hcp free-energy calculations, we choose a value that is a few hundred kelvin above the ex- perimental hcp-bcc phase transition temperature (marked by the hollow circles). Since the eventually computed PBE tran- sition temperatures (filled black circles) are lower than the experimental ones, the choice of the highest temperature for the grid adequately covers the thermodynamic stability region of hcp also including a decent extension into the metastable regime. As for the volume range of the grid, we start from slightly below the 0-K equilibrium volume and extend it well beyond the high-temperature equilibrium volume (cf. the light green line indicating the thermal expansion in Fig. 1). Note that an optimized volume-temperature grid may require a self-consistent adaption after a first estimate of the thermal expansion has been obtained. The same applies for the density of the grid points (which do not need to be equidistant). The MTPs are used to efficiently implement such a self-consistent procedure. For comparison with the bcc phase, we extend the analysis of the stability of the hcp phase beyond the grid, in particular up to the experimental melting temperature (at which the MTPs have been fitted). Close to the melting temperature we observe the onset of Frenkel defect formation and dif- fusion (also observed in the bcc phase as discussed below), the probability of which is marked as red lines in Fig. 1. Figure 3 illustrates such defect formation during an MD run at 1665 K in hcp Ti. In this particular scenario, a Frenkel defect is formed and it migrates within the first 30 ps of the MD run. Specifically, the accompanying self-interstitial migrates upward and reaches a basal-octahedral configuration [45] (in- dicated by the red ball and the red bond). The vacancy that is left behind attracts a nearby atom that relaxes into an in- terstitial position. During the next 5 ps of the MD run, defect migration leads to another configuration of the self-interstitial in a split dumbbell along the c axis [45–47]. After additional 5 ps, the self-interstitial and the vacancy reach each other and annihilate so that the Frenkel defect vanishes. Such defects are most commonly observed for Ti and least for Hf within 184107-4 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) FIG. 2. (a) Root-mean-square errors (RMSEs) in energies and atomic forces predicted by the MTPs (level 22-24) trained in this work with respect to ab initio values for hcp and bcc Ti, Zr, and Hf. For comparison, the RMSEs of the following EAM/MEAM potentials (hollow symbols) from literature are also shown: Ti [42], hcp Zr (potential 3 in Ref. [43]), bcc Zr (potential 2 in Ref. [43]), and Hf [44]. (b) Zoomed-in image of the MTP RMSEs from (a). Each group IV element has a set of three MTPs—fitted separately to the hcp and bcc phase, and a single combined MTP (labelled by the element symbol only). For comparison, we also show the RMSEs of previously trained MTPs for bcc refractory elements [5]. The RMSEs are calculated on test sets that contain 96-128-atom supercell configurations at the highest T in the corresponding (V, T ) grid, and at the equilibrium volume. the investigated temperature range. In all cases, however, the grid for the hcp free-energy calculations does not overlap with the diffusion contours, as observed in the top row of Fig. 1. For the dynamically stabilized bcc phase, it becomes more complicated to obtain the stable temperature window for per- forming thermodynamic integration. There is a constraint on both the lower and upper ends of the temperature regime, as shown by the grids in the bottom row of Fig. 1. Below a certain temperature, the bcc phase starts to transform into the ω phase (for the here applied constant volume boundary condition). The onset of the ω-phase formation is denoted by blue lines in Fig. 1. In order to distinguish the ω phase from the bcc phase, we use a structural descriptor that was earlier developed and used in Refs. [9,48] (see also Ref. [19]). In favor of our calculations, the transformation to the ω phase FIG. 3. Frenkel defect formation and kinetics in initially defect- free hcp Ti during an MTP MD run at 1665 K and V = 18.62 Å3/atom. (a)–(c) show MD snapshots and (d)–(f) the cor- responding relaxed structures, viewed from the [1120] direction. The arrows indicate displacement vectors of the atoms resulting in migration of the defects. Self-interstitials are marked by the red balls, while the squares indicate vacancies. occurs at temperatures below the transition temperature to the hcp phase, allowing us to construct a reasonable grid for the bcc phase. Similarly as for hcp, diffusion occurs also within the bcc structure, however at slightly lower temperatures, in partic- ular for Hf. Figure 4 illustrates the corresponding formation of Frenkel defects in bcc Ti at the experimental melting point. During the first 12 ps, one Frenkel defect is formed, followed by another during the next 3 ps. The correspond- ing self-interstitials form 〈110〉 dumbbell configurations with FIG. 4. Frenkel defect formation and kinetics in initially defect- free bcc Ti during an MTP MD run at 1943 K and V = 18.39 Å3/atom. Snapshots corresponding to MD times of (a) 12 and (b) 15 ps are shown, viewed from the [100] direction. The arrows indicate displacement vectors of the atoms resulting in migration of the defects. Self-interstitials are marked by red balls, while the squares indicate vacancies. The thermal fluctuations are reduced for visualization by averaging the atomic positions for 0.2 ps. 184107-5 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) neighboring atoms and they migrate via the interstitialcy mechanism (i.e., by pushing neighboring atoms into intersti- tial positions). Vacancy migration along the 〈111〉 direction is also observed. The migration of the atoms, and the formation of diffusion channels along the 〈111〉 direction are in agree- ment with previous works [9,49]. Unlike for hcp, the grid for the bcc free-energy calculations overlaps with the red lines that mark the onset of Frenkel defect formation and diffusion, as seen in the bottom panel in Fig. 1. This renders the free-energy calculations for bcc more challenging than for hcp. To tackle the challenge, we increase the statistics for computing the thermal averages. Specifically, in cases when there is an onset of diffusion, we repeat the thermodynamic integration with different initial random seeds until sufficient MD steps without defect formation are col- lected for the averaging. This effectively excludes the impact of the Frenkel defects on the free-energy values. To analyze whether there is a potential impact of diffusion, we have utilized the alternate temperature-integration approach [4] at these (V, T ) points and found a small effect on the free energy (see Ref. [19]). To further stabilize the free energy surface and the required derivatives, in particular along the volume axis, we extend the grid to smaller and larger volumes than what is covered by the thermal expansion. The extension to larger volumes has to be balanced against the increasing diffusion at high temperatures. With respect to the simulations at high temperatures and volumes, it is important to note that they are likely to reside in an overheated regime for the here employed PBE exchange- correlation functional. Studies for other systems have revealed a systematic underestimation of the experimental melting temperatures by the PBE functional [50,51]. Therefore it is expected that the PBE melting points of Ti, Zr, and Hf are also similarly lower than the experimental melting points. This is reflected in the highest temperature of the Ti and Hf bcc grid in Fig. 1 being lower than the experimental melting point (green dots). Further studies of the melting properties are required to analyze this aspect. C. Finite-temperature effective quasiharmonic reference Within the direct upsampling methodology, an effective quasiharmonic potential is used as the reference for the ther- modynamic integration to the MTP. The purpose of utilizing a harmonic potential as a reference is to enable the calculation of an absolute free energy, since thermodynamic integration and the subsequent perturbative upsampling only provide the free-energy differences to DFT accuracy. In general, a system can be approximated by various (quasi)harmonic references which lead to different ways of splitting the vibrational free energy. Nevertheless, the sum of F qh(V, T ) and F ah(V, T ) [in Eq. (1)] always provides the same total vibrational free energy, coming from atomic vibrations. Hence, for a particular system at a given (V, T ) point, based on the choice of the reference, the anharmonic part gets adjusted so that the total vibrational free energy remains the same. For the stable phases considered in earlier works [4,5], the effective quasiharmonic potential was fitted to low-temperature forces, specifically at 20 K. For the dynamically stable hcp phases of the present systems, we continue with this approach and fit the quasiharmonic matrix to DFT forces of MD snapshots at 20 K. However, for the bcc phase, a low-temperature quasihar- monic dynamical matrix exhibits imaginary phonon frequen- cies due to the dynamical instability and thus cannot be used as a reference for the thermodynamic integration. An apparent solution seems to be a fit of the quasiharmonic potential to forces within the stable temperature regime of the bcc phase. Such a procedure inhibits imaginary phonon frequencies in the effective dynamical matrix. However, inhibiting the imag- inary frequencies is not sufficient to produce a physically reasonable and numerically stable reference potential. In fact, irrespective of the temperature chosen for fitting the effective quasiharmonic potential within the stability regime for the bcc phase, problems occur during the thermodynamic integration to the MTP. By choosing a temperature from the lower temperature range within the stability regime—for example, the temper- atures indicated by the green crosses in the lower panel of Fig. 1—the resulting quasiharmonic potential is too soft as compared to the true, fully anharmonic interactions. Such a soft potential leads to physically incorrect short bond dis- tances during the thermodynamic integration for coupling constants λ close to zero. These short bonds bring the MTP far out of its original phase space region in which it was trained and which contains physically relevant configurations derived from DFT. As a consequence, the integration becomes numerically unstable. Choosing instead a higher temperature from the stability regime—for example, the experimental melting points indi- cated by the green dots on the bcc free-energy surfaces in Fig. 1—leads to an unphysical and strongly curved volume dependence of the corresponding quasiharmonic free energy. In consequence, the anharmonic free energy obtained in the thermodynamic integration step has a similarly (but oppo- sitely) curved volume dependence, which leads to numerical convergence issues during the parametrization of the anhar- monic free-energy surface. For the present purpose of obtaining a starting point with an analytical description of the free energies for the dynami- cally stabilized phase, we provide the following solution. We choose an effective dynamical matrix fitted to DFT forces of MD snapshots corresponding to the experimental melting temperature (green dots in the lower panel of Fig. 1) for the thermodynamic integration. This leads to a matrix with stiffer phonon frequencies as also observed experimentally [52–54]. We then switch the effective quasiharmonic and the respective anharmonic free energies to correspond to a dynamical ma- trix fitted to a lower temperature within the stability regime (the green crosses in the lower panel of Fig. 1). With this change, we ensure that the volume dependencies of both free- energy contributions [F qh(V, T ) and F ah(V, T )] can be well parametrized. There is no loss in numerical accuracy within this modified procedure. The anharmonic free energies and quasiharmonic thermo- dynamic properties shown in the following for the bcc phases thus correspond to the “intermediate” temperatures shown by the green crosses on the bcc free-energy surfaces (900 K for Ti and Zr, and 1700 K for Hf). This point should be kept in mind when comparing thermodynamic quantities with those for the 184107-6 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) hcp phase which are based on a 20-K harmonic reference (green crosses on the hcp surfaces in Fig. 1). D. Effective 0-K E-V curve and anharmonicity in the bcc phase Conventionally, the 0-K energy-volume curve E0K(V ) in Eq. (1) is calculated by routine DFT calculations for an ideal static lattice and subsequently parametrized by, for example, the Vinet equation of state [55,56] (refer to Ref. [19]). This procedure is well-justified for phases that are dynamically stable at 0 K and for which the E -V curve has a well-defined physical meaning as the starting point of the Taylor expansion of the potential energy surface [33]. For the hcp phases studied here, we take this conventional perspective of the 0-K energy- volume curve. In contrast, the 0-K static lattice E -V curve loses its mean- ing for the high-temperature dynamically stabilized phase. This is not directly evident from the curve itself, since it resembles a standard E -V curve as observed for a stable phase. The unphysical nature of the conventional calculation becomes apparent only when atoms are displaced from their ideal static lattice positions, leading to imaginary phonon fre- quencies as mentioned in the previous subsection (Sec. II C). The unphysical nature also gets reflected in the anharmonic free energy calculations of the dynamically stabilized phase, when referenced with respect to the static lattice E -V curve. The anharmonic free energy in the stable bcc regime—when extrapolated to 0 K as a function of temperature—leads to a large negative value at 0 K. Both of the above features express the unstable nature (absence) of the dynamically stabilized bcc at 0 K and low temperatures. Regardless, the standard static-lattice E0K(V ) values do not affect the final free-energy surface in the direct upsam- pling workflow, since the 0-K energies get counterbalanced inside the F ah(V, T ) term. Therefore, in principle, we could utilize the conventional 0-K energy-volume curve for direct upsampling to obtain the full free-energy surface F (V, T ) even for the bcc phase. However, the unphysical nature of this conventional definition of the 0-K energies has inconvenient practical consequences. Therefore, based on the obtained high-accuracy results within the stable temperature regime, we provide a more consistent description of the 0-K energies for dynamically stabilized phases. Specifically, we define an “effective” 0-K E -V curve as the extrapolation of the high- temperature free-energy parametrization to 0 K. Note that in this definition, we obtain the full volume dependence of the effective 0-K energies. A modified definition of the E0K(V ) behavior, consistent with the high-temperature free energies where all relevant contributions have been taken into account, is practically more beneficial. Firstly, it calibrates the anharmonic free energy to be zero at 0 K in consistency with the usual stable phases. Secondly, the definition is consistent with the CALPHAD ap- proach which traditionally extrapolates Gibbs energies from experimentally accessible temperature and composition re- gions [57]. Comparison of the 0-K enthalpy difference between the hcp and bcc phase calculated initially using the conventional static lattice 0-K E -V curve and then with the newly defined effective 0-K curve for the bcc phase reveals large differences of up to 80 meV/atom (as shown below in Fig. 13). These differences highlight the inconsistency of the conventional definition of the 0-K energies for the dynamically stabilized phase. With this in mind, we introduce the effective 0-K curve into the direct upsampling methodology instead of the conven- tional curve, and we recalculate thermodynamic and transition properties with the more consistent definition. Although the final free energies and thermodynamic properties remain the same, the values at intermediate stages of the calculation differ, and results for the bcc phase discussed in this paper correspond to the effective 0-K energies. Specifically, we reevaluate the anharmonic free energy in the bcc phase with respect to the effective 0-K E -V curve. In previous works outlining the direct upsampling method- ology [4,5], anharmonicity was defined as the free-energy contribution coming purely from vibrations beyond the effec- tive quasiharmonic reference, and calculated with respect to the conventional 0-K energy of the static lattice. As mentioned earlier, such a calculation leads to large negative values at 0 K for the dynamically stabilized phases. Hence, for the bcc phase, we instead use the effective 0-K energies to calibrate F ah(V, T ). This procedure provides a more consistent and physically meaningful comparison of anharmonicity between the dynamically stable and stabilized phases. The reader is referred to Ref. [19] for formulas related to the calibration of the anharmonic free energies. E. Vacancy contribution in the bcc phase The direct upsampling methodology can also be utilized to include the effect of vacancies on the bulk thermodynamics. To accurately capture the effect of thermal, noninteracting vacancies in the dilute limit, one needs the full DFT free- energy surface F (V, T ) [Eq. (1)] for a sufficiently large vacancy-containing supercell, in addition to the perfect bulk free-energy surface. A Legendre transformation of both surfaces gives the corresponding Gibbs energies, and the difference between them, after accounting for the relevant scaling factors, is the Gibbs energy of vacancy formation G f vac(p, T ). The contribution of the vacancies to the full ther- modynamics including their ideal configurational entropy is then exponentially dependent on G f vac(p, T ), as given by G(p, T ) = Gbulk (p, T ) − kBT exp [ −G f vac(p, T ) kBT ] , (2) where Gbulk (p, T ) is the Gibbs energy of the ideal perfect bulk without vacancies and kB is the Boltzmann constant. Obtaining the full free-energy surface F (V, T ) for the vacancy-containing supercell is computationally challenging despite the efficiency of direct upsampling. There is an onset of diffusion and severe vacancy migration at much lower temperatures than for the perfect bulk. As a consequence, it becomes extremely difficult to perform stable thermo- dynamic integration in the vacancy-containing supercell to obtain a converged free-energy surface, and eventually cal- culate G f vac(p, T ). Moreover, the error in G f vac scales with the number of atoms in the supercell (e.g., a 1 meV/atom error scales up to 127 meV/defect when utilizing the 127-atom bcc 184107-7 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) supercell with a single vacancy), adding to the complexity of the calculation. A practical solution which captures the main physical mechanism of vacancy formation was introduced in liter- ature [5]. The G f vac(T ) was approximated by a quadratic (non-Arrhenius) dependence on T at a given pressure p. Two data points were used to parametrize the formula—the 0-K vacancy formation enthalpy H f 0K, and the Gibbs energy of vacancy formation G f vac(p, T ′) at a single high temperature T ′. The method was used to calculate the vacancy contribution for different systems [5]. However, this approach is not applicable for the dynamically stabilized phase, with problems arising in the computation of both H f 0K and G f vac(p, T ′). In particular, conventional DFT calculations to obtain H f 0K are challenging due to problems arising from the instability during relaxation. Not only are the calculations challenging, but such an ap- proach where one uses the static lattice DFT values at 0 K has been rendered physically meaningless for the dynamically stabilized phase, as elaborated in the previous section. As for obtaining G f vac(p, T ′), it becomes difficult to choose a T ′ that is sufficiently high and where stable thermodynamic integration can be performed. For example, even as low as 1400 K, the thermodynamic integration becomes unstable in bcc Ti. Using a much lower temperature to parametrize the quadratic formula makes the extrapolation of G f vac to higher temperatures inaccurate. In the present work, we propose an alternate method to obtain G f vac(p, T ) also for dynamically stabilized phases. A related approach was utilized earlier to calculate high- temperature Helmholtz energies in vacancy-containing fcc Ni [58], and to corroborate free energies (calculated using conventional thermodynamic integration) in perfect bulk sys- tems that underwent diffusion at very high temperatures [4]. In the currently proposed approach, there is no longer a de- pendency on H f 0K. The method still involves calculation of G f vac(p, T ′), but the temperature T ′ can even be as low as the smallest temperature in the grid, where it is possible to perform stable thermodynamic integration for the perfect bulk and the vacancy-containing supercell. As the first step in the current approach, we calculate G f vac(p, T ′) at the smallest temperature in the grid Tmin (cf. Fig. 1). Importantly, we perform a modified thermodynamic integration, wherein, unlike conventional thermodynamic integration, we integrate the difference (between the vacancy- containing supercell and perfect bulk) in the total energy differences over λ. An Einstein solid (force constant of 8.75 eV/Å2) is used as the reference to avoid potential com- plications that could arise for an effective quasiharmonic model for the vacancy-containing supercell. The Gibbs energy of vacancy formation at T ′ = Tmin is thus given by G f vac(p, Tmin) = ∫ 1 0 �λ dλ + pV f vac, (3) where �λ = 〈 EMTP vac − EEins vac 〉 λ − x 〈 EMTP bulk − EEins bulk 〉 λ , (4) and where V f vac is the MTP volume of vacancy formation at Tmin, and λ the coupling constant between the Einstein refer- ence and MTP. Further, in Eq. (4), EMTP vac(bulk) and EEins vac(bulk) are the total energies of the MTP and the Einstein reference of the vacancy-containing supercell (perfect bulk), respectively. The factor x properly scales the energies to account for the difference in the number of atoms in the bulk and the vacancy- containing supercell. Two separate sets of MD runs at the smallest temper- ature in the grid and the corresponding MTP equilibrium volume are performed to obtain the two thermodynamic av- erages in Eq. (4) and thus �λ. The calculations are performed on medium-sized 4 × 4 × 4 (127/128-atom) supercells. A smaller 3 × 3 × 3 (53/54-atom) supercell is commensurate between the bcc and the ω phase and stabilizes the ω phase up to considerably higher temperatures [9]. The 3 × 3 × 3 super- cell should thus be avoided in order to enable stable MD at the lowest temperature in the grid. Since the Einstein reference is consistently applied to both the vacancy-containing supercell and the corresponding perfect bulk, the long-wavelength con- tributions to the phonons (the ones that are not accounted for in the 4 × 4 × 4 supercell) cancel each other. We employ the modified thermodynamic integration ap- proach in order to minimize computational errors at this stage of the calculation. The total energy difference between the Einstein reference and the MTP, calculated individually for the vacancy-containing supercell and the perfect bulk, is highly nonlinear and spans a wide range of energy values from λ = 0 (fully Einstein forces) to λ = 1 (fully MTP forces). Hence, it becomes difficult to fit a smooth function to accu- rately calculate the individual integrals over λ. However, by calculating the difference �λ [Eq. (4)], the integrand becomes much smaller with a smoother λ dependence, thus making it easier to parametrize. We fit a tangential function [4] to the �λ points to obtain the integral. An additional challenge arises during the MD run at λ values close to 1 in the vacancy-containing supercell. There are various instances of vacancy migration (despite utilizing the lowest temperature of the grid) making it cumbersome to obtain a statistically converged internal energy difference. Hence, we fit the above-described tangential curve up to a large enough λ value where there is no migration (λ = 0.97 in this work) and extrapolate it to λ = 1 to obtain a much more accurate value of the integral. Nonetheless, there is still marginal uncertainty which gets reflected in the final thermo- dynamic properties (shown below in Fig. 8). The calculated G f vac(p, Tmin) values are marked by crosses in Fig. 5, and the respective uncertainty in the value is represented by the error bars accompanying the crosses. The second step in the current approach enables a system- atic calculation of the temperature dependence of G f vac(p, T ) from the single value at Tmin up to the highest temperature in the grid. To this end, we perform an integration of the enthalpy of vacancy formation H f vac over the inverse tempera- ture according to G f vac(p, T )=G f vac(p, Tmin) T Tmin + T ∫ 1/T 1/Tmin H f vac(p, T̃ ) d ( 1 T̃ ) , (5) 184107-8 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) FIG. 5. Vacancy formation Gibbs energy for the bcc phase of Ti (blue), Zr (red), and Hf (green) calculated with thermodynamic integration at the lowest temperature in the grid (marked by the crosses), followed by temperature integration of the enthalpy up to the highest temperature in the grid. The error bars and the shaded regions represent statistical deviation during the two stages of the calculation. which derives from the Gibbs-Helmholtz equation and where H f vac(p, T ) = Hvac(p, T ) − x Hbulk (p, T ). (6) The enthalpy Hvac(bulk) of the vacancy-containing supercell (perfect bulk) is calculated from the statistical averages of the internal energy and the volume from an MD simulation using the MTP in an NpT ensemble (constant pressure at 100 kPa). The calculation is performed separately for the vacancy-containing supercell and the perfect bulk. Since there is a probability that melting occurs already at lower tem- peratures, the enthalpy is extrapolated using a second-order polynomial up to the highest temperature point in the grid. Additional details of the whole procedure can be found in Ref. [19]. The G f vac(p, T ) curves computed according to the de- scribed procedure are shown in Fig. 5, along with the statistical uncertainty represented by the shaded regions. The G f vac(p, T ) behavior is then introduced into Eq. (2) to include the effect of vacancies on the bulk thermodynamics. III. RESULTS AND DISCUSSION A. Free energies of the hcp and bcc phases The first column of Fig. 6 shows the calculated Gibbs en- ergy at ambient pressure G(p = 100 kPa, T ) for hcp and bcc Ti, Zr, and Hf. The Gibbs energy of the thermodynamically stable phase is plotted in the respective temperature interval and denoted by the strong colors. The faded colors indicate the continuation of the Gibbs energy into the metastable regime, i.e., to temperatures where the other phase is thermodynam- ically stable. The vertical solid lines indicate the calculated hcp-bcc transition temperatures (point of crossing of the full Gibbs energies), with the relevant contributions to the free energy included. The results are compared to CALPHAD data shown by blue circles, and to the experimental transition temperatures shown as vertical dotted lines. A systematic underestimation of the full Gibbs energy (blue curves) in comparison to CALPHAD is observed for both phases of all three elements. The underestimation is better visible in the insets where we plot the Gibbs energy at higher temperatures by taking the CALPHAD data as the reference (at each temperature value). An underestimation of the Gibbs energy by the PBE functional is consistent with pre- viously studied elements [4,5] except for vanadium in Ref. [5]. The PBE-predicted transition temperature is likewise under- estimated (by −81 to −163 K; cf. Table II) in comparison to the experimental value. The underestimation of the hcp-bcc transition temperature is analogous to the underestimation of the melting temperature (solid-liquid transition) predicted by the PBE functional for other systems [50,51]. In the second column of Fig. 6, we show the free energy at the PBE-predicted transition temperature F (V, T = T hcp→bcc PBE ) with the various, relevant contributions resolved individually. For all elements, the anharmonic contribution is negative in the hcp phase, and positive in the bcc phase. The magnitude of anharmonicity is highest in Hf, due to its higher transi- tion temperature, at which the plot is made. The electronic contribution (including the coupling from vibrations) is al- ways negative, and almost equal for both phases. A closer inspection reveals that the bcc phase has a marginally stronger electronic contribution (by a few meV/atom), the source of which will be analyzed and understood in Sec. III C in terms of the electronic density of states. The purely anharmonic free energy F ah is plotted as a function of volume, for different temperatures in the third column of Fig. 6. The temperatures of the effective quasi- harmonic reference for the anharmonic calculation are given alongside the plots. As discussed in Sec. II D, the anharmonic free energies plotted here are calibrated with respect to the effective 0-K energies for the bcc phase, in order to have a physically consistent definition with the hcp phase. For all the elements, it is observed that F ah decreases with increasing temperature in the hcp phase, indicating a relative softening in the structure with temperature. The effect of volume on the hcp anharmonic free energies is small. The values are similar in magnitude in all elements at similar absolute temperatures. Naturally, the anharmonic free energy in Hf reaches almost twice the magnitude (−43 meV/atom) in comparison to Ti and Zr, owing to the higher temperature regime of the hcp phase in Hf. The softening of the hcp phase with temper- ature is similar in nature to the previously studied hcp Mg system [4]. The trend is reversed in the effective F ah of the bcc phase, where a positive increase is observed with temperature. This indicates a relative stiffening of the bcc structure, analogous to the anharmonic behavior of group V bcc refractory elements (Ta and V) with temperature [5]. There is a considerable vol- ume dependence in the anharmonic contribution in bcc Zr and Hf, which is reflected in the thermodynamic properties to be discussed in Sec. III B. Even though we use an intermediate- temperature effective quasiharmonic reference, there is still a significant anharmonic contribution in the bcc phase, as high as 62 meV/atom in Hf near the highest temperature point. The strong anharmonicity and the opposite trends in hcp and bcc at the transition temperature are explained by an 184107-9 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) FIG. 6. Ab initio calculated thermodynamic potentials for hcp and bcc Ti, Zr, and Hf: Gibbs energy G(T ) at ambient pressure 100 kPa, free energy F (V ) at T hcp→bcc PBE , and anharmonic free energy F ah(V, T ) using the PBE exchange-correlation functional. The G(T ) values are referenced to the minimum energy of the hcp static lattice at 0 K. Calculations using the CALPHAD method [59] (aligned to the ab initio hcp values at room temperature) are shown in blue dots for comparison. The insets contain the full ab initio Gibbs energy at high temperatures taking the CALPHAD values as a reference at each temperature. For G(T ) and F (V ), the different excitation mechanisms are resolved, except for thermal vacancies as their contribution is negligible on this scale. In the first column, the solid and dashed vertical lines correspond to the PBE-predicted and experimental hcp-bcc transition temperature, respectively. In the second column, the dashed lines correspond to the PBE-predicted volume at 0 K and the transition temperature. In the third column, the vertical lines correspond to the equilibrium volume at the smallest and largest temperature in the grid. The temperature at which the effective quasiharmonic potential is fitted is mentioned alongside the third column. The plotted error bars (mostly hidden by the symbols) denote the 95% confidence interval. 184107-10 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) TABLE II. Hcp-bcc transition properties of Ti, Zr, and Hf from ab initio calculations in comparison with experimental values: transition temperature T hcp→bcc, enthalpy difference �Hbcc-hcp, entropy difference �Sbcc-hcp, and volume change �V bcc-hcp = V bcc − V hcp with the hcp and bcc volumes V hcp and V bcc. All differences and the absolute volumes correspond to the respective transition temperature. The ab initio errors are derived from the 95% confidence interval of the fitted anharmonic contribution to the free energy surface, which itself is derived from the estimated covariance matrix of the fitting coefficients. Ti Zr Hf GGA-PBE Experiment GGA-PBE Experiment GGA-PBE Experiment 1155a 1139a 2016a T hcp→bcc (K) 1013 ± 3 976 ± 6 1935 ± 6 1166 ± 10b 1135 ± 10b 2054 ± 50b �H bcc-hcp (kJ/mol) 3.76 ± 0.02 4.17 ± 0.13b 3.39 ± 0.04 4.02 ± 0.3b 4.72 ± 0.01 5.91 ± 0.2b (meV/atom) 38.9 ± 0.3 43.2 ± 1.3b 35.1 ± 0.4 41.6 ± 3b 48.9 ± 0.1 61.2 ± 2b �Sbcc-hcp [J/(mol K)]c 3.71 ± 0.01 3.58 ± 0.12b 3.47 ± 0.02 3.54 ± 0.27b 2.44 ± 0.02 2.88 ± 0.12b (kB/atom)c 0.446 ± 0.002 0.430 ± 0.014b 0.417 ± 0.002 0.426 ± 0.032b 0.293 ± 0.002 0.350 ± 0.014b �V bcc-hcp (Å3/atom) −0.059 ± 0.001 −0.03d −0.186 ± 0.001 −0.11d −0.008 ± 0.003 −0.2d, 0.1d Vhcp (Å3/atom) 17.780 ± 0.001 18.15e 23.842 ± 0.001 23.72g 23.577 ± 0.002 23.8i Vbcc (Å3/atom) 17.721 ± 0.001 18.13f 23.656 ± 0.001 23.61h 23.568 ± 0.002 23.62j, 23.9k aReference [59]. bReference [68]. c�Sbcc-hcp = �H bcc-hcp/T hcp→bcc. dThe experimental �V bcc-hcp values are obtained from differences of the absolute hcp and bcc volumes given in the respective rows. eThis Vhcp is calculated from the lattice constants at 28 ◦C [69] and the thermal expansion at the transition temperature [61]. fReference [70]. The same value is also obtained from the lattice constant at 28 ◦C [69] and the thermal expansion at the transition temperature [61]. gReference [71] at 1125 K. hReference [72]. iThe value is from an equation of state at 1960 K [73] which is linearly extrapolated to the experimental transition temperature. jThe value corresponds to 2073 K [74]. kFigure 5 in Ref. [66]. analysis of the first nearest-neighbor distributions, cf. Fig. 7. In the first two columns, we show the nearest-neighbor dis- tributions in bcc Ti, Zr, and Hf projected on the (110) plane, during a dynamics run at the corresponding transition temperature using the MTP and the effective quasiharmonic potential, respectively. Similar distribution plots for the hcp phase are shown in Fig. S9 in Ref. [19]. A strong anharmonic- ity is clearly visible from the MTP plots, as noticed by the wall-like shape of the distribution as the two atoms approach each other. An effective quasiharmonic potential fitted within the stable regime of the bcc phase is unable to replicate the physically correct distribution owing to the inherent symmet- ric nature of the model. This inability is reflected in the large anharmonic free energies. Such a dissimilarity in the shape of the distribution is observed both for the hcp and bcc phases, and is consistent with previous findings [5,18,60]. To pinpoint the difference in anharmonicity between bcc and hcp, we subtract the distributions shown in the first two columns in Fig. 7 for bcc and Fig. S9 for hcp (each normal- ized with respect to the density at the static nearest-neighbor distance) and plot the resulting difference in the third and fourth column of Fig. 7 for bcc and hcp, respectively. The bcc and hcp distribution-density differences are markedly distinct. Large blue patches around the nearest neighbor in bcc signify a negative density-distribution difference between the MTP, which includes anharmonicity, and the effective quasiharmonic model. This implies a smaller relative motion of the nearest-neighbor atoms and a stiffer potential (and positive free-energy contribution) when including anharmonic effects. On the other hand, patches of red signify a positive hcp density-distribution difference, implying a larger rela- tive motion of the nearest neighbor, and hence a softening (and negative free-energy contribution) when including an- harmonicity. The difference in trends seen in the anharmonic free energies in Fig. 6 (last column) can hence be directly corroborated with the first-neighbor distribution densities in the two phases. B. Thermodynamic properties of the hcp and bcc phases Figure 8 shows the entropy S, the isobaric heat capacity Cp, the coefficient of thermal expansion α, and the adiabatic bulk modulus BS , for all three investigated elements while including the different contributions to the free energy. For each quantity, the results for both the hcp and bcc phase are plotted together on the same plot, separated by the ab initio calculated transition temperature (represented by a vertical solid line). The entropy and heat capacity are compared to existing CALPHAD data represented with blue circles. The thermal expansion and bulk modulus are compared to existing experimental data represented with blue squares. Overall, we find a good agreement of the final thermo- dynamic properties (blue curves, including all contributions) with CALPHAD and experimental data. The agreement is, in fact, remarkable for Ti except for the bulk modulus of the bcc phase. The S and Cp curves are in close agreement with the CALPHAD data for all elements. Comparatively, the deviations observed in α and BS of Zr and Hf are larger. The overestimation of α and underestimation of BS can be 184107-11 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) FIG. 7. First nearest-neighbor (1NN) distributions in bcc Ti, Zr, and Hf projected on the (110) plane during an MD run at T hcp→bcc PBE (explicit values given alongside the plots) and at the corresponding equilibrium volume using the MTP (first column) and the effective harmonic potential (second column). The corresponding distributions for the hcp phase are shown in Fig. S9 in Ref. [19]. The difference in the normalized 1NN distributions during the MTP MD and the effective harmonic MD for bcc and hcp is shown in the third and fourth column, respectively. explained by the inherent limitation of the PBE exchange- correlation functional, as also observed in earlier works [4,5]. However, the disagreement to experiments is unusually high in BS for the bcc phase of Zr and Hf. For bcc Zr, we com- pare our results to two sets of experimental data—the second data point shown by a light blue square. The first set of experimental points are unnaturally decreasing with temper- ature, while the second experimental data point—which is an extrapolated value from higher pressures at a tempera- ture below the transition temperature—is wide apart from the other experimental data set. In bcc Hf, the single available experimental BS measurement is untypically higher than the corresponding values for the hcp phase, which drop with temperature. Hence, for the bulk modulus of the bcc phase, not only are there insufficiently few data points to make concrete comparisons, but the abnormalities suggest potential errors in existing experimental measurements. Based on the overall consistency of our results, the mostly good agreement with CALPHAD/experiment (also below for the hcp-bcc transi- tion properties), and the validation of the direct upsampling methodology with experiments in previous works [4,5], we can confidently state that our current ab initio results provide a better benchmark for future studies. For all three systems, for both phases, the electronic and anharmonic free-energy contributions play a crucial role in the thermodynamic properties. One exemplary observation is the considerable effect of the anharmonic free energy on α and BS in bcc Zr and Hf, which arises from the volume dependence of F ah discussed in Sec. III A (cf. third column of Fig. 6). In certain cases (for instance BS in bcc Hf), without the anhar- monic contribution, even the general trend with temperature is incorrect. Vacancies on the other hand have a noticeable effect only at the highest temperatures, in particular on Cp, due to the exponential dependence of the Gibbs energy on the vacancy formation energy [Eq. (2)]. The effect is strongest in Hf, somewhat smaller in Zr, and almost negligible in Ti. The impact of vacancies in Hf is indeed quite significant, even when compared to earlier results for bcc refractory el- ements with much higher melting points [5]. To understand this, we analyze the Gibbs energy of vacancy formation for various systems, as shown in Fig. 9. The results for Ti, Zr and Hf are shown only in the temperature range where bcc is the thermodynamically stable phase and the shaded regions represent the statistical uncertainty as discussed in Sec. II E, which subsequently propagates into the thermodynamic prop- erties. The Gibbs energies of vacancy formation of bcc Ti, Zr, and Hf are much smaller than for the bcc refractory elements, and similar in magnitude to those in fcc Al, Cu [67], and Ni [58]. Thus, the group IV elements are very susceptible to the formation of vacancies. Among the elements with a low 184107-12 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) FIG. 8. Ab initio calculated entropy S(T ), isobaric heat capacity Cp(T ), coefficient of thermal expansion α(T ), and adiabatic bulk modulus BS (T ) for hcp and bcc Ti, Zr, and Hf, in the temperature regime of their thermodynamic stability determined by the PBE transition temperature (vertical solid lines). Results considering different excitation mechanisms [effective quasiharmonic (qh), anharmonic (ah), electronic (el), and vacancies (vac; for bcc only)] are shown. The shaded region around the dashed-dotted lines marks the statistical uncertainty from the vacancy calculations. The calculation results are compared to the following CALPHAD (blue circles) or experimental data (blue squares): hcp Ti [59,61,62], bcc Ti [59,61,63], hcp Zr [59,61,64], bcc Zr [53,59,61,65], hcp Hf [59,61,64], and bcc Hf [54,59,66]. The BS value for bcc Zr labeled “Exp. 2 for bcc” is estimated from the experimental BT of Ref. [65] and our calculated data for the difference of BS − BT = 3 GPa (see Table S18 in Ref. [19]). The vertical dashed lines correspond to the respective experimental melting point. vacancy formation Gibbs energy, only Hf (and Zr to some extent) has a melting point that is high enough to really probe the effect of vacancies, which consequentially gets reflected in Cp (and also BS) of bcc Hf in Fig. 8. C. Hcp-bcc transition properties In Table II, the ab initio-calculated hcp-bcc transition tem- peratures and other transition properties are collected. The transition temperature of Hf is almost double that of Ti and Zr. This higher transition temperature arises from the larger effective 0-K energy difference between the hcp and the bcc phases, i.e., 102 meV/atom for Hf versus ≈60 meV/atom for Ti and Zr (see Table III). The higher transition temperature correlates with a higher enthalpy difference and a smaller en- tropy difference for Hf as compared to Ti and Zr (cf. Table II). These computed relations are in good agreement with exper- iment. With regards to the absolute transition temperatures, 184107-13 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) FIG. 9. Vacancy formation Gibbs energy for the bcc phase of Ti (blue), Zr (red) and Hf (green) plotted in the temperature range where bcc is the thermodynamically stable phase. The shaded regions represent the 95% confidence interval from statistical uncertainty. Vacancy formation Gibbs energies of previously studied elements (bcc refractories Mo, Ta, V, and W [5], fcc Al and Cu [67], and fcc Ni [58]) are plotted for comparison. for all three systems, the PBE exchange-correlation functional gives an underestimation as compared to the experimental val- ues, by −81 to −163 K, as already mentioned. The enthalpy differences are also consistently underestimated by the PBE functional, while there is no clear picture for the entropy and volume differences. For a more thorough comparison, it is better to analyze the transition properties as a function of temperature. To this end, Fig. 10 shows the properties on a relative temperature scale, defined as the temperature with respect to the hcp-bcc transition temperature estimated by the respective approach (PBE or experiment). There is a good agreement between the ab initio values and CALPHAD data for all elements, in particular for the differences in the thermodynamic potentials. Similarly as observed for the thermodynamic properties of each phase individually, the agreement in the thermodynamic transition properties is remarkable for Ti. For Zr, we observe a small deviation in the slope of the enthalpy �Hbcc-hcp(T ) and entropy �Sbcc-hcp(T ), while for Hf the ab initio and CAL- PHAD data are slightly shifted with respect to each other. In addition to the differences in the thermodynamic potentials, we also compare the ab initio computed and the experimental/CALPHAD temperature-dependent hcp-bcc volume difference as shown in the last row of Fig. 10. It is TABLE III. 0-K extrapolated Gibbs energy difference between the bcc and hcp phases, from ab initio (this work) and CALPHAD (solid and dashed blue curves in Fig. 13, respectively) for Ti, Zr, and Hf. Ti Zr Hf �Gbcc-hcp 0K (DFT, meV/atom) 64 58 102 �Gbcc-hcp 0K (CALPHAD, meV/atom) 72 79 133 Underestimation (meV/atom) –8 –21 –31 Relative underestimation (%) –11 –26 –23 important to note that the values observed for �V bcc-hcp are very small with respect to the absolute volumes of the phases (< 1%). This fact leads to a large experimental uncertainty as is most evident for Hf with two qualitatively different experi- mental results (+0.1 and −0.2 Å3/atom; cf. Table II). Given this experimental uncertainty, the agreement between the DFT values and experiment (the squares) is good for Ti and Zr. For these two elements, there is also a CALPHAD volume parametrization available for hcp and bcc [75] that leads to the �V bcc-hcp shown by the blue circles. For Zr, the temper- ature dependence is remarkably close to our DFT prediction, while there is a discrepancy for Ti at higher temperatures. The CALPHAD parametrization, in particular the extrapolation far above the transition temperature, suffers from the experimen- tal scatter in the accessible, thermodynamically stable region used for the fitting, since the errors are magnified strongly in the extrapolated regime. We expect our �V bcc-hcp results to be a more trustable reference, because both the hcp and the bcc phase can be well modelled at temperatures above the transition temperature. In general, as already discussed for the thermodynamic properties, our current results provide a better benchmark for comparisons in future. The agreement achieved with experiments/CALPHAD in the hcp-bcc transition properties is only possible by includ- ing the electronic and in particular anharmonic contributions to the free energy. We demonstrate this in Fig. 11 where we compare the Gibbs energy differences with the relevant contributions to the free energy resolved individually. The corresponding CALPHAD data are shown as blue circles for comparison. Using only the effective quasiharmonic approx- imation for the high-temperature free-energy contributions predicts a linear behavior of the Gibbs energy difference seen by the yellow dotted lines. A linear behavior using effective harmonic models was also observed earlier for Ti [17,76] (cf. Fig. S8 in Ref. [19]). The curvature of the CALPHAD Gibbs energy difference between the hcp and the bcc phases is captured only by including the fully anharmonic vibrational contribution as noticed by the orange dashed lines in Fig. 11. The inclusion of anharmonicity also significantly increases the transition tem- perature. The shift in transition temperature can be explained in terms of the anharmonic free energies in the third column of Fig. 6. A positive F ah for the bcc phase, and a negative contribution to the hcp phase (resulting in a difference of 26 to 88 meV/atom), shifts the hcp-bcc transition to a much higher temperature. The shift in the transition temperature for Hf is as high as 979 K and with this almost three times higher than for Ti and Zr. The stronger impact of anhar- monicity on the transition temperature in Hf arises from the relatively larger anharmonic free-energy difference between the bcc and hcp phases at the transition temperature, which in turn is a consequence of the higher transition temperature in Hf. Figure 11 also reveals the effect of the electronic con- tribution on the transition temperatures. The electronic free energies, which are obtained from the perturbative upsam- pling step, also include the coupling to atomic vibrations. There is only a gradual shift in the Gibbs energy differences caused by the electronic contribution, lowering the transition temperature almost equally in all three elements (by 83 to 184107-14 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) FIG. 10. Difference in Gibbs energy �Gbcc-hcp, enthalpy �Hbcc-hcp, entropy �Sbcc-hcp, and volume �V bcc-hcp between the bcc and hcp phases for Ti, Zr, and Hf at ambient pressure (100 kPa) from the current ab initio calculations (GGA-PBE) compared with CALPHAD data [59,75] and experiment as referenced in Table II [61,66,69–74]. The temperature is rescaled to the GGA-PBE or experimental hcp-bcc transition temperature, respectively. 172 K). The drop in transition temperature can be explained by analyzing the bcc and hcp electronic density of states (eDOS) in the three elements. Figure 12 shows the average eDOS of forty snapshots generated at the respective transition temperature for the hcp and bcc phases of each of the three elements. The eDOS is calculated using DFT while including the electronic temperature. In all elements, the value of the eDOS at the Fermi level is higher for the bcc phase. This indicates a marginally higher electronic free energy for bcc. Including the electronic contribution thus stabilizes the bcc phase, and thereby brings down the transition temperature. Compared to the anharmonic effects, however, the effect of the electronic contribution is much smaller, as reflected in the corresponding free-energy difference between the two phases. In Fig. 12, we note that the eDOS’s for both phases are more close to one another in Hf due to the increased vibrational smearing coming from the higher transition tem- perature of Hf. Despite the smaller difference in eDOS for Hf, its electronic contribution to the transition properties is stronger than for Ti and Zr, due to the higher transition temperature. The vacancy contribution to the bcc Gibbs energy at the transition temperature is negligible for Ti and Zr (well be- low 0.01 meV/atom) and small for Hf (0.4 meV/atom) as obtained from the MTP calculations (Sec. II E). The impact on the transition temperature can thus be expected to be negligible for Ti and Zr and small in Hf. To quantify this statement further studies are needed that include the effect of vacancies for the hcp phase and that go beyond the MTP level of accuracy. 184107-15 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) FIG. 11. Difference in Gibbs energy �Gbcc-hcp between the bcc and hcp phases for Ti, Zr, and Hf at ambient pressure (100 kPa) from ab initio calculations. Results considering different contributions [effective 0-K energy (0K) with effective quasiharmonic (qh), anharmonic (ah), and electronic (el)] are shown. The values are compared with CALPHAD data [59] shown as blue circles. The temperature is rescaled with the respective GGA-PBE or experimental hcp-bcc transition temperature. FIG. 12. Electronic density of states (eDOS) for hcp and bcc Ti, Zr, and Hf at the PBE transition temperature and at the corresponding equilibrium volume. In each case, the darker curve indicates an average over 40 snapshots, while the light shaded region represents the spread of the eDOS’s. The gray lines represent the Fermi-Dirac occupation function (ranging from 0 to 1) at the respective transition temperature. D. Resolving the ab initio-CALPHAD mismatch Thermodynamic databases of multicomponent systems critically rely on accurate descriptions of the unary end mem- bers in their thermodynamically stable, metastable, and even dynamically unstable regimes. Inaccuracies in the thermo- dynamics of the constituent elements will have a cascading effect on the multicomponent phase diagrams. Accurate phase stabilities of end members down to low temperatures are thus essential to guarantee reliable CALPHAD predictions, as needed by material scientists and engineers alike. Based on our high-accuracy free-energy surfaces obtained within the stable regime, we are able to reformulate the defi- nition of the 0-K energies for the dynamically stabilized bcc phase. By doing so, we address a long-standing mismatch in the comparison between ab initio and CALPHAD phase stabilities. In literature [6,31,32], comparisons were made between 0-K static-lattice ab initio results and CALPHAD values extrapolated to 0 K. Such comparisons for systems with a dynamically unstable 0-K phase gave a discrepancy in the phase stabilities, sometimes by a factor of two, with the DFT static-lattice stabilities being larger than the CALPHAD val- ues. However, this comparison is physically not meaningful, and the proper comparison is with the extrapolated, effective 0-K E -V energies of the dynamically stabilized phase. In Fig. 13, we thus extrapolate the full ab initio Gibbs energy difference between the hcp and bcc phase for all three elements down to 0 K (solid blue curves). For the extrap- olation, we use a second-order polynomial in temperature which well represents the curvature of the lattice stabilities both from our current results and CALPHAD data [59]. Since we extrapolate Gibbs energy differences, most higher-order CALPHAD coefficients of the bcc and hcp parametrizations cancel out, and the resulting data points are adequately fitted using the second-order polynomial. In fact, the second-order fitted curves fall on top of the curve obtained from the original parametrizations in Ref. [59] (see Fig. S4 in Ref. [19]). For the purpose of analysis, the extrapolated values are also compared to the unphysical 0-K static-lattice DFT values marked by orange crosses. 184107-16 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) FIG. 13. Comparison of the GGA-PBE (blue solid curves) and the CALPHAD extrapolations (blue dashed curves) from the stabilized bcc regime (blue circles, Ref. [59]), of the Gibbs energy difference �Gbcc-hcp between the bcc and hcp phase for Ti, Zr, and Hf at ambient pressure, 100 kPa. For a consistent comparison, the extrapolation for this analysis is done directly in �Gbcc-hcp(T ) using a second order polynomial in temperature for both our current results and CALPHAD data (cf. Ref. [19]). The 0-K GGA-PBE differences of the total energies of the static bcc and hcp lattices are indicated by the orange crosses. The temperature is rescaled to the GGA-PBE or experimental hcp-bcc transition temperature, respectively. For all three systems, we observe consistently that the finite temperature ab initio extrapolation provides a much smaller difference between hcp and bcc at 0 K than the static-lattice DFT result. The reduction in the energy differ- ence is indicated by the orange arrows and it amounts to −45 meV/atom (−41%) for Ti, −26 meV/atom (−31%) for Zr, and −80 meV/atom (−44%) for Hf. This result can be intuited by picturing the bcc static-lattice configuration (i.e., all atoms located on ideal bcc lattice positions) as a singular, high-energy point in the energy landscape. Thermal distortions of the atoms away from these ideal positions lead to a strong energy reduction. Finite-temperature vibrations therefore scan much lower energies as compared to the static bcc lattice. As an immediate consequence of the unphysical na- ture of the static-lattice energies, it is crucial to utilize the finite-temperature ab initio extrapolation for a consis- tent comparison with CALPHAD. For all three systems, the GGA-PBE extrapolation underestimates the CALPHAD ex- trapolation (see Fig. 13 and Table III). This underestimation can be traced back to a smaller curvature of the PBE Gibbs energy around the transition temperature as compared to CALPHAD. Ti shows the best agreement with CALPHAD with an underestimation of −8 meV/atom at 0 K. For Zr and Hf the underestimation is larger, −21 meV/atom and −31 meV/atom, respectively. These results are in line with our previous observation for the other thermodynamic and transition properties, i.e., that Ti shows the best agreement with experiment/CALPHAD, while there are larger discrep- ancies for Zr and Hf. It is important to stress that—irrespective of whether the CALPHAD extrapolation agrees better with the static-lattice 0-K energies or with the finite-temperature ab initio extrapolation—it should always be the latter that we utilize for a consistent and practically relevant comparison. More importantly, the extrapolated ab initio values provide a more physically motivated basis for future unary and multi- component phase diagram calculations, such as those under development for third-generation CALPHAD databases [77]. IV. CONCLUSIONS AND OUTLOOK Free-energy calculations of dynamically unstable phases have been long beset with severe obstacles arising from fundamental difficulties with 0 K or low-temperature approxi- mations. Even effective quasiharmonic Hamiltonians fitted to high temperature data fail to describe the correct temperature dependence of the free energy, as we have demonstrated in the present work. Worse yet, transition temperatures are off by several hundred, up to a thousand kelvin. Based on the developments put forth in the current article, the above issues have been resolved, and dynamically unsta- ble phases can now be computed as accurately as has been possible for stable phases. The approach treats the dynami- cally unstable phase primarily in the high-temperature regime, where it is stabilized by anharmonic vibrations. Thereby, the obstacles arising from the “dynamical instability” are averted and the phase is instead treated from the finite-temperature perspective as a “dynamically stabilized” phase. Such an ap- proach is possible due to the high accuracy achievable with the direct upsampling methodology in the description of the free-energy surface within the stable temperature regime. A crucial aspect in the refinement of the direct upsam- pling methodology has been the efficient determination of the stable temperature regime, in particular for the dynami- cally stabilized phase. The mainstay behind this development is high-accuracy machine-learning potentials, MTPs in this work, which enable a rapid scan of the entire volume- temperature space in search for the stability regime. The MTPs are trained on ab initio energies and forces of configura- tions generated at the experimental melting point of the target material. Phase-specific MTPs (separately fitted for the low- and high-temperature phase) and a single MTP spanning both phases show equal accuracy, outperforming classical inter- atomic potentials by an order in magnitude in the RMSEs. The stability boundaries in temperature and volume can thus be efficiently predicted using the MTPs for both conventionally stable and dynamically stabilized phases. 184107-17 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) For the dynamically stabilized phase, the lower bound of the stable temperature regime is given by the transition to a lower-symmetry phase (hexagonal ω), which is monitored by a structural descriptor. The upper bound is determined by the onset of Frenkel defect formation and diffusion. Due to the presence of low and high temperature bounds, the regime where one can perform stable thermodynamic inte- gration is found to be much narrower than for single-phase systems. Therefore a very carefully adjusted, dense volume- temperature grid for the explicit free-energy calculations is required to guarantee converged thermodynamic properties. The knowledge of the temperature bounds of the stable regime also facilitates the choice of an adequate temperature to fit a stable effective quasiharmonic reference for the dy- namically stabilized phase. A systematic strategy is employed where, at first, an effective quasiharmonic reference fitted at the melting point is used for thermodynamic integration. This is followed by a switch in the reference to an effective quasiharmonic potential fitted at an intermediate temperature which is at the lower end of the stable temperature regime of the bcc phase. Such a strategy assures stable thermodynamic integration across the entire (V, T ) grid, and simultaneously leads to free energies that can be well parametrized. After training the MTP, determining the stability window, and fitting the effective quasiharmonic reference, the tradi- tional direct upsampling methodology [4,5] to calculate the total free energy on a dense grid within the stable regime can be pursued for each phase separately. Thermodynamic properties of individual phases and transition properties can be then obtained from the free-energy surfaces. The effect of vacancies can be added on top of the perfect bulk thermo- dynamics, but it is challenging to do so for the dynamically stabilized phase. To tackle this challenge, we have proposed a novel two-step approach that involves (i) determination of the vacancy formation Gibbs energy at a single intermediate temperature followed by (ii) temperature-integration of the enthalpy at the MTP level of accuracy to extend the Gibbs energy to the highest temperature in the grid. Thus, in anal- ogy to the perfect bulk, vacancies are treated directly and efficiently in the relevant high temperature regime. Utilizing the high-accuracy free-energy surface, a new defi- nition of the 0-K E -V behavior for the dynamically stabilized phase has been established. The effective E -V curve is ob- tained by extending the high-temperature free-energy surface to 0 K, and it provides a more consistent description of the 0-K energetics for the dynamically stabilized phase. With this more meaningful description, one can disregard the conven- tional static-lattice E -V curve of the dynamically stabilized phase in any further analysis and comparison. Notably, the effective 0-K enthalpy should be utilized for a consistent comparison and integration with CALPHAD. As shown, such an approach corrects the long-standing mismatch in the com- parison between ab initio and CALPHAD 0-K enthalpies for dynamically stabilized phases. Using the redefined effective 0-K behavior, the anharmonic free energy of the dynamically stabilized phase has also been calibrated, in order to provide a definition of anharmonicity that is consistent with conven- tional stable phases. We have demonstrated the above methodological advance- ments on low-temperature hcp and high-temperature bcc phases of the prototype group IV elements Ti, Zr, and Hf. All three elements exhibit an opposite trend in anharmonic- ity between the two phases, with a negative and decreasing anharmonic free energy in the hcp phase and a positive and increasing effective anharmonic free energy in the bcc phase as temperature increases. For the dynamically stabilized bcc phase, even though a high-temperature effective quasihar- monic reference has been used, a substantial anharmonic contribution of around 50-60 meV/atom near the melting point is observed suggesting a strong asymmetric nature of vibrations that cannot be captured by harmonic models. We have illustrated this and further contrasted the anharmonic behavior of the two phases in terms of contour difference plots of the distribution of the first nearest neighbors. A comparison of the ab initio prediction to experimen- tal and CALPHAD data has revealed a good agreement for the heat capacities for all elements, and an overall strongest agreement for the thermodynamic properties of Ti. An- harmonicity, electronic excitations and their coupling to vibrations play a crucial role in the thermodynamic properties. Vacancies have an effect at high temperatures on Cp with a significant impact in Hf, and to a lesser degree in Zr. The thermal expansion coefficient and bulk modulus for Zr and Hf are systematically over- and underestimated, respectively, coming from the PBE exchange-correlation functional, as was observed earlier for other systems [4,5]. The hcp-bcc transi- tion temperatures are underestimated for all three elements (by −81 to −163 K), in line with PBE-predicted solid-liquid transitions [50]. The hcp-bcc transition properties agree very well with CALPHAD data on a relative temperature scale, once again with the strongest agreement for Ti. We have studied the effect of the different excitation mechanisms on the hcp-bcc transition properties. Including anharmonic free energies for both phases increases the hcp- bcc transition temperatures, by up to one thousand kelvin for Hf. Besides the mere shift, the anharmonic contribution is crucial in providing the correct curvature in the temperature dependence of the Gibbs energy difference between hcp and bcc. The electronic free energy (including the coupling to vibrations) has a smaller effect on the transition temperatures, bringing them down by as much as 172 K for Hf. Accurate thermodynamic properties are hence achievable only when the relevant free-energy contributions are properly included, which is possible with the direct upsampling methodology. Lastly, we have applied the newly established concept of the effective 0-K energy to the bcc phase of the three prototype elements. The effective 0-K enthalpies from DFT (GGA-PBE functional) underestimate the CALPHAD extrapolation for all three elements, by −11% to −26%. This underestimation is in strong contrast to the overestimation by up to 51% observed when utilizing the physically inconsistent bcc energies from static-lattice ab initio calculations. The physically motivated extrapolation shows the best agreement between DFT-GGA and CALPHAD for Ti, in consistency with the good agree- ment for other thermodynamic and transition properties for Ti. For Zr and Hf, a larger discrepancy of the extrapolated CALPHAD and DFT-GGA at 0 K is seen which is, how- ever, due to only a small discrepancy in the curvature of the Gibbs energy difference in the stabilized temperature window. 184107-18 DYNAMICALLY STABILIZED PHASES WITH FULL … PHYSICAL REVIEW B 108, 184107 (2023) The understanding achieved in this work on systems that feature dynamically stabilized phases paves the way for studying more complex materials in the future. For exam- ple, low-Ta-concentration TiZrHfTa high-entropy alloys were found to be dynamically unstable at low temperatures with a transformation to the ω phase [78]. Based on the procedure outlined here, such alloys can be studied up to ab initio accu- racy across the entire temperature range. The high-accuracy calculation of the free-energy surface of the dynamically stabilized phase described here can also be combined with liquid free-energy calculations in the Two-Optimized Refer- ences Thermodynamic Integration using Langevin Dynamics (TOR-TILD) scheme [50,51] to make accurate melting-point estimations. In addition to perfect bulk, the newly proposed methodology for the vacancy formation Gibbs energy also enables such calculations on vacancy-containing supercells, that will lead to more realistic melting-point predictions in the future. One of the most crucial conceptions is the ef- fective 0-K curve to describe the low temperature behavior of the dynamically stabilized phase. The redefined effective 0-K behavior is better suited for use in future CALPHAD databases for describing lattice stabilities of unaries and multicomponent alloys. Such a physically more consistent representation of dynamically stabilized phases is, in particu- lar, of high relevance to third-generation CALPHAD database development [77]. In addition, access to the full metastable temperature region, as demonstrated in the present work, offers further systematic improvement of the CALPHAD databases. In utilizing the present results, e.g., within CALPHAD, one should keep in mind that the thermodynamic proper- ties of the individual phases obtained with direct upsampling are limited by the accuracy of the DFT exchange-correlation functional (GGA in the present case). In order to moderate the discrepancy arising from the GGA functional, one can perform similar calculations also using the LDA functional. The underbinding and overbinding nature of the functionals can serve as an “ab initio confidence interval” for the re- sulting thermodynamic properties [51,79]. Alternatively, one can also resort to more advanced and accurate exchange- correlation functionals, which can be done efficiently within the upsampling step. Moreover, in an earlier work [5], we have suggested a homologous temperature scale (with respect to the corresponding melting point) for removing the dis- crepancy arising from the exchange-correlation functional on thermodynamic properties of individual phases. For transition properties and lattice stabilities computed in the present work, we have likewise performed a comparison to CALPHAD on a relative temperature scale (with respect to the correspond- ing transition temperature). Further studies are necessary to elucidate the inherent limits of DFT. In any case, the thermodynamic properties obtained in this study for systems with dynamically stabilized phases con- stitute an important contribution to the development of an ab initio-based thermodynamic database, providing a substan- tial extension beyond the previously investigated single-phase elements [4,5]. All data supporting the findings of this study are available within the paper, its Supplemental Material [19], and the DaRUS Repository [80]. The repository contains the training sets (VASP OUTCAR files), the low-MTPs and high-MTPs, the effective quasiharmonic potentials, the Gibbs energies of vacancy formation, and the final thermodynamic database (properties) for hcp, bcc Ti, Zr, and Hf. ACKNOWLEDGMENTS We appreciate fruitful discussions with Yuji Ikeda, Xi Zhang, Konstantin Gubaev, and Nikolay Zotov. We appre- ciate the help of Xiao-Gang Lu and Malin Selleby for providing an updated CALPHAD parametrization of molar volumes. This project has received funding from the Euro- pean Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 865855). The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 40/575-1 FUGG (JUSTUS 2 cluster). B.G. acknowledges the support by the Stuttgart Center for Simulation Science (SimTech). P.S. would like to thank the Alexander von Hum- boldt Foundation for their support through the Alexander von Humboldt Postdoctoral Fellowship Program. The authors gratefully acknowledge the scientific support and HPC re- sources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander- Universität Erlangen-Nürnberg (FAU) under the NHR Project No. a102cb. NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG) Grant No. 440719683. All authors designed the project, discussed the results, and wrote the manuscript. B.G. provided scripts; J.H.J., A.F., and P.S. performed the calculations. The authors declare no competing interests. [1] Q. Zheng, Y. Zhang, M. Montazerian, O. Gulbiten, J. C. Mauro, E. D. Zanotto, and Y. Yue, Understanding glass through differ- ential scanning calorimetry, Chem. Rev. 119, 7848 (2019). [2] M. van Schilfgaarde, I. A. Abrikosov, and B. Johansson, Origin of the Invar effect in iron–nickel alloys, Nature (London) 400, 46 (1999). [3] H.-Y. Chung, M. B. Weinberger, J. B. Levine, A. Kavner, J.-M. Yang, S. H. Tolbert, and R. B. Kaner, Synthesis of ultra-incompressible superhard rhenium diboride at ambient pressure, Science 316, 436 (2007). [4] J. H. Jung, P. Srinivasan, A. Forslund, and B. Grabowski, High- accuracy thermodynamic properties to the melting point from ab initio calculations aided by machine-learning potentials, npj Comput. Mater. 9, 3 (2023). [5] A. Forslund, J. H. Jung, P. Srinivasan, and B. Grabowski, Thermodynamic properties on the homologous tem- perature scale from direct upsampling: Understanding electron-vibration coupling and thermal vacancies in bcc refractory metals, Phys. Rev. B 107, 174309 (2023). 184107-19 https://doi.org/10.1021/acs.chemrev.8b00510 https://doi.org/10.1038/21848 https://doi.org/10.1126/science.1139322 https://doi.org/10.1038/s41524-022-00956-8 https://doi.org/10.1103/PhysRevB.107.174309 JUNG, FORSLUND, SRINIVASAN, AND GRABOWSKI PHYSICAL REVIEW B 108, 184107 (2023) [6] G. Grimvall, B. Magyari-Köpe, V. Ozoliņš, and K. A. Persson, Lattice instabilities in metallic elements, Rev. Mod. Phys. 84, 945 (2012). [7] A. van de Walle, Q. Hong, S. Kadkhodaei, and R. Sun, The free energy of mechanically unstable phases, Nat. Commun. 6, 7559 (2015). [8] Z.-K. Liu, Computational thermodynamics and its applications, Acta Mater. 200, 745 (2020). [9] D. Korbmacher, A. Glensk, A. I. Duff, M. W. Finnis, B. Grabowski, and J. Neugebauer, Ab initio based method to study structural phase transitions in dynamically unstable crystals, with new insights on the β to ω transformation in titanium, Phys. Rev. B 100, 104110 (2019). [10] W.-S. Ko, B. Grabowski, and J. Neugebauer, Development and application of a Ni-Ti interatomic potential with high predictive accuracy of the martensitic phase transition, Phys. Rev. B 92, 134107 (2015). [11] N. Antolin, O. D. Restrepo, and W. Windl, Fast free-energy calculations for unstable high-temperature phases, Phys. Rev. B 86, 054119 (2012). [12] Z.-G. Mei, S.-L. Shang, Y. Wang, and Z.-K. Liu, Density- functional study of the thermodynamic properties and the pressure–temperature phase diagram of Ti, Phys. Rev. B 80, 104116 (2009). [13] C.-E. Hu, Z.-Y. Zeng, L. Zhang, X.-R. Chen, L.-C. Cai, and D. Alfè, Theoretical investigation of the high pressure structure, lattice dynamics, phase transition, and thermal equation of state of titanium metal, J. Appl. Phys. 107, 093509 (2010). [14] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301(R) (2011). [15] O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Temperature dependent effective potential method for accurate free energy calculations of solids, Phys. Rev. B 87, 104111 (2013). [16] P. Liu, C. Verdi, F. Karsai, and G. Kresse, α−β phase transi- tion of zirconium predicted by on-the-fly machine-learned force field, Phys. Rev. Mater. 5, 053804 (2021). [17] S. Kadkhodaei, Q.-J. Hong, and A. van de Walle, Free energy calculation of mechanically unstable but dynamically stabilized bcc titanium, Phys. Rev. B 95, 064101 (2017). [18] P. Srinivasan, A. Shapeev, J. Neugebauer, F. Körmann, and B. Grabowski, Anharmonicity in bcc refractory elements: A detailed ab initio analysis, Phys. Rev. B 107, 014301 (2023). [19] See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevB.108.184107 for additional material includ- ing detailed computational parameters, convergence tests, figures on the absolute temperature scale, and the thermody- namic database (properties), which also includes Refs. [81–84]. [20] P. Souvatzis, O. Eriksson, and M. I. Katsnelson, Anomalous thermal expansion in α-titanium, Phys. Rev. Lett. 99, 015901 (2007). [21] Y. Nie and Y. Xie, Ab initio thermodynamics of the hcp metals Mg, Ti, and Zr, Phys. Rev. B 75, 174117 (2007). [22] P. Souvatzis and O. Eriksson, Ab initio calculations of the phonon spectra and the thermal expansion coefficients of the 4d metals, Phys. Rev. B 77, 024110 (2008). [23] C.-E. Hu, Z.-Y. Zeng, L. Zhang, X.-R. Chen, and L.-C. Cai, Density functional study of the phase diagram and ther- modynamic properties of Zr, Comput. Mater. Sci. 50, 835 (2011). [24] S. Zhang, X. Zhang, Y. Zhu, S. Zhang, L. Qi, and R. Liu, First-principles investigations on elastic and thermodynamic properties of zirconium under pressure, Comput. Mater. Sci. 61, 42 (2012). [25] H. Wang and K. Konashi, Investigation on electronic, mechani- cal and thermal properties of Hf-H system, J. Nucl. Mater. 443, 99 (2013). [26] C.-B. Zhang, W.-D. Li, P. Zhang, and B.-T. Wang, First- principles calculations of phase transition, elasticity, phonon spectra, and thermodynamic properties for hafnium, Comput. Mater. Sci. 157, 121 (2019). [27] A. V. Shapeev, E. V. Podryabinkin, K. Gubaev, F. Tasnádi, and I. A. Abrikosov, Elinvar effect in β-Ti simulated by on-the- fly trained moment tensor potential, New J. Phys. 22, 113005 (2020). [28] Y. Wang, Z. Chen, W. Windl, and J.-C. Zhao, Nonlinear Ar- rhenius behavior of self-diffusion in β-Ti and Mo, Phys. Rev. Mater. 6, 063402 (2022). [29] A. V. Shapeev, Moment tensor potentials: A class of system- atically improvable interatomic potentials, Multiscale Model. Simul. 14, 1153 (2016). [30] I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, The MLIP package: Moment tensor potentials with MPI and active learning, Mach. Learn.: Sci. Technol. 2, 025002 (2021). [31] H. L. Skriver, Crystal structure from one-electron theory, Phys. Rev. B 31, 1909 (1985). [32] Y. Wang, S. Curtarolo, C. Jiang, R. Arroyave, T. Wang, G. Ceder, L.-Q. Chen, and Z.-K. Liu, Ab initio lattice stability in comparison with CALPHAD lattice stability, Calphad 28, 79 (2004). [33] B. Grabowski, T. Hickel, and J. Neugebauer, Formation ener- gies of point defects at finite temperatures, Phys. Status Solidi B 248, 1295 (2011). [34] R. W. Zwanzig, High-temperature equation of state by a per- turbation method. I. Nonpolar gases, J. Chem. Phys. 22, 1420 (1954). [35] N. D. Mermin, Thermal properties of the inhomogeneous elec- tron gas, Phys. Rev. 137, A1441 (1965). [36] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117, 1 (1995). [37] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996). [38] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999). [39] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994). [40] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996). [41] Y. Zhou, P. Srinivasan, F. Körmann, B. Grabowski, R. Smith, P. Goddard, and A. I. Duff, Thermodynamics up to the melting point in a TaVCrW high entropy alloy: Systematic ab initio study aided by machine learning potentials, Phys. Rev. B 105, 214302 (2022). [42] R. G. Hennig, T. J. Lenosky, D. R. Trinkle, S. P. R