When the potential becomes increasingly anharmonic, however, both Eqs 12 and 13 become increasingly insufficient accounts of anharmonic free energy. Eq. 12 can only trace out real frequency shifts, as mentioned, and Eq. 13 is based on the lowest-order perturbation theory. The correction due to 3rd-order interaction ("bubble") would serve as the next-order correction to the latter, but this is more complicated to compute and still incomplete for very strongly anharmonic materials . Preparatory steps and required manual inputs are colored in yellow. Note that bulk modulus as the input can either be retrieved from the Materials Project database or obtained from DFT calculations using the previously implemented workflow in atomate. The core steps of the workflow are colored in fog blue. Decision points are colored in orange. Quantities of interest stored in the database are colored in ocean green. Other quantities (intermediate data or targeting results), which are not stored in the database, are colored in magenta. lists more detailed data on LTC. R 2 = 0.98 from the linear regression fitting here among all the materials. Figure : Performance speedup via estimated VASP runtime (in node hours) using this workflow compared with the conventional finitedisplacement method for 3rd-order anharmonic IFCs extraction. The testing structures are binaries and ternaries in seven crystal systems (cubic, hexagonal, trigonal, tetragonal, orthorhombic, monoclinic, and triclinic). The VASP runtime is the dominant portion of workflow runtime (85% of the total time) based on our tests as shown in Figure . Supercell configurations for the finite-displacement method are generated by Phono3py package (usually ≥ 1000 configurations). The estimated runtime is calculated by multiplying the configuration number and the averaged static calculations from 3-5 random configurations. More details are given in Table . : Vibrational free energy differentials (∆F) determine phase transition temperatures for Zr, ZrO 2 and GeTe (a) HCP-to-BCC phase transition in Zr, (b) tetragonal-to-cubic phase transition in ZrO 2 , (c) rhombohedral-to-cubic phase transition in GeTe. The blue lines represent free energies calculated using the quasi-harmonic approximation, denoted as ∆F(QHA), whereas the green lines correspond to free energies that include an additional anharmonic correction, indicated as ∆F(QHA+anh). The experimental phase transition temperatures (Expt. T c ) are marked in magenta squares. We separate all temperatures into three regions depending on whether the phonon modes are imaginary or not according to Figure . Temperatures are categorized into three distinct zones based on the phonon modes, as illustrated in Figure . The "imaginary zone" encompasses the temperature range where phonon modes are imaginary. Conversely, the "real zone" refers to where phonon modes are entirely real. The region between the highest imaginary-mode temperature and the lowest real-mode temperature is defined as the "transition zone" here. In this zone, the imaginary modes gradually vanish and become real. The width of this "transition zone" is determined by the interval of renormalization temperatures.