Anharmonic effects on the dynamical stability of Ce–Co–Cu intermetallic ternary compounds

Abstract

Ce-based intermetallic compounds are of growing interest for their potential applications in energy-efficient permanent magnets. While recent machine learning and DFT studies predicted several new Ce–Co–Cu ternary compounds to be stable at T = 0 K, their dynamical stability requires further investigation. We show that first-principles harmonic phonon calculations predict imaginary vibrational modes for some structures, suggesting they are dynamically unstable at 0 K. However, ab initio molecular dynamics (AIMD) simulations reveal that these structures are stable at finite temperatures. Vibrational density-of-states and phonon modes calculated at finite temperature through the AIMD simulations suggest that anharmonic interactions are important in stabilizing these predicted Ce–Co–Cu intermetallic compounds.

---

Introduction

The Ce–Co–Cu ternary system is of strong interest for exploration due to its potential to address both practical and fundamental challenges in materials design. First, cerium (Ce) is a more abundant and cheaper alternative to critical rare-earth elements like neodymium (Nd),¹ samarium (Sm), and dysprosium (Dy), which are commonly used in high-performance permanent magnets. Replacing these critical elements with Ce could improve the sustainability and economic viability of magnet technologies. Second, Ce exhibits rich and tunable electronic and magnetic behaviors due to its ability to fluctuate between Ce³⁺ and Ce⁴⁺ valence states,² making it an ideal platform for studying strongly correlated phenomena. Finally, despite the known binary phases in the Ce–Co and Ce–Cu systems,³ there is a notable absence of experimentally confirmed ordered ternary Ce–Co–Cu compounds. This lack of structural knowledge presents an opportunity for computational discovery and insight into the stability mechanisms that may govern their formation.

Discovery of new materials with emergent physical properties, such as superconductivity and exotic topological or magnetic states, has been a subject of intensive research for decades.⁴ Traditionally, most new materials were discovered through trial-and-error experiments, which require substantial time and effort. Rapid advances in modern computational hardware/software, machine learning, and advanced synthesis methods offer exciting opportunities for accelerating the pace of materials discovery. In the past few decades, several efficient new material searching methods have been developed, including random structure search,⁵ CALYPSO,⁶ genetic algorithm,^7,8 adaptive genetic algorithm (AGA),^9–15 machine-learning guided framework,^15–19 etc. These computational approaches have contributed significantly to the advances in new materials discovery. However, it should be noted that these computational methods only predicted the thermodynamic stability of the phases/structures at T = 0 K. Further validations of the dynamical stability for the predicted structures are often necessary.

Currently, most computational verification of dynamical stability of the predicted structures is based on the calculations of phonon dispersion relations using harmonic (or quasi-harmonic) approximations.^14,18 The appearance of imaginary phonon modes in such calculations is often used as a criterion to classify the structure as dynamically unstable.²⁰ However, many materials exhibit imaginary phonon modes at T = 0 K in harmonic calculations but are stable at finite temperatures. For example, group IV elements' body-centered cubic (BCC) structure, e.g., Ti, Zr, and Hf, have imaginary frequencies at 0 K but can exist at high-temperature conditions.²¹ Similarly, perovskites such as KTaO₃ ²² and CsPbBr₃ ²³ exhibit soft modes in harmonic calculations that are dynamically stabilized through strong anharmonic lattice fluctuations.

Numerous studies indicate a wide array of materials exhibiting anharmonicity. These include: (1) Perovskites and related frameworks, such as oxide and halide perovskites including KTaO₃,²² CaSiO₃,^24–26 MgSiO₃,²⁷ KNbO₃, NaNbO₃,²⁸ SrTiO₃,²⁹ CsPbBr₃,^23,30 CH₃NH₃PbI₃,³¹ and other cubic perovskites surveyed in systematic DFT studies.^32,33 These compounds exhibit soft phonon modes, quantum fluctuations, and temperature-driven instabilities. (2) Thermoelectric and chalcogenide systems like PbTe, SnTe, and SnSe,^34–38 where strong anharmonic phonon interactions suppress thermal conductivity and alter phonon dispersion. (3) Superconductors and quantum materials, including Y₂C₃,³⁹ Nb₂S,⁴⁰ H₃S,⁴¹ PdH,⁴² LuH₃,⁴³ and LaH₁₀,⁴⁴ in which anharmonic effects play a crucial role in mediating superconductivity. (4) High-pressure or high-temperature phases, such as Ti, Zr, and Hf,^21,45,46 simple-cubic Ca,⁴⁷ Be,⁴⁸ and Fe,^49,50 where lattice instabilities or quartic phonon potentials emerge under temperature or pressure changes. (5) Surfaces and low-dimensional materials, like Cu (110)⁵¹ and Be (0001),⁵² where surface phonon broadening is strongly temperature-dependent due to anharmonicity. These categories are not exhaustive but highlight the diversity of materials classes where anharmonic effects are known to play a critical role.

Several important factors beyond the harmonic interactions affect the dynamical stability of the crystalline structures.²⁰ Firstly, atoms do not sit in harmonic potential wells but in anharmonic potential wells in real materials. Because of the anharmonic interactions, the higher order expansions in the total energy cost with respect to the atomic displacement cannot be ignored, thus affecting the calculated phonon frequencies even at T = 0 K. Secondly, atoms do not sit in a potential created by static neighbors but instead by dynamic neighbors at finite temperatures. As temperature increases, atoms explore a larger part of the energy landscape through a larger vibrational amplitude, which means the atoms can access larger numbers of microstates and thus larger the system's entropy. A larger entropy with a higher temperature can significantly decrease the Gibbs free energy (G = U − TS + PV). Therefore, to assess the dynamical stability of a structure, one also needs to consider the effects of anharmonic interactions and the contribution of entropy at finite temperatures.

In this paper, we further demonstrate this point of view by studying the dynamical stability of the newly predicted Ce–Co–Cu ternary compounds¹⁵ through ab initio molecular dynamics (AIMD) simulations.^{27,28,53–57} Whereas harmonic phonon frequencies are often used as a stability filter, at the risk of discarding intrinsically anharmonic new materials, AIMD provides a more versatile approach for demonstrating the dynamic and thermal stability of newly discovered materials. By allowing large atomic displacements at finite-temperature, AIMD automatically samples the anharmonic effects and can be performed over a broad temperature range, yielding insights unavailable to harmonic approximations. AIMD is also computationally efficient for low symmetry systems: it avoids the many independent displacements required to build the full force-constant matrix in finite-displacement phonon calculations. From the resulting trajectories, we can obtain atomic probability densities; pair-correlation functions that reveal atomic neighbor relationships; and velocity autocorrelation functions (VACF), whose Fourier transforms give the vibrational density of states (DOS). Because the VACF formalism^{24,25,48,49,57–59} does not assume the underlying Hamiltonian, it captures anharmonic interactions to all orders, providing a powerful avenue to probe higher-order lattice effects beyond the reach of harmonic theory.

Methods and simulations details

The first-principles density functional theory (DFT) calculations were performed using Vienna ab initio simulation package (VASP)^60–62 with projector-augmented wave (PAW) pseudopotentials⁶³ and Perdew–Burke–Ernzerhof (PBE) exchange–correlation energy functional.⁶⁴ Converged results are obtained with a kinetic energy cutoff for the plane wave basis of 520 eV and a k-point grid with a spacing of 0.2 Å⁻¹. The equilibrium structure is optimized until all forces are below 0.01 eV Å⁻¹. DFT-PBE provides reasonable, low-cost ab initio description of structural and dynamical stability for many Cerium compounds. Since the main purpose of our present study is to compare dynamical stability from harmonic calculations versus anharmonic finite-temperature simulations, using the same DFT-PBE method for both calculations would be appropriate.

The harmonic phonon calculations are performed using the finite displacement method implemented in the Phonopy code,^65,66 based on spin-polarized forces obtained from VASP. The supercells, consisting of 128 atoms for CeCoCu₂, 200 atoms for Ce₁₀Co₁₁Cu₄, and 80 atoms for Ce₁₂Co₇Cu, are employed to capture the interatomic interactions. Due to computational constraints, these sizes may not be sufficient to achieve fully converged phonon spectra.

The ab initio molecular dynamics (AIMD) simulations were carried out using VASP without spin polarization since the compounds studied have been reported to be nonmagnetic at 0 K.¹⁵ The same supercells used in the harmonic phonon calculations are used for CeCoCu₂, Ce₁₀Co₁₁Cu₄, and Ce₁₂Co₇Cu. A lower kinetic energy cutoff of 300 eV was used to expedite the simulations. A time step of 3 femtoseconds (fs) was employed in the AIMD simulations. To validate this choice, convergence tests were performed using time steps of 1 fs, 2 fs, and 3 fs for Ce₁₀Co₁₁Cu₄ (see Fig. S3 in the SI). The resulting atomic trajectories exhibited no significant qualitative or quantitative differences across key indicators of dynamical stability, including total energy fluctuations, temperature control, atomic displacements, and the vibrational density of states (VDOS). Note that the system contains heavy elements (Ce, Co, and Cu); therefore, there are no high-frequency vibrational modes typically associated with light atoms such as hydrogen, allowing for the use of a larger time step (3 fs) in the AIMD simulations without compromising accuracy.

The supercells used in the AIMD simulations are the same as those used in the harmonic phonon calculations. For the Ce₁₂Co₇Cu, the 80-atom supercell with lattice dimensions of approximately 11.5 × 11.5 × 12.8 ų. Although this supercell contains a relatively modest number of atoms, its linear dimensions exceed those commonly employed in ab initio molecular dynamics simulations (typically around 10 × 10 × 10 ų with periodic boundary conditions). Finite-size effects on vibrational properties may exist, particularly on long wavelength dynamical behavior, but the effects on the short wavelength dynamical behavior should be minimal. Similar finite size effects are also in the harmonic phonon calculations. It is reasonable to compare the phonons from harmonic calculations and anharmonic AIMD simulations using the same supercell size and periodic boundary conditions to elucidate the anharmonic effects on the stability of the structures.

For each temperature, the system was first equilibrated for 3 ps in an NVT (constant number of particles, volume, and temperature) ensemble with temperature controlled by the Langevin thermostat. Then, the simulation is switched to an NVE (constant number of particles, volume, and energy) ensemble for more than 10 ps to collect the trajectories of the atoms for structural and dynamical analysis, including the analysis of atomic distribution functions, pair correlation functions, and velocity autocorrelation functions. The probability distribution of atomic displacement relative to the perfect structure at various temperatures provides insight into the preferred atomic positions at finite temperatures. These preferred positions are directly influenced by the shape of the potential well, as atoms tend to occupy lower-energy positions. To mimic conditions closer to experimental reality, we have also checked the results using an NPT (constant number of particles, pressure, and temperature) ensemble with temperature controlled by the Langevin thermostat.

The pair correlation function g_AB is defined as a measure of the probability of finding particles B at r away from A, relative to that for an ideal gas:

Here, N_A is the number of particles A, ρ_B is the average density of particles B, . In general, g_AB = 1 describes data with no structure (ideal gas). We then compare 3 pair correlation functions (PCF): PCF of the initial structure, PCF of the averaged structure, and time-averaged PCF. PCF of the initial structure represents the structural properties of the original structure at the initial time step. The purpose of this PCF is to provide a baseline reference, reflecting the structural order and symmetry at the beginning of the simulation. PCF of the averaged structure is calculated using the averaged atomic positions obtained from all structures at each time step. If the atom is vibrating around an equilibrium position, the positive displacement will cancel the negative displacement, and the averaged structure should be identical to the initial structure. The PCF of the average structure serves as an indicator of structural stability and anharmonic effects; if it closely resembles the initial structure, the material retains its original symmetry and order, whereas deviations suggest significant atomic rearrangements or anharmonic effects. Time-averaged PCF is calculated by averaging the PCFs of all configurations across the entire simulation trajectory. It reflects the time-averaged pair correlation over the whole simulation. Time-averaged PCF captures dynamic changes of the structure, including both equilibrium and transient states.

Velocity trajectories can be extracted from the AIMD simulation at each temperature and can generally be transformed into the phonon density of state (DOS) and phonon dispersion relation. The phonon DOS at each temperature can be calculated by taking the Fourier transform of the velocity autocorrelation function (VACF),^{24,25,48,49,57–59} defined as,

where 〈〉 is an ensemble average, and ν_n(t) is the velocity trajectory of atom n at time step t. Further projection of the velocity autocorrelation function onto each k point in the Brillouin zone, we can get the phonon spectral intensity at that temperature:where R_n is the lattice position of atom n.

Results and discussions

In our previous works,¹⁵ we predicted five stable compounds (Ce₃Co₃Cu, CeCoCu₂, Ce₁₂Co₇Cu, Ce₁₁Co₉Cu, Ce₁₀Co₁₁Cu₄) with formation energies below the convex hull, and two Co-rich low-energy compounds (Ce₄Co₃₃Cu, Ce₄Co₃₁Cu₃) with formation energies slightly above the convex hull. No imaginary vibrational frequencies were observed in the phonon dispersions of Ce₃Co₃Cu, CeCoCu₂, Ce₁₁Co₉Cu, Ce₄Co₃₃Cu, and Ce₄Co₃₁Cu₃, indicating the dynamic stability of these structures at 0 K. For Ce₁₀Co₁₁Cu₄ and Ce₁₂Co₇Cu, we obtain imaginary phonons from the harmonic calculations, which suggests that the structures might be dynamically unstable. However, harmonic calculations alone may not fully reflect the actual stability of the structure, as anharmonic interactions can play an important role in dynamical stability, especially at finite temperatures.

To see the effects of anharmonic interactions and temperature on the vibrational properties and dynamical stability, we select three representative structures, i.e., CeCoCu₂, Ce₁₀Co₁₁Cu₄, and Ce₁₂Co₇Cu, from our previous work for more detailed comparative studies. Under the harmonic approximation, CeCoCu₂ is dynamically stable without any imaginary phonon modes, Ce₁₀Co₁₁Cu₄ exhibits some imaginary modes around high symmetry k-points (Z, D, E and C₂), while Ce₁₂Co₇Cu has many imaginary modes around several k-points.

CeCoCu₂ phase

The total and projected phonon density of states (DOS) of the stable CeCoCu₂ phase obtained from harmonic calculation at T = 0 K¹⁵ and AIMD simulations at T = 500 K are shown in Fig. 1(a)–(d). The two DOS profiles are in good agreement: both exhibit phonon states extending up to 8 THz and share a prominent peak at 4 THz corresponding to the Cu atoms. The close similarity between the results demonstrates the consistency between the harmonic and AIMD-VACF approaches for stable materials. The discrepancies between the two DOS profiles (e.g., the sharp peaks in the harmonic DOS at T = 0 K are washed out in the AIMD simulations at 500 K) can be attributed to the inclusion of anharmonic effects in the VACF-based calculation, which are not accounted for within the harmonic approximation framework.

Fig. 1 Total and elemental-projected vibrational density of states (DOS) of CeCoCu₂ calculated using the harmonic approximation at 0 K are compared with those from AIMD simulations at 500 K using the velocity autocorrelation function (VACF) method. (a) Total DOS. (b–d) Partial DOS projected onto Ce, Co, and Cu atoms, respectively.

Ce₁₀Co₁₁Cu₄ phase

Imaginary phonon frequencies are observed around the Z, D, E, and C₂ points for this compound from the harmonic calculation as shown in Fig. 2(a). To examine whether the imaginary vibrational modes originate from the force convergence criterion of 0.01 eV Å⁻¹, we re-relax the structure using tighter force tolerance (≤0.001 eV Å⁻¹). The refined structure is then used for the phonon calculation. The imaginary phonon modes around the Z, D, E, and C₂ high-symmetry points persist even with the tighter force convergence criterion. However, in the projected phonon density of states (PDOS), computed using the velocity autocorrelation function (VACF) approach from AIMD simulations at 500 K, as shown in Fig. 2(b), no signs of diffusive behavior are observed. The PDOS obtained from the AIMD simulation suggests that this structure is dynamically stable at 500 K, presumably due to the anharmonic effects. The PDOS obtained from AIMD closely resembles that from the harmonic approximation, with the key difference being the absence of the imaginary phonon mode in the AIMD results. When atoms are displaced along the eigenvectors corresponding to the imaginary vibrational mode at the Z point, a shallow double-well potential is observed, with a well depth of approximately 0.003 meV per atom and minima located around 0.07 Å, as illustrated in Fig. 2(c), indicating the presence of weak anharmonicity. In this plot, the energy versus displacement is calculated by performing static total energy calculations for a series of structures generated by displacing atoms in both positive and negative directions along the mode Z eigenvector to evaluate the system's energy response to the distortion. As the anharmonic interactions are included through the AIMD, no soft modes are observed at the Z point at 500 K, as can be seen from the vibration spectra at the Z point from AIMD simulations shown in Fig. 2(d).

Fig. 2 (a) Phonon dispersion relation of Ce₁₀Co₁₁Cu₄. (b) Projected density of states (PDOS) at 500 K (NVT-NVE) obtained from AIMD. (c) The energy of the imaginary phonon mode at the Z point is shown along with its corresponding eigenvector, and the displacement magnitude reflects the largest atomic contribution within the eigenvector. (d) The corresponding phonon spectral intensity of the Z mode from the AIMD simulation at 500 K.

The dynamical stability of this structure is also confirmed by the statistics of the atomic position in the AIMD simulations, as shown in Fig. 3. As shown in the probability distribution in Fig. 3(a), atomic vibrations are centered around equilibrium positions. At lower temperatures, atomic displacements remain minimal, while increasing temperature leads to larger vibrational amplitudes due to enhanced thermal energy. To further validate structural stability, we performed an NPT simulation at 500 K. The pair correlation function (PCF) of the averaged structure (orange line) and the time-averaged PCF (green line), shown in Fig. 3(b), closely match the initial structure's PCF (blue line). These findings collectively confirm that the structure remains stable up to 800 K despite the presence of imaginary phonon modes at 0 K.

Fig. 3 (a) Probability distribution of atomic displacements for the ideal Ce₁₀Co₁₁Cu₄ at various temperatures (NVT-NVE). (b) Pair correlation function (PCF) of each element pair at 500 K (NPT). The blue line is the PCF of the initial structure, the orange line is the PCF of the averaged structure, and the green line is the time-averaged PCF of all time steps.

Ce₁₂Co₇Cu phase

Fig. 4(a) displays the phonon dispersion relation of Ce₁₂Co₇Cu from T = 0 K harmonic calculation, where imaginary phonon modes indicate the structure is dynamically unstable under the harmonic approximation. Again, the PDOS obtained from AIMD simulation at 500 K, as shown in Fig. 4(b), suggests that this structure is dynamically stable at 500 K. We note that strong imaginary modes would lead to a lower-energy structure if atoms are displaced along the corresponding eigenvector of the soft modes. To explore this, we displaced the atoms along the eigenvector at the Γ point, as illustrated in the insert of Fig. 4(c). The associated energy profile with displacement amplitude obtained from our DFT calculation is also plotted in Fig. 4(c). The resulting energy landscape reveals a double-well potential along the y-direction, with the minima located approximately 0.2 Å from the original atomic positions and a well depth of about 0.25 meV per atom. While the double well potential is responsible for the soft modes at T = 0 K, the projected density of states (PDOS) in Fig. 4(b) and the phonon spectral intensity of the X mode in Fig. 4(d) from the AIMD simulation at 500 K exhibit no evidence of diffusive behavior, indicating this structure is stable at finite temperatures. Notably, the spectral intensity of the X mode remains positive across all frequencies, indicating that this phonon mode—initially imaginary at 0 K—can be stabilized by anharmonic interactions. Furthermore, the AIMD and harmonic PDOS match well, except that the AIMD result lacks the imaginary mode.

Fig. 4 (a) Phonon dispersion relation of Ce₁₂Co₇Cu. (b) Projected density of states (PDOS) at 500 K (NVT-NVE) obtained from AIMD. (c) The energy of the imaginary phonon mode at the Γ point is shown along with its corresponding eigenvector, and the displacement magnitude reflects the largest atomic contribution within the eigenvector. (d) The corresponding phonon spectral intensity of the X mode from the AIMD simulation at 500 K.

The dynamical stability of this structure at finite temperature can also be seen from the statistical analysis of temperature-dependent atomistic trajectories from the AIMD simulations. Fig. 5(a) shows the probability distribution of atomic displacements relative to the equilibrium positions at various temperatures. At 300 K, 500 K, and 800 K, atoms are distributed symmetrically around a peak at the initial equilibrium positions, respectively. At 100 K, vibrations along the x- and z-directions also peak at the initial equilibrium position. However, along the y-direction, atoms Ce, Co, and Cu exhibit bimodal distribution, with split peaks located at 0.22 Å, 0.07 Å, and 0.45 Å, respectively, away from the initial equilibrium position, indicating occupation of two distinct positions within a double-well potential. At low temperatures, the thermal energy is insufficient to overcome the barrier between the wells, so atoms are trapped in either the left or right potential minimum.

Fig. 5 (a) Probability distribution of atomic displacements for the ideal Ce₁₂Co₇Cu at various temperatures (NVT-NVE). (b) Pair correlation function (PCF) of each element pair at 500 K (NPT). The blue line is the PCF of the initial structure, the orange line is the PCF of the averaged structure, and the green line is the time-averaged PCF of all time steps.

By examining the time-dependent atomic trajectories from the AIMD simulation at 100 K and 200 K, we observed that the Ce and Cu atoms are relaxed away from the high symmetry saddle point and randomly reside in one potential well or another with almost equal probability. Atomic hopping between the two wells is rare (quantum tunneling is not considered in AIMD simulations). Therefore, the structure below 200 K can be characterized as thermally average mixtures of locally distorted configurations. As the temperature increases, hopping events become increasingly frequent due to the very small energy barrier (about 0.25 meV per atom) between the two wells. At 300 K and above, the atomic distributions show a single peak at the high symmetry lattice points, the system thus can be considered as a dynamically stabilized single phase. Therefore, crossover from thermally average mixtures of locally distorted configurations to a dynamically stabilized single phase at a temperature between 200 and 300 K is observed for this complex compound.

This bimodal behavior at low temperature aligns with predictions from harmonic analysis, wherein modes involving y-direction displacements appear as imaginary. As the temperature rises, atoms gain sufficient kinetic energy to overcome the wells, effectively averaging the potential landscape into a single-well shape. Additional NPT simulations confirm these findings, yielding pair correlation functions (PCFs) consistent with those obtained under the NVT-NVE ensemble, as shown in Fig. 5(b). These results demonstrate that the structure remains dynamically stable, with the observed imaginary phonons at 0 K arising from the limitations of the harmonic approximation in capturing double-well potentials. Similar double-well potential and temperature-driven stabilization mechanisms have been reported in perovskites and thermoelectric,^{22,23,25,28,32,35,50,67,68} providing context for the observed behavior in Ce-based intermetallics.

Conclusions

In conclusion, we have studied the dynamical stability of three proposed Ce-based ternary compounds (CeCoCu₂, Ce₁₀Co₁₁Cu₄, and Ce₁₂Co₇Cu) using harmonic phonon calculations at 0 K and ab initio molecular dynamics (AIMD) simulations at various temperatures.

While the CeCoCu₂ phase is stable under both methods, the Ce₁₀Co₁₁Cu₄ and Ce₁₂Co₇Cu phases show imaginary modes in harmonic approximation, suggesting instability. For Ce₁₀Co₁₁Cu₄, the potential energy landscape exhibits a shallow double-well potential with a depth of only 0.003 meV per atom and minima located at ±0.07 Å, suggesting very weak anharmonicity, close to harmonic stability. For Ce₁₂Co₇Cu, a more pronounced double-well is observed with a depth of 0.25 meV per atom and minima at ±0.2 Å, indicating stronger anharmonic effects and a significant breakdown of the harmonic approximation. These observations are further supported by AIMD simulations: for Ce₁₀Co₁₁Cu₄, the atomic displacement distributions remain unimodal even at 100 K, and for Ce₁₂Co₇Cu, clear bimodal distributions emerge at both 100 K and 200 K. AIMD results confirm that all three phases are dynamically stable up to 800 K due to anharmonic effects. In both Ce₁₀Co₁₁Cu₄ and Ce₁₂Co₇Cu phase, double-well potential leads to 0 K instabilities that are mitigated at higher temperatures. Thermal fluctuations allow atoms to transition between the wells, effectively averaging the potential landscape and restoring dynamical stability.

While our discussion of anharmonic effects is qualitative, we have supplemented it with quantitative characterizations of atomic behavior that directly reflect the impact of anharmonicity. Specifically, we present atomic displacement distributions extracted from AIMD simulations, which offer a statistical description of atomic motion at finite temperatures. These distributions provide direct insight into the atomic fluctuations associated with anharmonic potential energy landscapes. In addition, we analyze real-space atomic trajectories, vibrational densities of states (VDOS), and pair correlation functions to further support the characterization of dynamical behavior. Together, these data provide a physically grounded picture of how anharmonicity contributes to the finite-temperature stability of the predicted materials.

The finite-temperature dynamical behavior of the studied compounds is evaluated using several complementary indicators derived from the AIMD simulations. First, the atomic trajectories do not exhibit long-range diffusion, suggesting that the systems remain structurally bounded over 10 ps. Second, both the instantaneous temperature and total energy fluctuate stably around their target values over simulations exceeding 10 ps, indicating numerically stable finite-temperature sampling.¹⁵ Third, the short-range structural correlations are maintained throughout the simulations, as reflected by well-defined pair correlation functions, implying the preservation of local atomic order. Finally, the vibrational spectra obtained from velocity autocorrelation functions exhibit positive spectral intensity and their close similarity to the harmonic results demonstrates the overall consistency between the harmonic and AIMD-VACF approaches. Collectively, these observations support the conclusion that the investigated structures exhibit finite-temperature dynamical stability within the timescales and conditions considered here.

These results demonstrate that imaginary phonon modes at 0 K do not necessarily indicate instability at finite temperatures. AIMD provides a more comprehensive framework for evaluating thermal stability, as it inherently captures both anharmonic effects and thermal fluctuations. Overall, our findings emphasize the importance of going beyond harmonic phonon calculations when searching for new materials to avoid overlooking promising materials that exhibit anharmonic stability at finite temperatures.

As a complementary approach, the temperature-dependent effective potential (TDEP) method^23,37,69 can be used to study anharmonic systems by fitting renormalized force constants from finite-temperature AIMD trajectories. It should be noted that AIMD becomes less reliable as the temperature approaches 0 K,⁴⁷ as its classical Newtonian mechanics framework fails to capture the quantum effects that dominate in the low-temperature limit. In these cases, the stochastic self-consistent harmonic approximation (SSCHA) is a more suitable method,^41,70 as it incorporates quantum statistics and can accurately capture anharmonic behavior beyond the harmonic approximation, even as the temperature approaches zero.

While anharmonic effects and dynamical stability are essential for understanding the physical behavior of materials, they alone do not guarantee experimental synthesizability.⁷¹ The synthesis of materials is governed by both thermodynamics and kinetic factors. Thermodynamic factors – such as having a formation energy on or below the convex hull (e.g., ΔE_hull = 0), the absence of imaginary phonon modes, and dynamic stability confirmed by molecular dynamics simulations – are commonly used to determine whether a material is energetically and dynamically stable. However, even if these thermodynamic conditions are satisfied, synthesis may still be challenging due to unknown formation pathways.⁷² Kinetic factors – such as reaction pathways,⁷³ activation barriers,⁷⁴ and experimental constraints – play a crucial role in synthesizing materials. In many cases, multiple phases may compete during synthesis, and without detailed kinetic information, it can be challenging to predict which phase will form. Most current computational studies focus on thermodynamic stability and do not explicitly mention kinetic behavior, leaving a gap in our predictive capabilities for materials discovery. In addition, materials databases are inherently biased toward successful syntheses, often omitting failed attempts, further skewing discovery efforts. A detailed treatment of kinetic obstacles is beyond the scope of this study. Nevertheless, identifying low-energy, dynamically stable structures remains crucial in materials discovery. These candidates define the landscape of potentially accessible compounds and provide valuable guidance for experimental efforts. Even without full kinetic models, exploring these stable structures allows us to uncover new materials and expand the boundaries of known chemical space.

Conflicts of interest

There are no conflicts to declare.

Data availability

Processed datasets, along with VASP INCAR, POSCAR, and XDATCAR files are publicly available in Zenodo at https://zenodo.org/records/17836343. The archive includes the following: two INCAR files used for the NVT and NVE simulations, structure files (POSCAR) for the three systems, the corresponding MD trajectories file (XDATCAR), and their displacement relative to the equilibrium positions at each time step (Displacement_###K).

Supplementary information (SI): additional analyses of the vibrational and dynamical properties of Ce–Co–Cu compounds. Fig. S1 and S2 compares harmonic and AIMD-derived vibrational density of states (DOS) for Ce₁₀Co₁₁Cu₄ and Ce₁₂Co₇Cu. Fig. S3 examines the effect of AIMD time-step size on atomic displacement distributions and DOS, while Fig. S4 evaluates the influence of spin polarization on the dynamical properties. See DOI: https://doi.org/10.1039/d5ra09680d.

References

1. A. K. Pathak, M. Khan, K. A. Gschneidner Jr, R. W. McCallum, L. Zhou, K. Sun, K. W. Dennis, C. Zhou, F. E. Pinkerton, M. J. Kramer and V. K. Pecharsky, Cerium: An Unlikely Replacement of Dysprosium in High Performance Nd–Fe–B Permanent Magnets, Adv. Mater., 2015, 27, 2663–2667.
2. R. Xie and H. Zhang, Mixed valence nature of the Ce 4f state in CeCo₅ based on spin-polarized DFT+DMFT calculations, Phys. Rev. B, 2022, 106, 224411.
3. D. Girodin, C. H. Allibert, F. Givord and R. Lemaire, Phase equilibria in the CeCo₅–CeCu₅ system and structural characterization of the Ce(Co_1-xCu_x)₅ phases, J. Less-Common Met., 1985, 110, 149–158.
4. P. C. Canfield, New materials physics, Rep. Prog. Phys., 2020, 83, 016501.
5. C. J. Pickard and R. J. Needs, Ab initio random structure search, J. Phys.: Condens. Matter, 2011, 23, 053201.
6. Y. Wang, J. Lv, L. Zhu and Y. Ma, CALYPSO: A method for crystal structure prediction, Comput. Phys. Commun., 2012, 183, 2063–2070.
7. M. P. Lourenco, J. Hostas, L. B. Herrera, P. Calaminici, A. M. Koster, A. Tchagang and D. R. Salahub, GAMaterial – A genetic-algorithm software for material design and discovery, J. Comput. Chem., 2023, 44, 814–823.
8. P. C. Jennings, S. Lysgaard, J. S. Hummelshoj, T. Vegge and T. Bligaard, Genetic algorithms for computational materials discovery accelerated by machine learning, npj Comput. Mater., 2019, 5, 46.
9. S. Q. Wu, M. Ji, C. Z. Wang, M. C. Nguyen, X. Zhao, K. Umemoto, R. M. Wentzcovitch and K. M. Ho, An adaptive genetic algorithm for crystal structure prediction, J. Phys.: Condens. Matter, 2014, 26, 035402.
10. Z. Yang, S. Wu, X. Zhao, M. C. Nguyen, S. Yu, T. Wen, L. Tang, F. Li, K. M. Ho and C. Z. Wang, Structures and magnetic properties of iron silicide from adaptive genetic algorithm and first-principles calculations, J. Appl. Phys., 2018, 124, 073901.
11. M. Sakurai, R. Wang, T. Liao, C. Zhang, H. Sun, Y. Sun, H. Wang, X. Zhao and S. Wang, et al., Discovering rare-earth-free magnetic materials through the development of a database, Phys. Rev. Mater., 2020, 4, 114408.
12. R. Wang, Y. Sun, V. Gvozdetskyi, X. Zhao, F. Zhang, L.-H. Xu, J. V. Zaikina, Z. Lin, C. Z. Wang and K. M. Ho, Theoretical search for possible Li-Ni-B crystal structures using an adaptive genetic algorithm, J. Appl. Phys., 2020, 127, 094902.
13. S. Yu, X. Zhao, S. Wu, M. C. Nguyen, Z.-Z. Zhu, C. Z. Wang and K. M. Ho, New structures of Fe₃S for rare-earth-free permanent magnets, J. Phys. D: Appl. Phys., 2018, 51, 075001.
14. C. Zhang, Y. Sun, F. Zhang, K. M. Ho and C. Z. Wang, An ultra-incompressible Mn₃N compound predicted by first-principles genetic algorithm, J. Appl. Phys., 2020, 128, 055112.
15. W. Xia, W. S. Tee, P. C. Canfield, R. Flint and C. Z. Wang, Search for stable and low-energy Ce-Co-Cu ternary compounds using machine learning, Inorg. Chem., 2025, 64, 10161–10169.
16. W. Xia, M. Sakurai, B. Balasubramanian, T. Liao, R. Wang, C. Zhang, H. Sun, K. M. Ho, J. R. Chelikowsky, D. J. Sellmyer and C. Z. Wang, Accelerating the discovery of novel magnetic materials using machine learning-guided adaptive feedback, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2204485119.
17. W. Xia, L. Tang, H. Sun, C. Zhang, K. M. Ho, G. Viswanathan, K. Kovnir and C. Z. Wang, Accelerating materials discovery using integrated deep machine learning approaches, J. Mater. Chem. A, 2023, 11, 25973.
18. Y. Zhao, M. Al-Fahdi, M. Hu, E. M. D. Siriwardane, Y. Song, A. Nasiri and J. Hu, High-throughput discovery of novel cubic crystal materials using deep generative neural networks, Adv. Sci., 2021, 8, 2100566.
19. T. Pakornchote, N. Choomphon-anomakhun, S. Arrerut, C. Atthapak, S. Khamkaeo, T. Chotibut and T. Bovornratanaraks, Diffusion probabilistic models enhance variational autoencoder for crystal structure generative modeling, Sci. Rep., 2024, 14, 1275.
20. I. Pallikara, P. Kayastha, J. M. Skelton and L. D. Whalley, The physical significance of imaginary phonon modes in crystals, Electron. Struct., 2022, 4, 033002.
21. P. Souvatzis, O. Eriksson, M. I. Katsnelson and S. P. Rudin, Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory, Phys. Rev. Lett., 2008, 100, 095901.
22. L. Ranalli, C. Verdi, L. Monacelli, G. Kresse, M. Calandra and C. Franchini, Temperature-dependent anharmonic phonons in quantum paraelectric KTaO₃ by first principles and machine-learned force fields, Adv. Quantum Technol., 2023, 6, 2200131.
23. T. Lanigan-Atkins, X. He, M. J. Krogstad, D. M. Pajerowski, D. L. Abernathy, G. N. M. N. Xu, Z. Xu, D.-Y. Chung, M. G. Kanatzidis, S. Rosenkranz, R. Osborn and O. Delaire, Two-dimensional overdamped fluctuations of the soft perovskite lattice in CsPbBr₃, Nat. Mater., 2021, 20, 977–983.
24. Z. Zhang and R. M. Wentzcovitch, Ab initio anharmonic thermodynamic properties of cubic CaSiO₃ perovskite, Phys. Rev. B, 2019, 103, 104108.
25. T. Sun, D.-B. Zhang and R. M. Wentzcovitch, Dynamic stabilization of cubic CaSiO₃ perovskite at high temperatures and pressures from ab initio molecular dynamics, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 89, 094109.
26. Y. Shin, E. D. Lucente, N. Marzari and L. Monacelli, The thermodynamics of CaSiO₃ perovskite in Earth's lower mantle, Phys. Rev. B, 2025, 112, 174113.
27. Z. Zhang and R. M. Wentzcovitch, Anharmonic thermodynamic properties and phase boundary across the postperovskite transition in MgSiO₃, Phys. Rev. B, 2022, 106, 054103.
28. K. Kim, W. Hwang, S.-H. V. Oh and A. Soon, Exploring anharmonic lattice dynamics and dielectric relations in niobate perovskites from first-principles self-consistent phonon calculations, npj Comput. Mater., 2023, 9, 154.
29. F. Libbi, A. Johansson, B. Kozinsky and L. Monacelli, Nonequilibrium quantum dynamics in SrTiO3 under impulsive THz radiation with machine learning, Sci. Adv., 2025, 11, eadw1634.
30. E. Fransson, P. Rosander, F. Eriksson, J. M. Rahm, T. Tadano and P. Erhart, Limits of the phonon quasi-particle picture at the cubic-to-tetragonal phase transition in halide perovskites, Commun. Phys., 2023, 6, 173.
31. L. D. Whalley, J. M. Skelton, J. M. Frost and A. Walsh, Phonon anharmonicity, lifetimes, and thermal transport in CH₃NH₃PbI₃ from many-body perturbation theory, Phys. Rev. B, 2016, 94, 220301(R).
32. R. X. Yang, J. M. Skelton, E. L. da Silva, J. M. Frost and A. Walsh, Assessment of dynamic structural instabilities across 24 cubic inorganic halide perovskites, J. Chem. Phys., 2020, 152, 024703.
33. J. Yang and S. Li, An atlas of room-temperature stability and vibrational anharmonicity of cubic perovskites, Mater. Horiz., 2022, 9, 1896–1910.
34. P. T. Jochym, J. Łażewski and W. Szuszkiewicz, Phonon mode potential and its contribution to anharmonism, Sci. Rep., 2020, 10, 19783.
35. J. M. Skelton, L. A. Burton, S. C. Parker, A. Walsh, C.-E. Kim, A. Soon, J. Buckeridge, A. A. Sokol, C. R. A. Catlow, A. Togo and I. Tanaka, Anharmonicity in the High-Temperature Cmcm Phase of SnSe: Soft Modes and Three-Phonon Interactions, Phys. Rev. Lett., 2016, 117, 075502.
36. G. A. Ribeiro, L. Paulatto, R. Bianco, I. Errea, F. Mauri and M. Calandra, Strong anharmonicity in the phonon spectra of PbTe and SnTe from first principles, Phys. Rev. B, 2018, 97, 014306.
37. Y. Lu, T. Sun and D.-B. Zhang, Lattice anharmonicity, phonon dispersion, and thermal conductivity of PbTe studied by the phonon-quasiparticle approach, Phys. Rev. B, 2018, 97, 174304.
38. U. Aseginolaza, R. Bianco, L. Monacelli, L. Paulatto, M. Calandra, F. Mauri, A. Bergara and I. Errea, Phonon Collapse and Second-Order Phase Transition in Thermoelectric SnSe, Phys. Rev. Lett., 2019, 122, 075901.
39. N. K. Nepal, P. C. Canfield and L.-L. Wang, Imaginary phonon modes and phonon-mediated superconductivity in compressed yttrium hydride, Phys. Rev. B, 2024, 109, 054518.
40. R. Lucrezi and C. Heil, Superconductivity and strong anharmonicity in novel Nb–S phases, J. Phys.: Condens. Matter, 2021, 33, 174001.
41. I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma and F. Mauri, Quantum hydrogen-bond symmetrization in the superconducting hydrogen sulfide system, Nature, 2016, 532, 81–84.
42. I. Errea, M. Calandra and F. Mauri, First-Principles Theory of Anharmonicity and the Inverse Isotope Effect in Superconducting Palladium-Hydride Compounds, Phys. Rev. Lett., 2013, 111, 177002.
43. R. Lucrezi, P. P. Ferreira, M. Aichhorn and C. Heil, Temperature and quantum anharmonic lattice effects on stability and superconductivity in lutetium trihydride, Nat. Commun., 2024, 15, 441.
44. I. Errea, F. Belli, L. Monacelli, A. Sanna, T. Koretsune, T. Tadano, R. Bianco, M. Calandra, R. Arita, F. Mauri and J. A. Flores-Livas, Quantum crystal structure in the 250-kelvin superconducting lanthanum hydride, Nature, 2020, 578, 66–69.
45. O. Hellman, I. A. Abrikosov and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B:Condens. Matter Mater. Phys., 2011, 84, 180301(R).
46. N. Antolin, O. D. Restrepo and W. Windl, Fast free-energy calculations for unstable high-temperature phases, Phys. Rev. B:Condens. Matter Mater. Phys., 2012, 86, 054119.
47. I. Errea, B. Rousseau and A. Bergara, Anharmonic stabilization of the high-pressure simple cubic phase of calcium, Phys. Rev. Lett., 2011, 106, 165501.
48. Y. Lu, T. Sun, P. Zhang, P. Zhang, D.-B. Zhang and R. M. Wentzcovitch, Premelting hcp to bcc transition in beryllium, Phys. Rev. Lett., 2017, 118, 145702.
49. L. Sandoval, H. M. Urbassek and P. Entel, Solid-solid phase transitions and phonon softening in an embedded-atom method model for iron, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 80, 214108.
50. Y. Ikeda, A. Seko, A. Togo and I. Tanaka, Phonon softening in paramagnetic bcc Fe and its relationship to the pressure-induced phase transition, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 90, 134106.
51. P. D. Ditlevsen, P. Stoltze and J. K. Nørskov, Anharmonicity and disorder on the Cu(110) surface, Phys. Rev. B:Condens. Matter Mater. Phys., 1991, 44, 13002.
52. J. Yang, W. Hu and S. Xiao, Anharmonic effects on Be(0001): A molecular dynamics study, Comput. Mater. Sci., 2006, 37, 607–612.
53. O. I. Malyi, K. V. Sopiha and C. Persson, Energy, phonon, and dynamic stability criteria of two-dimensional materials, ACS Appl. Mater. Interfaces, 2019, 11, 24876–24884.
54. T. Lan and Z. Zhu, Renormalized phonon microstructures at high temperatures from first-principles calculations: Methodologies and applications in studying strong anharmonic vibrations of solids, Adv. Condens. Matter Phys., 2016, 2016, 2714592.
55. C. Z. Wang, C. T. Chan and K. M. Ho, Tight-binding molecular-dynamics study of phonon anharmonic effects in silicon and diamond, Phys. Rev. B:Condens. Matter Mater. Phys., 1990, 42, 11276.
56. K. Mizushima, S. Yip and E. Kaxiras, Ideal crystal stability and pressure-induced phase transition in silicon, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 14952.
57. C. Z. Wang, C. T. Chan and K. M. Ho, Molecular-dynamics study of anharmonic effects in silicon, Phys. Rev. B:Condens. Matter Mater. Phys., 1989, 40, 5.
58. J. M. Dickey and A. Paskin, Computer simulation of the lattice dynamics of solids, Phys. Rev., 1969, 188, 5.
59. Z. Zhang, D.-B. Zhang, T. Sun and R. M. Wentzcovitch, phq: A Fortran code to compute phonon quasiparticle properties and dispersions, Comput. Phys. Commun., 2019, 243, 110–120.
60. G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15–50.
61. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
62. S. Steiner, S. Khmelevskyi, M. Marsmann and G. Kresse, Calculation of the magnetic anisotropy with projected-augmented-wave methodology and the case study of disordered Fe_1-xCo_x alloys, Phys. Rev. B, 2016, 93, 224425.
63. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
64. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
65. A. Togo, L. Chaput, T. Tadano and I. Tanaka, Implementation strategies in Phonopy and Phono3py, J. Phys.: Condens. Matter, 2023, 35, 353001.
66. A. Togo, First-principles phonon calculations with Phonopy and Phono3py, J. Phys. Soc. Jpn., 2023, 92, 012001.
67. A. N. Beecher, O. E. Semonin, J. M. Skelton, J. M. Frost, M. W. Terban, H. Zhai, A. Alatas, J. S. Owen, A. Walsh and S. J. L. Billinge, Direct Observation of Dynamic Symmetry Breaking above Room Temperature in Methylammonium Lead Iodide Perovskite, ACS Energy Lett., 2016, 1, 880–887.
68. S.-W. Kim, K. Wang, S. Chen, L. J. Conway, G. L. Pascut, I. Errea, C. J. Pickard and B. Monserrat, On the dynamical stability of copper-doped lead apatite, npj Comput. Mater., 2024, 10, 16.
69. T. Lan, C. W. Li, O. Hellman, D. S. Kim, J. A. Muñoz, H. Smith, D. L. Abernathy and B. Fultz, Phonon quarticity induced by changes in phonon-tracked hybridization during lattice expansion and its stabilization of rutile TiO₂, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 92, 054304.
70. L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea and F. Mauri, The stochastic self-consistent harmonic approximation: Calculating vibrational properties of materials with full quantum and anharmonic effects, J. Phys.: Condens. Matter, 2021, 33, 363001.
71. C. J. Bartel, Review of computational approaches to predict the thermodynamic stability of inorganic solids, J. Mater. Sci., 2022, 57, 10475–10498.
72. A. Wustrow, B. Key, P. J. Phillips, N. Sa, A. S. Lipton, R. F. Klie, J. T. Vaughey and K. R. Poeppelmeier, Synthesis and characterization of MgCr₂S₄ thiospinel as a potential magnesium cathode, Inorg. Chem., 2018, 57, 8634–8638.
73. A. J. Martinolich and J. R. Neilson, Toward Reaction-by-Design: Achieving Kinetic Control of Solid State Chemistry with Metathesis, Chem. Mater., 2017, 29, 479–489.
74. M. Aykol, J. H. Montoya and J. Hummelshoj, Rational Solid-State Synthesis Routes for Inorganic Materials, J. Am. Chem. Soc., 2021, 143, 9244–9259.

---