Approximating the impact of nuclear quantum effects on thermodynamic properties of crystalline solids by temperature remapping

When computing ﬁnite-temperature properties of materials with atomistic simulations, nuclear quantum effects are often neglected or approximated at the quasiharmonic level. The inclusion of these effects beyond this level using approaches like the path integral method is often not feasible due to their large computational effort. We discuss and evaluate the performance of a temperature-remapping approach that links the ﬁnite-temperature quantum system to its best classical surrogate via a temperature map. This map, which is constructed using the internal energies of classical and quantum harmonic oscillators, is shown to accurately capture the impact of quantum effects on thermodynamic properties at an additional cost that is negligible compared to classical molecular dynamics simulations. Results from this approach show excellent agreement with previously reported path integral Monte Carlo simulation results for diamond cubic carbon and silicon. The approach is also shown to work well for obtaining thermodynamic properties of light metals and for the prediction of the fcc to bcc phase transition in calcium.


I. INTRODUCTION
Many phenomena and properties in materials sensitively depend on the temperature at which the material is processed or operating, making it essential for modern simulation techniques to accurately account for temperature dependence. For solids, lattice vibrations beyond those predicted by the harmonic model significantly impact key material properties such as the heat capacity, thermal expansion, and free energy [1]. Well-established approximations are routinely used to account for the effect of temperature on these properties. At low temperatures, the quasiharmonic approximation (QHA) is the method of choice as it provides a computationally efficient schema to analytically derive the free energy from the easy-tocalculate force constant matrix [2,3]. A significant advantage of this approach is that it includes (harmonic) quantum effects by construction. At higher temperatures, classical molecular dynamics (MD) simulations are an established tool [4].
Using one of these two approaches to describe the temperature dependence of any material poses several limitations: (1) Anharmonic contributions, which can exist even at T = 0 K, are not captured by the QHA. (2) Classical MD simulations based on Newtonian dynamics do not include any quantum effects. (3) The asymptotic behaviors of the two approaches when going to high/low temperatures are fundamentally different. (4) Key material quantities such as the free energy and * dsouza@mpie.de entropy contain information from T = 0 K all the way up to the target temperature [see Eq. (6)]. Such quantities will be affected even at high temperatures by quantum nuclear effects. Important examples for which this issue becomes critical are the computation of accurate phase transitions [5] and the generation of thermodynamic phase diagrams, both of which are crucial for materials design.
One can address many of the above issues by splitting the atomic degrees of freedom into one or a few collective coordinates that capture the anharmonic behavior, while the remaining majority are (approximately) harmonic. In the QHA, for example, the lattice constant (volume) of the crystal captures the anharmonic behavior of a perfect bulk crystal. This separation allows one to perform for any fixed volume (the anharmonic degree of freedom) the free energy calculation over the remaining (close to harmonic) degrees of freedom in the harmonic approximation. The remaining classical anharmonic contribution (not accounted for by the QHA) in these degrees of freedom can then be computed by thermodynamic integration (TI) [6][7][8][9][10]. Since the classical anharmonic contribution vanishes at T = 0 K by construction, this approach smoothly interpolates between the quasiharmonic quantum mechanical and high-temperature classical solutions. However, this approach captures explicit anharmonicity, i.e., beyond the anharmonic degree of freedom, solely classically. It also requires a priori knowledge of the dominant anharmonic degree(s) of freedom, which is straightforward for simple bulk systems such as cubic fcc or bcc structures in which a single bulk lattice constant is often sufficient. For more complex bulk systems and systems with defects, making this decomposition is challenging and may not even be systematically possible.
A framework that overcomes all of the above-listed limitations, i.e., in which both quantum effects and anharmonicity are fully captured, was first introduced by Feynman and Hibbs in the mid-20th century [11,12]. It takes action-minimizing path integrals (PIs) over different configurations of the system and can be realized computationally with a hyperdimensional polymer approach [13] in which the quantum-classical isomorphism underlying Feynman's insights is represented by the concurrent and coupled evolution of multiple images of the system. By its nature, the approach of path integral Monte Carlo (PIMC) or path integral molecular dynamics (PIMD) is computationally demanding. Markland and Ceriotti [14] thoroughly reviewed methodological work up to the current state of the art. They focused on various acceleration approaches which allow quantum effects to be fully accounted for in systems of hundreds of atoms. These advances are significant to model systems of large sizes such as multicomponent alloys and the accompanying defects in their microstructure (e.g., dislocations, grain boundaries, phase boundaries). However, the fact that these acceleration schemes have not been used to simulate such large systems indicates that their implementations need to be more accessible and that there is a need for more direct approaches that do not require such complex implementations.
In this work, we revisit an idea for approximating nuclear quantum effects (NQEs) first introduced by Wang et al. in 1990 [15] that is easy to implement and introduces only an O(1) cost compared to running classical MD simulations. In this approach, quantum effects at a given temperature are approximated by finding a classical system at a different temperature that best represents the quantum system. While most previous applications of this approach have focused on computing conductivity and transport properties (e.g., Refs. [16][17][18][19][20][21][22][23][24][25]), Xu et al. recently used a combination of this approach with thermodynamic integration to compute the free energy of Ar [26]. In the present study, we systematically explore and analyze this concept, which we call the temperature-remapping approximation (TRA), starting with its fundamentals and discussing how it is effectively the Einstein approximation to the Debye "colored noise" approach for accelerating PIMD [27,28]. We make a rigorous comparison of the TRA with existing PI results in the literature for C and Si and use the approach to calculate material properties like the internal energy, heat capacity, thermal expansion, and free energy of metals (Al, Mg, and Ca). We thus demonstrate that the TRA captures the explicitly calculated quantum effects to within a few meV/atom at a negligible additional cost compared to running classical MD simulations.

A. Theory
Unlike classical systems, at low temperatures, quantum mechanical systems have a quantized phonon spectra and zero-point energy; that is, even at T = 0 K, the kinetic energy of the system does not drop to exactly zero. At high temperatures, quantum and classical behaviors converge, but quantum effects can contribute to the overall behavior even up to moderate temperatures for typical nuclei. This quantum/classical difference is illustrated in Fig. 1(a), which compares the probability density P of finding a single particle at a point r in a simple Lennard-Jones potential V (r) at different temperatures T , classically by Boltzmann inversion [P C ∝ exp(−V/k B T )] and quantum mechanically by solving the time-independent Schrödinger equation P QM . The classical single-particle density P C becomes increasingly sharply peaked at low temperatures, ultimately reducing to a δ spike at the minimum-energy configuration as the temperature goes down to zero. Although this one-dimensional figure is purely pedagogical, Fig. 1(b) illustrates how the probability density of the quantum system P QM can be almost perfectly matched by a classical surrogate system at a "rescaled" higher temperature. For the example potential used in Fig. 1(b), the quantum mechanical T = 0 K density is very closely matched by a classical density at 86 K.
In the case of a symmetric double-well potential where the zero-point energy is comparable to the height of the barrier like in Fig. 2(a), the quantum mechanical distribution at temperatures close to T = 0 K cannot be accurately mapped to a classical surrogate. However, the mapping improves significantly if the barrier is high enough to include more energy states in the two wells of the potential, as shown in Fig. 2(b). We can thus exploit this almost perfect match between the quantum mechanical density and the temperature-rescaled FIG. 2. Probability density P for a particle in a symmetric double-well potential V (r) (solid black line) with (a) a shallow barrier where the zero-point energy level is comparable to the height of the barrier and (b) a higher barrier at T = 0 K. The quantum mechanical density is represented by the solid blue line, and its best classical surrogate is represented by the orange dashed line with circles in both plots. The lowest energy levels are indicated by the dashed green lines. As the height of the barrier increases, more energy levels can be accommodated in the two wells of the potential, largely improving the agreement between the quantum mechanical and mapped classical densities. classical one to build a temperature map T C (T ), which provides the optimal temperature at which the classical system can best reproduce the probability density or a property in the quantum system at the target temperature T .
In principle, any temperature-dependent property can be used to construct a temperature map as long as it has both classical and quantum mechanical solutions. In addition, for the TRA to be computationally efficient, the estimation of this property must be inexpensive across the relevant range of temperatures, i.e., from T = 0 K up to the melting temperature. A property that meets these criteria is the internal energy of a system of harmonic oscillators. For an atomic system, the harmonic approximation (HA) treats all the atoms in the system as harmonic oscillators whose temperature-dependent internal energy can be obtained analytically. The internal energy for the classical system of oscillators is and the internal energy for the quantum mechanical system of oscillators is where h is the reduced Planck constant, k B is the Boltzmann constant, N is the number of atoms, and D(ν) is the phonon density of states for phonon frequencies ν obtained from the dynamical matrix. We can then replace T by T C (T ) in Eq. (1) and solve for T C (T ) by setting U QM = U C . More generally, for any quantum mechanical property X QM , the temperature map can be built using the following implicit definition for T C : This temperature map allows us to obtain for any temperaturedependent property from classical MD simulations the corresponding quantum mechanical property by simply calculating the property at T C (T ) instead of T . The choice of the internal energy of the harmonic system as the mapping property is meaningful, as it can be obtained directly as an output result of MD simulations and is critical for correctly predicting important thermodynamic properties such as the heat capacity and free energy. A limitation of building the temperature map in the harmonic limit is that quantum anharmonicity, i.e., anharmonic contributions that are probed by a quantum mechanical system even at T = 0 K, is not included. However, the final results still include anharmonic effects, as we ultimately apply this map to properties from classical MD simulations at finite temperatures that fully include the (classical) anharmonic contributions.
Among existing methods for capturing NQEs, the approach outlined above is closest to the technique of colored noise thermostats [27,28]. These thermostats reproduce the impact of quantization on nuclei by giving the thermostat target temperature an explicit frequency dependence. When combined with PIMD, these thermostats give access to a highly accurate evaluation of NQEs using fewer images than would be possible with PI alone [29]. The TRA is a radical simplification of colored noise thermostats, analogous to the relationship between the Einstein and Debye models for specific heat: Instead of the Debye-like treatment of adjusting the thermostat in a frequency-dependent way, TRA corresponds to an Einsteinlike approach that makes a single, frequency-independent (temperature-dependent) adjustment to the thermostat. Although only approximate, the TRA removes the need for any replicas, such that the only additional computational cost compared to simple, classical MD simulations is the one-time overhead of building the temperature map T C (T ), which, in this case, is simply an HA at T = 0 K.

B. Computational details
All classical MD simulations were performed using interatomic potentials within the LAMMPS simulator [30]. Phonon frequencies necessary for the HA and the QHA were obtained using the PHONOPY package [31]. Both of these were used within PYIRON [32], an IDE which facilitates reproducible and readable workflows. We chose five chemical species for our simulations: C and Si (diamond cubic), Al (fcc), Mg (hcp), and Ca (fcc with a phase transition to bcc). C and Si were chosen because PIMC studies have been conducted for these species [33][34][35] and can be used as a benchmark for the TRA. From the figures reported in these studies, we extracted numeric data using WEBPLOTDIGITIZER [36]. In order to perform the benchmarking, identical empirical potentials as used in the original PIMC studies were used: For the simulations of C, we used the same Tersoff-type potential [37] as Ref. [33], which includes a modification of parameters A and B to 1387.3 and 348.3 eV, respectively, to give better agreement of the lattice parameter with experimental values than the original parameters A and B. For Si, we used the original Stillinger and Weber potential [38] as used in Refs. [34,35]. Embedded-atom method (EAM) potentials with Finnis-Sinclair [39] parametrization were used for Al [40], Mg [41], and Ca [42].
Internal energies required to build the temperature map were obtained from the HA performed on bulk structures of cell size 4 × 4 × 4, from T = 0 K up to the maximum PIMC reference temperatures for C and Si (3000 and 1000 K, respectively), from T = 0 K up to the melting temperatures for Al and Ca (926 and 997 K, respectively, obtained from their potentials using the methodology described in Ref. [43]) and from T = 0 K up to 600 K for Mg. The Debye temperatures T Debye were also obtained for the species from the HA at T = 0 K and are shown in Table I. The temperature map was obtained following Eqs. (1) through (3).
The cell sizes used for the MD simulation are shown in Table I. They were obtained after a careful convergence test, which involved comparing the internal energies and lattice parameters simulated using increasing cell sizes (2 × 2 × 2 to 7 × 7 × 7). The optimum cell size was the size for which a further increase in cell size did not change the internal energy by more than 0.5 meV and the lattice parameter by more than 0.5 mÅ. This was followed by the same test using the optimum cell size with decreasing MD time steps (1 to 0.1 fs). The converged time step values are also shown in Table I. MD data were sampled using NPT simulations for 100 temperatures over the temperature range with five independent runs for each temperature. Each single MD simulation was run for a total simulation time of 1 ns, including 50 ps of thermal equilibration time for the optimum cell size and time step for each of the species, for a pressure of 1 atm for C and Si (as used in the PIMC studies) and zero pressure for Al, Mg, and Ca. FIG. 3. Accuracy of the polynomial fit to describe the temperature dependence of the internal energy for varying polynomial degrees. The polynomial degree is given in the legend. The example shown here has been performed for Al. For polynomials of degree 5 and above, only high-frequency noise is observed, implying that a fifth-order polynomial is sufficient.
Since the internal energy and lattice parameter vary smoothly as a function of temperature, we used a polynomial to fit these data. The associated fit error was analyzed and is shown for the example of Al in Fig. 3 for a variety of polynomials. While for the second-, third-, and fourth-order polynomials low-frequency information is evident, for polynomials of fifth order and higher only high-frequency noise is visible. This behavior was observed for all the species studied here. We therefore used a fifth-order polynomial, which guarantees fitting errors that are well below 1 meV/atom. Having this polynomial fit allows us to perform a continuous remapping of the classical results in the TRA.
For the QHA calculations with cubic symmetry, 15 volumes around the T = 0 K equilibrium volume were used. For hcp Mg, we used a 7 × 7 grid of a and c lattice parameter values.
TI simulations were performed for Al, Mg, and Ca using Frenkel-Ladd nonequilibrium thermodynamic integration. The calculations were performed according to the methodology in Ref. [44] and provided the classical Helmholtz free energy with the center of mass correction. The forward and backward sampling steps were scaled up to 100 000 steps, following a 20 000-step thermal equilibration. The resulting Helmholtz free energy was converted to Gibbs free energy following the volume-distribution approach described in Ref. [45]. The TI simulations were performed for 100 temperatures over the temperature range with five independent runs for each temperature.

A. Temperature map
We first look at the classical to quantum temperature maps obtained from the TRA for all five species. Figure 4(a) shows the calculated temperature maps, and Fig. 4(b) shows their respective ratios of the difference between the remapped and For C (solid blue line), we find that the lowest remapped classical temperature, i.e., the one describing the quantum mechanical T = 0 K temperature, is 823 K. This value is surprisingly large: To approximate zero-point vibrations, a classical simulation well above room temperature is required. As may be expected, this value is well below the Debye temperature (T Debye = 2313 K) and the sublimation temperature of diamond cubic C. It is interesting to note that even at temperatures as high as 3000 K (i.e., close to sublimation), the behavior is not purely classical: The remapped temperature still sits about 89 K above the true temperature. It indicates that even at the sublimation/melting temperature, NQEs cannot be neglected and may have to be included to obtain accurate melting temperatures. This is supported by the temperature maps of Si and Al, which also show a similar offset at the highest simulated temperature (30 K offset for Si at 1000 K, 32 K offset for Al at the melting temperature). Figure 4(a) also shows that for Si, Al, Mg, and Ca, the lowest classical temperature corresponding to T = 0 K is significantly lower (274, 138, 145, and 75 K, respectively). For the metallic systems, these values are well below their melting temperatures. On a relative temperature scale with respect to the Debye temperature, the temperature maps look much more similar, as shown in Fig. 4(b).

B. Diamond cubic C
To test the performance of the TRA for approximating NQEs, we compare structural and thermodynamic properties obtained from the TRA to those computed using PIMC simulations by Herrero and Ramírez for bulk diamond cubic C [33]. We first consider the internal energy. Next to the resulting internal energy from the TRA, Fig. 5(a) also shows the internal energies from the MD simulations and QHA, as well as PIMC literature values. By construction, the TRA results show excellent agreement with the QHA at low temperatures. With increasing temperature, the TRA internal energy smoothly approaches the classical one from the MD simulations, showing that the approach can smoothly and accurately connect the low-and high-temperature limits.
To evaluate how the TRA compares to PIMC, we look directly at the energy differences between the two approaches, which are on the meV scale, in Fig. 5(b). On this scale, the impact of size convergence becomes visible. Using the same computational cell as Ref. [33] (2 × 2 × 2), the TRA overestimates the full PIMC-calculated internal energy by about 7 meV at T = 0 K, increasing to about 12 meV at 3000 K. A significant part (about 4 meV) of the difference between the PIMC and TRA results occurs already at T = 0 K, which is visible from the difference between the QHA and PIMC internal energies. This implies that zero-point vibrations already experience anharmonicity that the HA does not capture. This inability to account for T = 0 K anharmonicity extends to the TRA since we use HA as a reference for the temperature mapping.
The computationally inexpensive nature of the TRA allows us to systematically test size convergence. As shown in Fig. 5(b), the internal energy for the 2 × 2 × 2 cell is not well converged with respect to size. Taking the 4 × 4 × 4 cell (which is size converged with respect to 5 × 5 × 5 within 1 meV/atom), we see that the absolute value and the temperature dependence change, resulting in deviations of more than 6 meV. Thus, at high temperatures, the size convergence error of the 2 × 2 × 2 cell is comparable to the error induced by the TRA. Overall, at meV accuracy, the errors associated with the TRA must be weighed against the need for size-converged cells, which are inexpensive and easy to obtain with the TRA compared to more expensive PIMC calculations. Next, we turn to other physical properties to which the temperature map can be applied. The temperature dependence of the lattice parameter is shown in Fig. 6(a). We obtain excellent agreement between the PIMC lattice parameters and those obtained at the remapped temperatures. As with the internal energy, Fig. 6(a) highlights the ability of the TRA to smoothly connect the quantum-dominated regime to the classically dominated regime in a single framework, resulting, like PIMC, in a lattice parameter between the purely classical (high temperature) and quasiharmonic (low temperature) representations.
We next calculate the isobaric heat capacity, a derived thermodynamic property which we estimate from the remapped enthalpy H of the system, using the relation where H = U + PV . This quantity can also be calculated directly from the enthalpy fluctuations of the system at the remapped temperatures using the equation The derivation of this equation is included in the Appendix. The remapped heat capacities, along with those calculated from classical MD, the QHA, Eq. (5), and the PIMC literature values, are shown in Fig. 6(b). The size-converged 4 × 4 × 4 cell shows excellent agreement with low-temperature quasiharmonic and PIMC heat capacities and lies slightly closer to the classical MD heat capacity at high temperatures than the PIMC result. We thus conclude that the TRA allows us to approximate NQEs for large (e.g., size converged) simulation ensembles and that it works very well for obtaining derived thermodynamic quantities like the heat capacity with nuclear quantum effects directly from MD simulations.

C. Diamond cubic Si
We next look at bulk diamond cubic Si and compare results from the TRA to those computed by Ramírez and Herrero [34] and Noya et al. [35] using PIMC simulations.
The internal energies from the QHA, MD, and the TRA are shown in Fig. 5(c). The TRA, PIMC, and the QHA agree almost perfectly on this scale. Figure 5(d) shows the differences with respect to PIMC as well as the cell size dependence of the internal energy. As can be seen, when using the same cell size (2 × 2 × 2 cell), the TRA internal energy agrees with the PIMC internal energy, showing only tiny deviations on the order of 1 meV. These deviations are substantially smaller than the already small errors observed for C. The reason why C shows larger deviations is the significantly higher classical temperature (823 K) compared to Si (274 K) to optimally describe zero-point vibrations. The higher this temperature is, the more the system will probe the anharmonic contributions of the potential energy surface. Thus, for C, even at T = 0 K, the zero-point energy contains anharmonic contributions that are not captured by the harmonic approximation, which is used as a reference for the TRA.
Like for C, we see that the errors due to size convergence are comparable in magnitude to those introduced by the TRA and that size convergence errors become larger at higher temperatures. This reinforces our conclusion that accuracy must be considered holistically, including not only inherent approximations but also the constraints that methodologies may impose on parameters like system size.
We also study the thermal lattice expansion of Si. While the TRA underestimates the temperature-dependent lattice parameters slightly when compared to the PIMC result at low temperatures, we obtain good convergence with the hightemperature MD lattice parameters, as can be seen in Fig. 6(c). The slight remaining deviation of the PIMC lattice parameters from the MD result at high temperatures is likely related to the fact that the Monte Carlo simulations do not accurately capture system size variations. This expectation is based on the agreement of our MD lattice parameters with MD literature values [46] for the same potential [see blue squares in Fig. 6(c)]. Finally, as with C, we obtain the isobaric heat capacity as a function of temperature using Eqs. (4) and (5). From Fig. 6(d), we see excellent agreement between the TRA and PIMC heat capacities. As before, the TRA accurately follows the trend in PIMC to sit between QHA and purely classical results.

D. Metals: Al (fcc) and Mg (hcp)
Having shown that the TRA gives results that agree well with those from PIMC studies, we use the TRA to calculate properties of bulk metals, in particular Al and Mg, for which we are not aware of any PI studies in the literature. The results are shown in Fig. 7.
For Al [ Fig. 7(a)], the TRA lattice parameter agrees with the QHA to within milliangstroms, even though the temperature map is constructed purely in the harmonic limit and no explicit information regarding volume expansion is included. For Mg [ Fig. 7(f)], the agreement is similar, and we again see the trend of a smooth connection between the quantum and classical regimes. Unlike the other materials studied here, Mg has the added complication of two independent lattice constants. Despite this additional degree of freedom in the atomic structure, the second lattice parameter, indicated by the c/a ratio, follows the same qualitative quantum-to-classical trend as the lattice parameters from the cubic systems. Figures 7(b) and 7(c) show that for Al, the desired behavior is obtained for the internal energy and the isobaric heat capacity: In both cases, the TRA accurately and smoothly connects the low-temperature quantum region with the hightemperature classical region. We further use the isobaric heat capacity to compute the Gibbs free energy using the relation G = H + T S, where H = U (because we consider a zeropressure system) and the entropy S is given by Additionally, we also add the center of mass correction described in Ref. [45]. The TRA Gibbs free energy, together with the ones computed using the QHA and fully classical TI, is shown in Fig. 7(d). We observe that the TRA Gibbs free energy curve smoothly transitions from agreement with the QHA free energy to agreement with free energy from TI. Deviations from the corresponding low-and high-temperature references are within about 7 meV/atom, as can be seen in Fig. 7(e). For Mg, very similar trends are observed for the internal energy [ Fig. 7(g)], isobaric heat capacity [ Fig. 7(h)], and the Gibbs free energy [Figs. 7(i) and 7(j)].

E. Phase transitions: fcc to bcc Ca
Finally, we consider Ca, which shows a prominent structural phase transformation from fcc at low temperatures to bcc at high temperatures. A set of results similar to those obtained for Al and Mg is shown in Fig. 8 for the fcc and bcc phases. Qualitatively, the same trends as observed for the other species are found. For the potential that we use [42], classical TI simulations predict a phase transition from the fcc to the bcc phase at ≈347 K, as can be seen in Fig. 9(a).
Note that the exact value of the transition temperature is highly sensitive to even minute changes in the Gibbs free energies of the two phases: A difference of 1 meV between the Gibbs free energies of the two competing phases shifts the transition temperature by about 30 K. Because of this sensitivity, even small changes in the relative free energies arising from NQEs may play a significant role. Our calculations show excellent agreement with the QHA difference in free energies at very low temperatures and, like before, gradually become more classical with increasing temperature. With the TRA, we predict a phase transition temperature of ≈311 K, slightly below the purely classical result of ≈347 K but significantly above the QHA prediction of ≈221 K. Due to the nature of the empirical potential used, none of the transition temperatures obtained from the QHA, TI, and the TRA accurately reflect the experimental transition temperature (716 K) or  [5]. However, this example highlights the importance of including NQEs even at elevated temperatures and demonstrates how these effects can be inexpensively estimated and included in such temperature-sensitive properties using the TRA.

IV. CONCLUSIONS
In the present study, we have systematically explored the suitability and performance of the TRA to accurately describe NQEs and anharmonicity from T = 0 K over the entire tem- perature range. An easy-to-obtain temperature map was used which was constructed in the harmonic limit from the relation between the internal energy of quantum oscillators and classical oscillators that best represent those quantum oscillators. This map was then applied to material properties obtained from classical MD simulations to obtain the same properties but with NQEs at negligible extra computational effort. The TRA accurately reproduced the low-temperature quantum harmonic and high-temperature classical anharmonic limits. Systematic comparisons with available PIMC data for C and Si, which can be regarded as a gold standard, demonstrated that the TRA can describe the temperature dependence of key material properties over the entire relevant temperature range. We then applied the TRA to various bulk metals for which PIMC literature data are unavailable. We found that NQEs have an effect even at higher temperatures and should be included when computing temperature-dependent material properties like heat capacity and free energy, which in turn are used to predict other important properties such as phase transitions.
The fundamental concept of temperature remapping can be straightforwardly extended to chemically more complex systems relevant in materials science such as binary, ternary, and multicomponent alloys. A prominent use case would be the calculation of Gibbs free energies of these material systems, needed as an input for CALPHAD (CALculation of PHAse Diagrams)-based approaches to construct bulk phase diagrams.
A scenario in which a computationally inexpensive approach like the TRA can break is when the system under consideration has multiple inequivalent bonds with very different strengths, molecular crystals, and crystals containing light nuclei or other defects. In these systems, in which two or more principal frequencies can exist, an approach based on a single Einstein mode reaches its limits, and a careful analysis of the vibrational response becomes mandatory. Extensions in the spirit of a colored noise thermostat, which allows control of individual frequency modes, would be an interesting option to explore, with applications for bulk metal or alloy structures with defects, whose quantum approximated thermodynamic properties can be used to construct defect phase diagrams [47].

APPENDIX: DERIVATION OF ISOBARIC HEAT CAPACITY FROM ENTHALPY FLUCTUATIONS
The energy (enthalpy) of a canonical ensemble at a temperature T fluctuates randomly about a fixed mean value H that can be expressed as where Z C is the canonical partition function and β = 1/k B T . The isobaric heat capacity C P at the remapped temperature T C (T ) can be obtained by substituting Eq. (A1) in Eq. (4) and using the chain rule: