Large-scale molecular dynamics simulations of negative thermal expansion in Bi(1-x)La(x)CoO(3) using a pretrained machine learning force field.

Large-scale molecular dynamics simulations of negative thermal expansion in Bi1−xLaxCoO3 using a pretrained machine learning force field - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Search Log in Primary site navigation Search Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Sci Rep . 2026 Apr 10;16:16919. doi: 10.1038/s41598-026-47983-9 Large-scale molecular dynamics simulations of negative thermal expansion in Bi 1− x La x CoO 3 using a pretrained machine learning force field Shogo Wakazaki 1, ✉ , Yuki Sakai 1, 2, 4 , Kazuki Takahashi 1 , Hena Das 1, 2, 3 , Masaki Azuma 1, 2, 3, ✉ Author information Article notes Copyright and License information 1 Materials and Structures Laboratory, Institute of Integrated Research, Institute of Science Tokyo, Yokohama, 226–8501 Kanagawa Japan 2 Kanagawa Institute of Industrial Science and Technology, Ebina, 243-0435 Japan 3 Research Center for Autonomous Systems Materialogy (ASMat), Institute of Integrated Research, Institute of Science Tokyo, Yokohama, 226- 8501 Kanagawa Japan 4 Neutron Science and Technology Center, Comprehensive Research Organization for Science and Society (CROSS), Tokai, Naka, 319- 1106 Ibaraki Japan ✉ Corresponding author. Received 2025 Dec 22; Accepted 2026 Apr 6; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13230641  PMID: 41963474 Abstract Negative Thermal Expansion (NTE) is a key property for achieving dimensional stability in advanced technologies, yet understanding its underlying mechanisms, particularly the dynamic processes of phase transitions, remains a formidable challenge. In this study, we applied large-scale molecular dynamics simulations using the Crystal Hamilton Graph Network (CHGNet), a pretrained universal machine learning force field (MLFF), to elucidate the atomic-level mechanisms of the NTE phenomenon in the perovskite oxide system Bi 1− x La x CoO 3 . CHGNet achieved a computational speed approximately 20,000 times faster than first-principles calculations and qualitatively reproduced the chemical trend of the NTE transition temperature systematically decreasing as La substitution increased. Leveraging this computational efficiency, we successfully visualized the dynamic process of the NTE phase transition at the atomic level, which had previously been unobservable. The results revealed a nucleation and propagation mechanism wherein the phase transition initiates heterogeneously in La-rich regions and propagates to Bi-rich regions. These findings demonstrate that universal MLFFs are powerful tools for more quickly elucidating complex phase transition phenomena and open up new possibilities for computational science-driven exploration of new materials. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-026-47983-9. Subject terms: Chemistry, Materials science, Mathematics and computing, Physics Introduction While most materials expand upon heating, a small but technologically critical class exhibits the counterintuitive behavior of contracting with temperature, known as negative thermal expansion (NTE). Such anomalous thermal expansion is highly sought after in applications demanding extreme dimensional stability—most notably precision optical systems and semiconductor manufacturing—where nanometer-scale thermally induced distortions can directly degrade alignment, imaging performance, and overlay accuracy. In particular, the continued scaling of device features and the increasing reliance on ultra-high-precision lithography (including EUV optics and masks) amplify the need for near-zero thermal expansion (ZTE) structural components. A widely adopted strategy is to engineer ZTE composites by combining NTE phases with conventional positive-thermal-expansion materials, enabling tailored, application-specific thermal expansion coefficients 1 – 6 . The mechanisms driving NTE are diverse, but a type driven by a first-order phase transition from a large-volume, low-temperature stable phase to a small-volume, high-temperature phase is of particular interest due to the colossal volume changes often involved 1 . Research on high-pressure synthesized perovskite oxides provides representative examples of this phase-transition-driven NTE. For instance, in the PbVO₃ derivatives with PbTiO 3 type structure with a giant c / a ratio, a colossal volumetric NTE exceeding 8% has been reported, induced by a phase transition from a polar tetragonal structure to a nonpolar cubic structure 5 , 7 , 8 . In the BiNiO 3 derivatives, NTE occurs through a phase transition from a divalent Ni phase with a charge-disproportionated state on the Bi site (Bi 3+ /Bi 5+ , Ni 2+ ) to a uniform trivalent Ni phase (Bi 3+ , Ni 3+ ) involving charge transfer between Bi and Ni. By optimizing elemental substitution in this system, the phase transition temperature has been controlled to achieve giant NTE near room temperature 9 , 10 . Thus, phase-transition-driven NTE is a complex phenomenon where multiple degrees of freedom—lattice, charge, and spin—are intricately coupled, making the atomic-level elucidation of its mechanism a critical challenge from both fundamental and applied perspectives. The perovskite oxide BiCoO 3 , the focus of this study, belongs to this family of phase-transition-type NTE parent materials. Under ambient pressure, the stereochemically active 6 s 2 lone pair of Bi 3+ and the Jahn-Teller effect of high-spin Co 3+ lead to a large polar displacement of Co within the CoO 6 octahedra as schematically illustrated as the inset of Fig. 1 (b), resulting in a pyramid-like coordination and a ferroelectric tetragonal structure (space group P 4 mm ) with a significantly elongated c -axis as c / a = 1.27 11 . When a pressure of approximately 3 GPa is applied, it undergoes a spin-crossover phenomenon where Co 3+ ions transform from high-spin to low-spin states 12 . This spin-crossover is accompanied by a transition to a nonpolar paraelectric orthorhombic phase with tilted, isotropic CoO 6 octahedra, resulting in volume contracting colossally by about 13%. This small-volume phase, stable under high pressure, can be stabilized at ambient pressure through elemental substitution. Although irreversible, temperature-induced NTE has been reported in systems where Bi and Co are substituted for Ba and Ti, respectively 15 . Similarly, in the Bi 1− x La x CoO 3 system studied here, substituting La 3+ for Bi 3+ is expected to stabilize the collapsed, high-pressure-like phase at ambient pressure, thereby enabling a thermally driven first-order transition accompanied by NTE 12 , 16 . Fig. 1. Open in a new tab Pressure dependence of the lattice parameters (a) and volume (b) and of BiCoO 3 calculated by MLFF-MD reproducing the experimentally observed polar to non-polar transition accompanied by a volume collapse. From a computational perspective, BiCoO 3 has served as a prototypical system in which a large tetragonal distortion is stabilized by the stereochemically active Bi 6 s 2 lone pair and the high-spin configuration of Co 3+ . First-principles studies have reproduced the pressure-induced collapse toward a lower-volume, nonpolar perovskite structure and emphasized that the structural transformation is intimately coupled to a Co 3+ spin-state change 12 . Ravindran et al. used DFT and fixed-spin-moment analyses to demonstrate magnetic instability (metamagnetism) and a strong magnetovolume effect in BiCoO 3 , emphasizing the intimate coupling between ferroelectric distortion and the Co spin state 13 . Ming et al. further investigated tetragonal BiCoO 3 under pressure by first-principles calculations and reported a pressure-induced Co 3+ magnetic moment collapse (HS LS), accompanied by a first-order isosymmetric structural transition, volume collapse, and an insulator-to-semimetal transition 14 . These studies established BiCoO 3 as a canonical lattice–spin-coupled oxide and provide a crucial microscopic basis for interpreting volume-collapse phenomena in BiCoO 3 -based systems. More recently, Hubberd U corrected density functional theory (DFT + U ) studies on Bi 1− x La x CoO 3 suggested that La substitution progressively stabilizes the low-spin (or reduced-spin) state and lowers the energetic preference for the large-volume polar phase, consistent with an enhanced tendency toward volume collapse upon external stimuli 16 . However, these computational works predominantly address near-0 K energetics or limited structural pathways and therefore do not directly resolve the finite-temperature, real-space dynamics of the first-order transition—such as heterogeneous nucleation, propagation across chemically inhomogeneous regions, and the microscopic origin of partial reversibility/phase coexistence observed experimentally. Bridging this gap requires simulations that simultaneously access large length scales and picosecond–nanosecond dynamics while retaining a first-principles-quality description of the coupled lattice–spin landscape. In particular, ab initio molecular dynamics (AIMD) sampling of nucleation events and interfacial migration is prohibitively expensive for the supercell sizes and time scales required to capture first-order transition kinetics in disordered perovskites. First-principles calculations based on DFT provide a rigorous framework to analyze complex phase transitions at the atomic level, yet their computational cost severely limits accessible time and length scales for first-order transitions involving nucleation and interfacial motion. To overcome this bottleneck, machine-learning interatomic potentials—whose modern development dates back at least to the mid-2000s—have evolved into a practical route for accelerating atomistic simulations by learning DFT-quality potential energy surfaces 17 – 19 . A particularly important recent advance is the emergence of pretrained “universal” or “foundation” models trained on very large, chemically diverse reference datasets, enabling transfer to new materials with little or no system-specific retraining 20 – 22 . In this context, a critical open question is whether such pretrained universal models can faithfully describe transition-metal oxides that undergo discontinuous first-order transitions with colossal volume changes and strong coupling between lattice and spin/electronic degrees of freedom. Here, we address this challenge using CHGNet, a pretrained universal machine-learning force field (MLFF) that incorporates magnetic information, to perform large-scale molecular dynamics simulations of Bi 1− x La x CoO 3 without any additional fine-tuning 22 . This approach not only enables direct visualization of the phase-transition dynamics as a “computational microscope,” but also allows a consistent, apples-to-apples comparison across compositions under an identical simulation protocol—an advantage that is difficult to achieve with conventional first-principles workflows. Recently, the application of MLFFs to NTE materials has gained attraction, offering atomic-scale insights that are difficult to access experimentally. For instance, MLFFs based on the Deep Potential scheme have successfully elucidated the role of polyhedral vibrations in the pressure-induced amorphization of ZrW 2 O 8 23 , the “lateral and tilt” vibrational modes in Zn(CN) 2 24 , and the impact of octahedral tilting on the NTE magnitude in fluorides such as ScF 3 and CaZrF 6 25 . Furthermore, modular development strategies for MLFFs have extended these capabilities to solid solutions like Pb x Sr 1− x TiO 3 , reproducing complex topological textures 26 . However, these pioneering studies typically rely on constructing system-specific models, which necessitate the generation of extensive DFT datasets and dedicated training processes for each target material system. This “train-from-scratch” paradigm remains a bottleneck for high-throughput exploration, particularly when searching for new compositions or investigating substitution effects across a wide range of elements (e.g., exploring various dopants in BiCoO 3 ). This is where “universal” MLFFs pretrained on massive databases, such as CHGNet, offer a distinct advantage by eliminating the need for ad hoc DFT calculations. Yet, the applicability of these universal models to transition metal oxides exhibiting colossal volume changes (> 5%) and discontinuous first-order phase transitions, phenomena deeply coupled with electronic and spin states, has not been sufficiently verified. Unlike the continuous NTE mechanisms seen in framework structures like ZrW 2 O 8 , the phase-transition-driven NTE in Bi-based perovskites involves drastic structural rearrangements and electronic structure changes that challenge the transferability of pretrained models. Therefore, demonstrating that a universal MLFF can reproduce these complex behaviors in a solid solution system without additional training represents a critical step toward the democratization of computational materials design. In this study, we focus on the Bi 1− x La x CoO 3 system and perform Molecular Dynamics (MD) simulations using both CHGNet and first-principles calculations to reproduce the finite-temperature phase transition and elucidate its underlying principles at the atomic level. Results Figure 1 shows the pressure dependence of the lattice parameters and the resulting unit cell volume of BiCoO₃, as calculated by the universal MLFF, CHGNet. In order to reproduce the structural transition from P 4 mm to Pnma phases and the random distribution of substituted elements, the symmetry is lowered to P 1 as detailed in Tables S1–S10 in the Supplementary Information. To clarify the origin of the pressure-induced phase transition, we first analyzed the changes in lattice parameters in detail. In the low-pressure phase, the structure, which is experimentally known to be tetragonal ( P 4 mm , a = b ), showed a slight deviation from this symmetry in our simulation, resulting in distinct a - and b -axis lattice parameters (see Fig. 1 a). To clarify the origin of the symmetry breaking ( a ≠ b ) observed in the MLFF simulation at 0 GPa, we performed static energy calculations for pure BiCoO 3 using both CHGNet and DFT (VASP, Vienna Ab initio Simulation Package). We compared two structures: the experimentally reported tetragonal phase ( a = b ) and the slightly distorted monoclinic-like phase ( a ≠ b ) observed in our MLFF simulation. The DFT energies were corrected using the MaterialsProject2020 Compatibility scheme implemented in pymatgen to ensure consistency with the training data 27 . The calculations revealed a subtle discrepancy in the ground state stability. While DFT correctly predicts the tetragonal phase to be more stable by 13.31 meV/f.u., CHGNet predicts the distorted phase to be more stable by 18.44 meV/f.u. compared to the tetragonal one. This small energy inversion (≒18 meV) in the MLFF’s potential energy surface explains why the simulation settles into a lower-symmetry structure with different a - and b -lattice parameters. However, since the energy difference is minute, the global trend of the pressure-induced phase transition is preserved. Despite this slight symmetry breaking, a significantly elongated c -axis ( c / a = 1.27 at 0 GPa) was clearly observed, reflecting the effect of the 6 s 2 lone pair in Bi 3+ . A very sharp first-order phase transition is observed at 3.7 GPa. At this transition pressure, the c -axis length dramatically contracts from about 4.5 to 4.1 Å, while the a - and b -axes lengths slightly expand. This anisotropic lattice deformation is qualitatively consistent with the experimentally observed structural phase transition to a GdFeO 3 -type structure involving tilting of the CoO 6 octahedra 12 . This selective compression along the c -axis is the primary cause of the colossal volume contraction. As a result of these lattice parameter changes, the volume discontinuously decreases from approximately 62 Å 3 to about 56 Å 3 at the transition pressure. The total volume change, Δ V / V 0 , reaches − 9.7%, In the initial structure optimization, a slight deviation from the experimentally observed tetragonal symmetry ( a = b ) was noted, resulting in a minor discrepancy between the a - and b -axis lattice parameters. This is likely due to the inherent limitations of the universal MLFF. While CHGNet is remarkably powerful, the potential energy surface it predicts for a complex system like BiCoO 3 , which is strongly influenced by Bi 3+ lone pairs and the Jahn-Teller effect of Co 3+ , may contain shallow local minima. This can lead to “systematic softening,” which will be discussed later. Compared to the experimental results, the transition pressure obtained in our simulation (3.7 GPa) overestimates the experimental value (approx. 3.0 GPa), while the volume change (−9.7%) underestimates the experimental value (approx. −13%) 12 . These quantitative discrepancies may be attributed to inherent errors in the universal MLFF and the simulation conditions, as discussed later. However, the crucial point is that the universal MLFF, without any additional training for this complex system, was able to qualitatively reproduce the essential physics of a pressure-induced first-order phase transition in a solid accompanied by volume greatly shrinking. This result strongly suggests that CHGNet can describe the structural energy landscape of BiCoO 3 -derivatives reasonably accurately, providing a solid foundation for investigating the more complex temperature-induced NTE phenomenon. Next, to verify the occurrence of a phase transition upon La substitution, we performed Ab Initio AIMD simulations for Bi 0.5 La 0.5 CoO 3 . Table 1 shows the average volume and the c / a ratio, an indicator of tetragonality, after a 2-ps simulation at each temperature to ensure structural equilibration. In the temperature range from 300 to 450 K, the volume slightly increased from 71.21 to 71.34 Å 3 , exhibiting normal positive thermal expansion. During this interval, the c / a ratio remained large, increasing from 1.265 to 1.284, indicating the stability of the low-temperature tetragonal phase. However, a dramatic change is observed between 450 and 500 K. The volume sharply decreases by approximately 7.0%, from 71.34 to 66.34 Å 3 , confirming a clear NTE. Simultaneously, the c / a ratio approaches 1, decreasing from 1.284 to 1.069, which suggests a transformation to a more isotropic crystal structure. Table 1. Volume and c/a ratio of at finite temperatures from AIMD calculations. Temperature (K) Volume (Å 3 ) c/a ratio 300 71.21 1.265 400 71.23 1.269 450 71.34 1.284 500 66.34 1.069 Open in a new tab This abrupt change in lattice parameters is due to a local atomic coordination change, as shown in Fig. 2 . In the low-temperature phase at 450 K, the Jahn-Teller distortion of Co 3+ persists, resulting in a polar, near-tetragonal structure where Co ions adopt a five-coordinated pyramidal configuration. In contrast, in the high-temperature phase at 500 K, the thermal energy overcomes the stereochemical activity of the 6s 2 lone pairs of Bi 3+ and the Hund coupling stabilizing the high-spin Co 3+ , and the structure transforms into a smaller-volume, higher-symmetry phase where Co ions are centered in six-coordinated octahedral environments. Fig. 2. Open in a new tab Crystal structures of Bi 0.5 La 0.5 CoO 3 from AIMD simulations (at 450 and 500 K). These AIMD calculations confirm at a first-principles level that Bi 0.5 La 0.5 CoO 3 undergoes a first-order phase transition accompanied by distinct local structural changes in the 450–500 K temperature range, resulting in NTE. Next, using the NTE phenomenon confirmed by AIMD as a benchmark, we extended our investigation to a systematic compositional search of Bi 1− x La x CoO 3 using MLFF-MD. Figure 3 (a) shows the temperature dependence of the lattice volume from MLFF-MD simulations with the La substitution systematically varied from 0 to 100%. This comprehensive simulation clearly revealed three important trends: Fig. 3. Open in a new tab Temperature dependence of (a) the lattice volume and (b) the c -axis length of Bi 1− x La x CoO 3 calculated by MLFF-MD. Occurrence of NTE : For compositions with La substitution ranging from 25% to 75%, a first-order phase transition, i.e., NTE, characterized by volume discontinuously and sharply decreasing upon heating, is clearly reproduced. Systematic change in transition temperature : The most significant finding is that the NTE transition temperature systematically shifts to a lower side as the La substitution increases. This suggests that La substitution dilutes and weakens the effect of the Bi 3+ lone pair, enabling the transition to the smaller-volume high-temperature phase to be induced with lower thermal energy. Disappearance of NTE : When the La substitution was 12.5% or less, the volume no longer sharply contracts upon heating; only monotonic positive thermal expansion is seen. Furthermore, these compounds were found to decompose at elevated temperatures instead of undergoing the NTE transition, highlighting the critical role of La substitution in increased stability. To identify the specific symmetry and octahedral tilting pattern of the high-temperature phase, we performed a detailed structural analysis using the MLFF-MD trajectories. Direct observation of snapshots at high temperatures is often hindered by significant thermal fluctuations. Therefore, we calculated the time-averaged structure over a window of 100 fs (100 frames) to extract the equilibrium atomic positions. Figure 4 presents the projected structures of Bi 0.75 La 0.25 CoO 3 before and after the transition. The low-temperature phase (Fig. 4 (a)) clearly maintains the pyramidal coordination of Co with no observable tilting of the polyhedra ( a 0 a 0 c 0 ), preserving the tetragonal P 4 mm -like local symmetry. In contrast, the high-temperature phase (Fig. 4 (b)) exhibits a cooperative rotation of the CoO 6 octahedra. By analyzing the rotation patterns along each axis, we identified an antiphase rotation along the a - and c -axes ( a − , c − ) and an in-phase rotation along the b -axis (b + ). This a − b + c − (or simplified as a − b + a − ) tilting system is the distinct signature of the GdFeO 3 -type orthorhombic perovskite structure (space group Pnma ). This result provides robust computational evidence that the thermally induced phase transition in La-substituted BiCoO 3 involves a symmetry change from P 4 mm to Pnma , consistent with the high-pressure phase of the parent BiCoO 3 . Fig. 4. Open in a new tab Visualization of octahedral tilting and structural phase transition in Bi 0.75 La 0.25 CoO 3 analyzed by time-averaged MLFF-MD trajectories. (a) The low-temperature phase exhibiting a polar structure with 5-coordinated pyramidal Co ions and no octahedral tilting (Glazer notation a 0 a 0 c 0 ), consistent with the tetragonal P 4 mm symmetry. (b) The high-temperature phase exhibiting a non-polar structure and 6-coordinated octahedra with tilting pattern of a − b + c − (or a − b + a − assuming a ≒ c ), corresponds to the Pnma space group (GdFeO 3 -type). Furthermore, a detailed analysis of the temperature and composition dependence of the c -lattice parameter (Fig. 3 (b)) reveals that, similar to the pressure-induced phase transition, the temperature-induced NTE is also primarily driven by the selective contraction of the c -axis. Notably, a clear relationship is observed between the La substitution amount and the lattice parameters of the low-temperature phase. In compositions exhibiting NTE, the a - and b -axis length did not significantly change after the La substitution (see Figures S1 and S2 in the Supplementary Information), whereas the c -axis length markedly tended to contract even below the transition as La content increased. This suggests that substitution with La 3+ relaxes the structural anisotropy caused by the 6 s 2 lone pair of Bi 3+ . To clarify the electronic origin of the observed NTE, we focused on the magnetic moments of the Co ions, a secondary output from the simulations. The CHGNet model used in this study can predict not only the forces and energy but also the magnetic moment at each atomic site. Figure 5 shows the temperature and composition dependence of the average magnetic moment on the Co sites obtained from the MLFF-MD simulations. Fig. 5. Open in a new tab Temperature and composition dependence of the average magnetic moment on Co sites obtained from MLFF-MD simulations. The results clearly show a strong correlation between the NTE and the spin state of the Co ions. In compositions where NTE is clearly observed (0.25 ≤ x ≤ 0.75), the magnetic moment decreases discontinuously from approximately 3.0 µ B to about 2.6–2.7 µ B , perfectly coinciding with the transition temperatures where the volume sharply contracts. This behavior corresponds to the spin-crossover phenomenon, where the high-spin state Co 3+ ions, stable at low temperatures, transition to a low-spin or intermediate state upon heating. The difference in the magnetic moment below and above the transition temperature decreases as the La substitution decreases, most probably because the thermally activated low- (high-) spin state is mixed to the ground state high- (low-) spin state. The magnetic moment of approximately 2.6–2.7 µ B in the post-transition phase is too large to be attributed to a low-spin (LS, S = 0) state. Instead, this value strongly suggests a spin-crossover to a phase dominated by the intermediate-spin (IS, S = 1) configuration of Co³⁺. This interpretation is consistent with detailed studies on the related perovskite LaCoO₃. Through GGA + U calculations, Pandey et al. concluded that the electronic structure of LaCoO₃ at room temperature is best described by a dominant IS state configuration 28 . They also highlighted the strong hybridization between Co 3 d and O 2 p orbitals within the perovskite octahedra, which differentiates the system from a purely ionic picture. Therefore, we interpret the observed magnetic moment not as a covalent LS state, but as clear evidence of a transition to a high-temperature phase with a dominant IS character. We note that the “magnetic moments” provided by CHGNet are model-predicted local moments and do not uniquely map onto an integer spin quantum number in a covalent transition-metal oxide. Co–O hybridization and finite-temperature fluctuations can partially quench and redistribute spin density, such that the moment magnitude should be interpreted as an indicator of a reduced-spin character rather than a strict HS/IS/LS label. Accordingly, we describe the change as an HS → IS-like (reduced-spin-dominated) transition consistent with the concurrent structural collapse. This change in spin state directly impacts the lattice structure. High-spin Co 3+ ion is Jahn-Teller active, which distorts the octahedron and stabilizes the pyramidal coordination. In the intermediate- or low-spin state, however, the Jahn-Teller effect vanishes, allowing for a more symmetric octahedral coordination. This change in local structure is the driving force behind the lattice volume contraction, i.e., NTE. On the other hand, in compositions where NTE was not observed ( x ≤ 0.125 and x = 1.0), the absence of such a clear spin-crossover transition supports this mechanism. From the analysis above, our MLFF-MD simulations strongly suggest that the NTE in the Bi 1− x La x CoO 3 system is a phenomenon in which the lattice and spin degrees of freedom are strongly coupled, accompanied by a temperature-induced spin-crossover of Co 3+ ions. This demonstrates that a universal MLFF is not merely a tool for structural prediction but also an extremely effective method for capturing the essential physics, such as changes in electronic states, behind complex phase transitions. To investigate the reversibility of the phase transition and the hysteresis behavior, we performed cooling simulations for the x = 0.25 composition following the heating run. Figure 6 displays the evolution of total energy and unit cell volume during the thermal cycle. Interestingly, the simulation revealed a “one-way” irreversible phase transition. As shown in Fig. 6 (a) , once the system transforms into the small-volume phase (GdFeO 3 -type) around 600 K during heating, it does not revert to the original large-volume phase upon cooling. Instead, it maintains the compact, octahedral-coordinated structure down to low temperatures. The energy profile (Fig. 6 (b) ) provides a clear thermodynamic explanation for this behavior: the potential energy of the post-transition phase is lower than that of the initial pre-transition phase at the same temperature. This indicates that, within the potential energy surface described by CHGNet, the experimentally observed low-temperature PbTiO 3 -type phase (tetragonal) corresponds to a metastable local minimum, while the GdFeO 3 -type phase (orthorhombic) represents the global minimum for this composition. This irreversibility may be partly exaggerated by the finite simulation time and supercell size, which can hinder the nucleation of the reverse transition. As will be discussed in the next section, experimental measurements on La-substituted samples ( x = 0.2) exhibit a reversible NTE, while the fraction of the orthorhombic phase changes only 25 to 57% on heating and to 32% on cooling, not from 0 to 100%. We note that the effective heating/cooling rates in MD are orders of magnitude faster than experimental protocols, which can strongly bias apparent hysteresis and reversibility in first-order transitions. Therefore, quantitative comparison of hysteresis widths should be made with caution, whereas the qualitative nucleation–growth mechanism remains robust. In addition, reproducing persistent two-phase coexistence in real polycrystalline samples likely requires explicit treatment of mesoscale domain structures, interfacial strain accommodation, and defect/grain-boundary pinning, which are not represented in the present periodic atomistic supercells. Accordingly, the current simulations are better suited to resolving intrinsic local transition pathways than to quantitatively reproducing the full extent of experimentally observed phase coexistence. Fig. 6. Open in a new tab Irreversibility of the phase transition in B 0.75 La 0.25 CoO 3 revealed by heating and cooling MLFF-MD simulations. (a) unit cell volume and (b) total potential energy for the x = 0.25 composition during heating (orange solid line) and cooling (blue dashed line) cycles. While a sharp first-order transition from the large-volume tetragonal phase to the small-volume orthorhombic phase is observed around 600 K upon heating, the reverse transition does not occur upon cooling. The system remains in the small-volume phase down to low temperatures. The energy profile (b) indicates that the post-transition phase is energetically more stable (lower energy) than the initial phase at 0 K, suggesting the irreversible nature of this transition in the simulation. To verify the “one-way” transition behavior suggested by the simulation, we compared these results with experimental synchrotron X-ray diffraction data for a compositionally similar sample, Bi 0.8 La 0.2 CoO 3 ( x = 0.20). Figure 7 shows the temperature dependence of the unit cell volumes and the phase fraction. Fig. 7. Open in a new tab Experimental thermal evolution of lattice parameters and phase fraction in Bi 0.8 La 0.2 CoO 3 . Temperature dependence of the unit cell volume and phase fraction for Bi 0.8 La 0.2 CoO 3 obtained from synchrotron X-ray diffraction experiments. The blue squares and red triangles represent the volumes of the tetragonal (PbTiO 3 -type) and orthorhombic (GdFeO 3 -type) phases, respectively. The black circles indicate the phase fraction of the orthorhombic phase. The green circles show the weighted average volume ( V PC ). Arrows indicate the heating and cooling processes. The experimental system exists as a two-phase mixture of the large-volume tetragonal phase (PbTiO 3 -type) and the small-volume orthorhombic phase (GdFeO 3 -type) even at 100 K. Upon heating, the fraction of the orthorhombic phase (black line) increases, from 25 to 57%, leading to the NTE behavior observed in the average volume (green line). Crucially, during the subsequent cooling process, the phase fraction exhibits significant hysteresis. The system does not fully revert to its initial state; instead, the fraction of the orthorhombic phase (32%) remains higher than it was before heating. This experimental observation qualitatively aligns with the irreversibility predicted by our MLFF-MD simulation. The simulation indicated that the orthorhombic phase is energetically more stable than the tetragonal phase over a wide temperature range for La-substituted compositions. In the idealized simulation box with periodic boundary conditions, this led to a complete lock-in to the stable orthorhombic phase. In the real polycrystalline sample, the transformation is spatially constrained by domain walls and inter-granular strain, preventing a complete transformation and resulting in the observed hysteresis and mixed-phase state. Nevertheless, both simulation and experiment point to the same conclusion: La substitution fundamentally stabilizes the collapsed GdFeO 3 -type structure, making the NTE transition partially or wholly irreversible depending on the microstructural constraints. Figure 8 shows the dependence of the NTE transition temperature on the La substitution amount, as obtained from MLFF-MD. Experimental results shown in the Fig. 7 indicate that for a Bi 0.8 La 0.2 CoO 3 , the NTE transition is observed over a wide temperature range from approximately 480 K to 700 K. The closest composition in our simulation, x = 0.25, predicted a transition temperature of 725 K, which is considerably higher than the onset temperature observed experimentally. For the x = 0.5 composition, the MLFF prediction (approx. 350 K) was 100–150 K lower than the AIMD prediction (450–500 K). These results indicate that the universal MLFF has quantitative challenges in predicting the absolute values of transition temperatures. Nevertheless, it can be said that MLFF possesses sufficient capability to correctly capture the chemical trends, such as whether NTE will occur with La substitution and how the transition temperature changes as the substitution amount increases. Fig. 8. Open in a new tab NTE transition temperatures calculated by MLFF-MD. The true value of the MLFF-MD simulation used in this study extends beyond the prediction of thermodynamic properties. It functions as a “computational microscope” that enables the dynamic process of how a phase transition proceeds to be directly observed, with the movement of each individual atom tracked. Figure 9 shows representative snapshots of atomic configurations within ± 20 fs of the moment the NTE phase transition occurs in the heating simulation of Bi 0.5 La 0.5 CoO 3 . Fig. 9. Open in a new tab Visualization of the phase transition dynamics in the MLFF-MD simulation on heating. Green, purple, blue, and red atoms show La, Bi, Co, and O, respectively. (a) Before, (b) during, and (c) after the phase transition. Before the transition : In the low-temperature phase, the entire system adopts a near-tetragonal structure elongated along the c -axis, with Co 3+ ions in a 5-coordinated pyramidal configuration. During the transition : Most notably, the snapshot during the transition clearly shows that the transformation to the small-volume high-temperature phase (with Co 3+ in 6-coordinated octahedral structures) does not occur uniformly throughout the system. It clearly captures the transition propagating from the La 3+ -rich regions (which lack lone pairs) to the Bi 3+ -rich regions. After the transition : Ultimately, the entire system completely transforms into the high-temperature phase with octahedral coordination. This visualization provides a new perspective for understanding complex phase transition mechanisms. Experimental techniques, typified by in-situ transmission electron microscopy (TEM) observation, are powerful methods able to capture structural changes under heating with atomic resolution. However, the time resolution of conventional in-situ TEM observations, which is typically limited by the video frame rate to the order of milliseconds, is insufficient to resolve the femtosecond-to-picosecond timescale on which the elementary processes of a phase transition, such as atomic vibrations, occur. Therefore, the information obtained from such techniques is essentially a time-averaged structure over these ultrafast events. MLFF-MD simulation can track the dynamic behavior of individual atoms in this ultrafast time domain, which is difficult to observe experimentally. The heterogeneous propagation behavior, nucleating from La-rich regions, as revealed here, is a transient phenomenon that cannot be captured from the averaged information. This compositional bias is plausibly driven by local chemical-pressure and bond-stiffness effects introduced by La substitution, which preferentially stabilizes the collapsed (reduced-volume) coordination environment. In this picture, La-rich regions lower the effective nucleation barrier for the small-volume phase, providing a microscopic link between configurational disorder and the kinetics of the first-order transition. Furthermore, the ability to perform such systematic dynamic analyses for numerous compositions at a realistic computational cost is a significant advantage over experiments. Therefore, MLFF-MD not only enables high-speed screening but also plays a complementary role to experimental observations, demonstrating its potential to become a new analytical method for elucidating the fundamental elementary processes of complex solid-state reactions. To assess the influence of finite size effects on the phase transition, we performed additional simulations using a larger supercell (1080 atoms) and compared the results with those from the supercell (320 atoms) (see Supplementary Figure S3). The results confirm that the NTE behavior is robust against system size changes. Two notable features were observed in the larger system: Reduced Fluctuations : The standard deviation of the lattice parameters and volume is significantly smaller in the system, indicating that the thermodynamic properties are well-converged. Shift in Transition Temperature : The NTE transition temperature in the system was approximately 80 K lower than in the system. This shift is physically consistent with the nucleation and growth mechanism characteristic of first-order phase transitions. In a larger simulation cell, the probability of a critical nucleus forming within a given time at a given temperature increase, allowing the transition to proceed at a lower temperature during the heating process. Note that while the magnitudes of the a - and b -axis lattice parameters in the high-temperature phase appear swapped between the two system sizes, this is merely due to the crystallographic equivalence of these axes in the transition from the tetragonal ( P 4 mm ) to the orthorhombic ( Pnma ) structure, where the tilting direction can spontaneously select either axis. Discussion The quantitative deviation in the transition temperature observed in this simulation is thought to stem from the inherent limitations of universal MLFFs. Several recent studies have pointed out that many universal MLFFs, including CHGNet, exhibit a tendency known as “systematic softening,” where they underestimate the interatomic forces in high-energy regions 29 . The primary cause of this phenomenon is the training dataset for the universal model (in CHGNet’s case, the Materials Project DFT relaxation trajectories 30 biased towards near-equilibrium structures around energy minima, with insufficient sampling of high-energy regions such as transition states and highly strained structures. The phase transition temperature is determined by the point where the free energies of the two phases become equal, and the contribution of vibrational entropy is significant to the free energy. If the potential energy surface (PES) is soft (i.e., has a small curvature), phonon frequencies are underestimated, which in turn affects the evaluation of vibrational entropy and can shift the predicted transition temperature. The underestimation of the transition temperatures seen in this study can be regarded as a direct manifestation of this systematic softening. To further validate the reliability of the MLFF-MD simulation and assess the impact of this systematic softening, we performed single-point DFT calculations on 400 snapshots extracted from the MD trajectory of Bi 0.5 La 0.5 CoO 3 during the heating process. Figure 10 compares the potential energies predicted by CHGNet with those calculated by DFT (VASP). Crucially, both methods exhibit a clear and simultaneous drop in total energy around snapshot 4,000, which corresponds to the experimentally and computationally observed phase transition temperature. This agreement confirms that the phase transition reproduced by CHGNet is not a mere artifact of the machine learning model but is physically driven by the stabilization of the high-temperature phase, consistent with first-principles energetics. However, as shown in the high-temperature region (around and after snapshot 4,000), the energy difference between CHGNet and DFT gradually increases. This trend provides that MLFF tends to under- or overestimate the energy of structures with large thermal displacements. Nevertheless, the fact that the distinct energy jump associated with the first-order transition is qualitatively captured supports the validity of using this universal MLFF to elucidate the transition mechanism, even if quantitative discrepancies in the transition temperature persist. Fig. 10. Open in a new tab Validation of the potential energy surface along the MD trajectory. Comparison of the total energies predicted by CHGNet (blue circles) and calculated by DFT (red squares) for 400 snapshots extracted from the AIMD trajectory of Bi 0.5 La 0.5 CoO 3 during the heating process. The DFT energies were corrected using the Pymatgen MaterialsProject2020Compatibility module. To further quantify the discrepancies observed in Fig. 10 , we performed a systematic comparison of total energies and atomic forces between CHGNet and DFT for 400 snapshots extracted from the AIMD trajectory (see Supplementary Figure S4). The analysis yields a Mean Absolute Error (MAE) of 1.676 eV for the total energy of the 40-atom supercell (approx. 42 meV/atom) and 0.458 eV/Å for the atomic forces. As shown in the parity plots, CHGNet consistently predicts lower total energies (i.e., more stable structures) compared to DFT. This tendency can be attributed to the characteristics of the training dataset used for CHGNet (MPtrj from the Materials Project), which is predominantly composed of relaxed, near-equilibrium structures. Consequently, for high-temperature snapshots containing significant thermal distortions—configurations that are energetically unstable and under-represented in the training data—the model likely underestimates the energy penalty associated with such instability. Despite this offset in absolute values due to the training data bias, the model successfully captures the relative energy jump associated with the phase transition, confirming its capability to reproduce the thermodynamic trends required for this study. On the other hand, a likely reason for CHGNet’s ability to capture the qualitatively correct physics is that its architecture explicitly incorporates magnetic moment (spin) information 22 . Since the phase transition in the BiCoO 3 system is closely related to the change in the spin state of the Co ions, this feature likely enabled the model to somewhat learn the coupling between the electronic state and the structure, leading to its qualitative success. Furthermore, the computational efficiency of the MLFF-MD method employed in this study is remarkable. An AIMD calculation using a 40-atom 2 × 2 × 2 supercell required approximately 20 h to run a 2-ps simulation at a single fixed temperature. In contrast, an MLFF-MD calculation using a 320-atom 4 × 4 × 4 supercell (8 times larger than AIMD) and simulating a continuous temperature ramp from 1 K to 1000 K (1000 points) took only 8 h. A simple throughput comparison shows that the MLFF achieved a computational speed over 20,000 times faster than AIMD. This overwhelming speed-up makes comprehensive simulations (screenings) of material behavior over wide compositional and temperature ranges, which were previously impractical due to computational cost, a tangible reality. However, the current universal training data is insufficient to quantitatively and accurately describe the subtle energy changes associated with spin-crossover, which is considered another source of quantitative error. For future high-precision predictions, an approach that uses the universal model as a base and then fine-tunes it with a small amount of high-accuracy DFT data specific to the target system (such as the AIMD trajectories generated in this study) is suggested to be effective 31 . Conclusion In this study, we applied large-scale molecular dynamics (MD) simulations using the Crystal Hamilton Graph Network (CHGNet), a pretrained universal machine learning force field (MLFF), to elucidate the negative thermal expansion (NTE) phenomenon in the perovskite oxide Bi 1− x La x CoO 3 . The main achievements are as follows: Methodology Validation and NTE Reproduction : The universal MLFF, CHGNet, qualitatively reproduced the giant volume contraction induced by pressure in BiCoO 3 . Furthermore, it successfully reproduced the experimentally observed pressure-induced phase transition in BiCoO 3 and NTE in the Bi 1− x La x CoO 3 , and indicated the chemical trend, where the NTE transition temperature decreases as La substitution increases. Elucidation of Atomic-Level Dynamic Mechanism : By leveraging the overwhelming computational efficiency of MLFF-MD (over 20,000 times faster than AIMD), we visualized the dynamic process of the NTE phase transition at the atomic level. As a result, we uncovered the nucleation and propagation behavior, where the phase transition initiates heterogeneously in La-rich regions and propagates to Bi-rich regions. These achievements concretely demonstrate that universal MLFF-MD simulation can serve as a powerful “computational microscope” for predicting the behavior of functional inorganic materials with complex phase transitions and for elucidating atomic-level mechanisms that are unobservable in experiments. This research not only embodies a new paradigm for materials science research driven by computational and data sciences but also provides important feedback for developing next-generation, higher-precision foundational models. Methods Machine learning force field molecular dynamics simulations In this study, the CHGNet (ver. 0.4.0) was employed as the MLFF 22 .CHGNet is a universal MLFF pretrained on the extensive first-principles calculation database of the Materials Project 30 . It was used in this research without any additional training (fine-tuning). MD simulations were performed using the open-source physics simulation library ASE (ver. 3.26.0) 32 with PyTorch (ver. 2.7.1) as its backend. All MLFF-MD calculations were executed on a single NVIDIA RTX 4090 GPU. A 4 × 4 × 4 supercell (320 atoms) created by expanding the unit cell of BiCoO 3 and Bi 1− x La x CoO 3 solid solutions was used as the simulation cell. For Bi 1− x La x CoO 3 , Bi and La were distributed randomly at the A-sites at each target composition. We used one representative random configuration per composition to maintain a consistent protocol across the composition sweep, and the same A-site arrangement was retained throughout the corresponding heating/cooling simulations. Prior to the MD simulations, the initial structure for each composition was optimized by using the BFGS algorithm, a quasi-Newton method, until the residual force on each atom was less than 0.001 eV/Å. For the subsequent evaluation of thermal expansion behavior, the NpT ensemble (isothermal-isobaric ensemble) was employed to maintain constant pressure and temperature. Temperature was controlled using a Nosé-Hoover thermostat (time constant = 25 fs) 33 , and pressure was controlled with a Nosé-Hoover barostat. A time step of 1.0 fs was used. In the pressure-induced phase transition simulations, the temperature was kept at 10 K while the pressure was incrementally increased from ambient pressure (1 bar) to 5 GPa in steps of 0.1 GPa. The structure was held for 2 ps at each pressure point to ensure sufficient relaxation. For the temperature-induced phase transition simulations, the pressure was maintained at ambient pressure (1 bar), and the temperature was ramped up and down from 0 K to 1000 K in 1 K intervals. At each temperature point, an MD calculation of 2000 steps (2 ps) was performed to equilibrate the system, and the lattice parameters and volume were recorded. In the heating simulations, the final structure from the previous temperature step was used as the initial structure for the next step to ensure continuity. First-principles molecular dynamics simulations To compare and validate the MLFF-MD results, AIMD calculations were also performed using the VASP 34 . A 2 × 2 × 2 supercell of Bi 0.5 La 0.5 CoO 3 (40 atoms: Bi 4 La 4 Co 8 O 24 ) was used for the simulation cell. The Bi and La atoms were placed on the A-sites by a random assignment at the target composition. The same A-site arrangement was used consistently for all AIMD runs reported here. The initial structure for AIMD was obtained by DFT pre-relaxation, as detailed below, to ensure a well-defined starting point for NpT dynamics. The initial structure for AIMD was obtained by fully relaxing the same 2 × 2 × 2 supercell at the DFT(+ U ) level prior to NpT dynamics. The Perdew-Burke-Ernzerhof (PBE) functional was used to describe the exchange-correlation energy of electrons. To account for the strong electron correlation among the 3 d electrons of Co atoms, the GGA + U method (Dudarev’s approach) 35 , 36 was applied, with an effective U parameter ( U eff ) of 3.32 eV set for the Co d-orbitals. The cutoff energy for the plane-wave basis set was 520 eV, and the first Brillouin zone was sampled using a 3 × 3 × 2 k-point mesh centered at the Gamma point. The simulations were conducted in the NpT ensemble, with temperature and pressure controlled by using Langevin dynamics. Calculations were performed at temperatures of 300, 400, 450, and 500 K for 2 ps (2000 steps) with a time step of 1.0 fs. Structural information for the equilibrium state was obtained from the trajectory of the final 1 ps. Supplementary Information Below is the link to the electronic supplementary material. Supplementary Material 1 (1.9MB, pdf) Acknowledgements This work was partially supported by JSPS KAKENHI (grant no. JP24H00374), JST-CREST (JPMJCR22O1), the Kanagawa Institute of Industrial Science and Technology and Design and Engineering by Joint Inverse Innovation for Materials Architecture project of the Ministry of Education, Culture, Sports, Science and Technology. Author contributions S.W. performed the MD simulations and analyzed the data. Y.S, K.T. and M.A performed the high-pressure synthesis and synchrotron X-ray diffraction experiments. H.D. and M.A. conceived the project and supervised the work. All authors contributed to writing the manuscript. Funding This work was partially supported by JSPS KAKENHI (grant no. JP24H00374), JST-CREST (JPMJCR22O1), the Kanagawa Institute of Industrial Science and Technology and Design and Engineering by Joint Inverse Innovation for Materials Architecture project of the Ministry of Education, Culture, Sports, Science and Technology. Data availability All data generated or analyzed during this study are included in this published article and its supplementary information files. Declarations Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Contributor Information Shogo Wakazaki, Email: wakazaki.s.80d4@m.isct.ac.jp. Masaki Azuma, Email: mazuma@msl.titech.ac.jp. References 1. Takenaka, K. Negative thermal expansion materials: Technological key for control of thermal expansion. J. Phys. Condens. Matter 24 , 013001 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 2. Azuma, M. Phase transition type giant negative thermal expansion materials. J. Ceram. Soc. Jpn. 133 , 450–454 (2025). [ Google Scholar ] 3. Chen, J., Hu, L., Deng, J. & Xing, X. Negative thermal expansion in functional materials: Controllable thermal expansion by chemical modifications. Chem. Soc. Rev. 44 , 3522 (2015). [ DOI ] [ PubMed ] [ Google Scholar ] 4. Takenaka, K. Mater. Trans. 65 , 243 (2024). [ Google Scholar ] 5. Kanaoka, M. et al. Removal properties of low-thermal-expansion materials with rotating-sphere elastic emission machining. Sci. Technol. Adv. Mater. 8 , 170–172 (2007). [ Google Scholar ] 6. Carapella, A. P., Duran, C. A., Hrdina, K. E., Sears, D. J. & Tingley, J. E. ULE® glass for EUVL applications, a fictive temperature correlation. J. Non-Cryst. Solids 357 , 37–42 (2013). [ Google Scholar ] 7. Yamamoto, H. et al. Colossal negative thermal expansion in electron-doped PbVO₃. Angew Chem. Int. Ed. 57 , 8170–8173 (2018). [ DOI ] [ PubMed ] [ Google Scholar ] 8. Nishikubo, T. et al. Polar–Nonpolar Transition-Type Negative Thermal Expansion with 11.1% Volume Shrinkage by Design. Chem. Mater. 35 , 870–878 (2023). [ Google Scholar ] 9. 9. Azuma, M. et al . Giant negative thermal expansion in BiNiO₃. Nat. Commun. 2 , 347 (2011)., Nabetani K, Appl. Phys. Lett. , 106 , 061912 (2015). [ DOI ] [ PMC free article ] [ PubMed ] 10. Nishikubo, T. et al. Enhanced negative thermal expansion induced by simultaneous charge transfer and polar–nonpolar transitions. J. Am. Chem. Soc. 141 , 19397–19403 (2019). [ DOI ] [ PubMed ] [ Google Scholar ] 11. Belik, A. A. et al. Neutron Powder Diffraction Study on the Crystal and Magnetic Structures of BiCoO 3 . Chem. Mater. 18 , 798–803 (2006). [ Google Scholar ] 12. Oka, K. Pressure-induced spin-state transition in BiCoO₃. J. Am. Chem. Soc. 132 , 9438–9443 (2010). [ DOI ] [ PubMed ] [ Google Scholar ] 13. Ravindran, P. et al. Magnetic-instability-induced giant magnetoelectric coupling. Adv. Mater. 20 , 1353–1356 (2008). [ Google Scholar ] 14. Ming, X. et al. Pressure-induced magnetic moment collapse and insulator-to-semimetal transition in BiCoO 3 . J. Phys. Condens. Matter 21 , 295902 (2009). [ DOI ] [ PubMed ] [ Google Scholar ] 15. Pan, Z. et al. Pronounced negative thermal expansion in lead-free BiCoO₃-based ferroelectrics triggered by the stabilized perovskite structure. Chem. Mater. 31 , 6187–6192 (2019). [ Google Scholar ] 16. Benayad, N., Djermouni, M., Zaoui, A., Kacimi, S. & Boukortt, A. Spin state transition in Bi 1−x La x CoO 3 perovskite alloys: DFT + U study. Mater. Chem. Phys. 207 , 507–512 (2018). [ Google Scholar ] 17. Behler, J. & Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett. 98 , 146401 (2007). [ DOI ] [ PubMed ] [ Google Scholar ] 18. Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. Phys. Rev. Lett. 104 , 136403 (2010). [ DOI ] [ PubMed ] [ Google Scholar ] 19. Unke, O. T. et al. Machine learning force fields. Chem. Rev. 121 , 10142–10186 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 20. Chen, C. & Ong, S. P. A universal graph deep learning interatomic potential for the periodic table. Nat. Comput. Sci. 2 , 718–728 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 21. Batatia, I. et al. MACE: Higher order equivariant message passing neural networks for fast and accurate force fields. Adv. Neural Inf. Process. Syst. 35 , 11487–11501 (2022). [ Google Scholar ] 22. Deng, B. et al. CHGNet as a pretrained universal neural network potential for charge-informed atomistic modelling. Nat. Mach. Intell. 5 , 1031–1041 (2023). [ Google Scholar ] 23. He, R., Wu, H., Lu, Y. & Zhong, Z. Origin of negative thermal expansion and pressure-induced amorphization in zirconium tungstate from a machine-learning potential. Phys. Rev. B 106 , 174101 (2022). [ Google Scholar ] 24. Shi, G., Zhong, Z., Su, Y., He, R. & Zhang, C. Investigating the atomic-level mechanisms of negative thermal expansion in Zn(CN)2 using machine learning-enabled molecular dynamics. Phys. Scr. 100 , 066002 (2025). [ Google Scholar ] 25. Zhao, K. et al. Octahedral tilt distortion in negative thermal expansion in the fluorides CaZrF 6 and ScF 3 . Phys. Rev. B 110 , 064322 (2024). [ Google Scholar ] 26. Wu, J., Yang, J., Ma, L., Zhang, L. & Liu, S. Modular development of deep potential for complex solid solutions. Phys. Rev. B 107 , 144102 (2023). [ Google Scholar ] 27. Jain, A. et al. Formation enthalpies by mixing GGA and GGA calculations. Phys. Rev. B . 84 , 045115 (2011). [ Google Scholar ] 28. Pandey, S. K. et al. Investigation of the spin state of Co in LaCoO 3 at room temperature: Ab initio calculations and high-resolution photoemission spectroscopy of single crystals. Phys Rev. B 77 , 045123 (2008). 29. Deng, B. et al. Systematic softening in universal machine learning interatomic potentials. npj Comput. Mater 11 , 9 (2025). 30. Jain, A. et al. The materials project: A materials genome approach to accelerating materials innovation. APL Mater. 1 , 011002 (2013). [ Google Scholar ] 31. Chen, P.-Y., Shibata, K. & Mizoguchi, T. High precision machine learning force field development for BaTiO3 phase transitions, amorphous, and liquid structures. APL Mach. Learn. 3 , 036115 (2025). [ Google Scholar ] 32. Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. J. Phys.: Condens. Matter 29 , 273002 (2017). [ DOI ] [ PubMed ] [ Google Scholar ] 33. Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. J. Chem. Phys. 81, 511–519 ; Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985). (1984). [ DOI ] [ PubMed ] 34. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B . 54 , 11169–11186 (1996). [ DOI ] [ PubMed ] [ Google Scholar ] 35. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77 , 3865–3868 (1996). [ DOI ] [ PubMed ] [ Google Scholar ] 36. Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA + U study. Phys. Rev. B 57 , 1505–1509 (1998). [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplementary Material 1 (1.9MB, pdf) Data Availability Statement All data generated or analyzed during this study are included in this published article and its supplementary information files. Articles from Scientific Reports are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (2.7 MB) Cite Collections Permalink RESOURCES Similar articles Cited by other articles Links to NCBI Databases Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top

Related documents

• Higher-order nonadiabaticity governs the temperature dependence of the phonon spectrum#464617
• Conformational landscapes in cryo-electron tomography data based on molecular dynamics simulations#1158053
• Chem Lab Simulations, No. 1 and No. 2.#333279
• Observation of conformal field theory spectra in a quantum simulator#677093
• Low-Temperature Large Deviations for the Two-Dimensional Wick-Ordered Cubic Wave Equation#986457