The origin of the anomalous expansion of the first peak in the radial distribution function during the rapid solidification of tantalum metal

Abstract

The cause of the anomalous shift in the first maximum peak of radial distribution functions (RDFs) with decreasing temperature in metallic melts and glasses remains highly controversial. In this study, we show that the first RDF peak exhibits anomalous expansion as the temperature decreases during the non-equilibrium solidification (γ₁ = 1 × 10¹⁰ K s⁻¹ and γ₂ = 1 × 10¹¹ K s⁻¹) of liquid tantalum. This behavior is primarily due to alterations in both the geometric and electronic structures of the system. In terms of geometric structure, for example, at the cooling rate of γ₁, the system forms a significant number of cage-like icositetrahedral Voronoi polyhedra (0,0,12,2) and standard icosahedral Voronoi polyhedra (0,0,12,0) at low temperatures. These Voronoi polyhedra have longer bond lengths and lower binding energies compared to their high-temperature counterparts. Furthermore, these Voronoi polyhedra nest together, forming a stable Ta₂₆-C_2v atomic configuration with minimal changes in bond lengths. This unique geometric arrangement contributes fundamentally to the anomalous expansion of the first peak of the RDF. Regarding the electronic structure, the temperature influences the interactions between Ta atoms. At higher temperatures, the electronic localization functions (ELFs) and the Mulliken bond overlap populations (Q_i–j) are significantly increased, leading to stronger electronic interactions and a denser arrangement of nearest-neighbor atoms with shorter bond lengths. Consequently, the combined effects of geometric and electronic structural changes during non-equilibrium solidification could explain the anomalous expansion of the first peak of the RDF.

---

1. Introduction

It is well-established that most materials undergo thermal expansion when heated due to anharmonic vibrations in their atomic lattices.¹ This results in a positive temperature coefficient of volumetric expansion, meaning the material's volume increases as the temperature rises.¹ Thermal expansion is particularly pronounced during phase transitions, such as the transition from solid to liquid or gas.¹ However, the challenges posed by positive thermal expansion (PTE) are common in practical applications, especially in fields like electronics and precision instrumentation. For example, significant temperature fluctuations can lead to measurement inaccuracies, and mismatched thermal expansion coefficients between materials in electronic packaging can cause mechanical failure.^2,3 To mitigate these issues, materials exhibiting minimal or zero thermal expansion, or even negative thermal expansion (NTE), offer promising solutions for maintaining dimensional stability, with significant implications for advanced technological applications.^1–3 The study of NTE materials has a long history, beginning with the discovery of Invar alloys,⁴ which exhibit nearly zero thermal expansion over a wide temperature range. This discovery, made by C. E. Guillaume in 1897 with the Fe₆₅Ni₃₅ alloy,⁴ introduced the “Invar effect”, now synonymous with the exceptional NTE behavior of certain materials. NTE materials can be classified as functional materials, including Fe–Ni alloys,^1,4 oxides,⁵ nitrides,⁶ cyanides,⁷ fluorides,⁸ and metal–organic frameworks (MOFs).⁹ From a mechanistic perspective, NTE materials can generally be divided into two categories: those driven by electronic effects,¹ such as magnetic and ferroelectric materials,^1,4 and those influenced by phonon interactions,¹ such as open-framework structures.¹⁰

For structural materials, many studies have shown that, during the non-equilibrium or rapid solidification of liquid elemental metals, the volume–temperature relationship generally follows the trend of PTE rather than NTE. Despite the absence of NTE behavior (often referred to as the “Invar effect”),^1–10 an interesting anomaly emerges when analyzing the radial distribution function (RDF) of the atomic structure at different temperatures.^11–16 Specifically, as the temperature decreases, the first peak of the RDF unexpectedly shifts to a larger r value instead of contracting to a smaller r value. In other words, the distance between the first-nearest neighbor atoms does not follow the conventional pattern of “thermal expansion upon heating and thermal contraction upon cooling”, but instead exhibits a reverse trend: thermal contraction upon heating and thermal expansion upon cooling. For example, as early as a decade ago, Lou et al. reported an anomalous shift in the position of the first peak in the RDF for pure metallic melts of Al, Ni, Cu, Ag, Au, and In.¹¹ They suggested that this anomaly arises because high-coordination polyhedral clusters tend to transform into low-coordination clusters with increasing temperature.¹¹ Gangopadhyay et al. proposed that the anomalous thermal contraction of the first peak of RDF could be attributed to qualitative changes in coordination numbers and bond lengths as clusters in liquids break up with increasing temperature.¹² However, the conclusion that the contraction of the first peak of RDF is caused by the first coordination polyhedron remains highly controversial.^11–15 For instance, Sukhomlinov et al.¹³ reported that their simulations reproduced a shift of the first peak in the RDF to smaller r at higher temperatures, but a more refined analysis revealed that the ‘true’ bond lengths increased with temperature in all investigated systems. These results align with conclusions drawn by Ding et al., who convincingly demonstrated in 2014 that the contraction could be attributed to the asymmetry of the first maximum peak in the RDFs.¹⁴ Three years later, Ding et al.¹⁵ further explained that the contraction of the nearest-neighbor distance in metallic melts upon heating is due to an increase in the skewness of RDF and a thermally exacerbated asymmetric distribution of neighboring atoms.¹⁵ They identified three main effects influencing the first peak positions in RDFs for metallic materials: atomic structure, atomic volume, and the asymmetry of neighboring atom redistribution.¹⁵ Clearly, the cause of the anomalous shift in the first peak of the RDF remains highly controversial for both metallic melts and solids.

To resolve this controversy and further elucidate the underlying cause of the anomalous expansion of the first peak in the RDF of liquid metals, this study focuses on the refractory metal Ta. The choice of Ta is motivated by a landmark achievement in 2014 by Zhong et al.¹⁶ who demonstrated that several pure refractory body-centered cubic (BCC) metals, including Ta, V, W, and Mo, could be successfully vitrified into metallic glasses by achieving an exceptionally high liquid-quenching rate of 10¹³–10¹⁴ K s⁻¹.¹⁶ This pioneering work removed long-standing theoretical barriers to the study of non-equilibrium solidification in liquid metals and the atomic structure of metallic glasses, thus opening the door to research utilizing ultra-fast cooling rates. As a result, it has become almost inconceivable that anyone would question the experimental feasibility of cooling rates used in computational simulations, ranging from 10⁹ K s⁻¹ to 10¹⁴ K s⁻¹. Following this breakthrough, a series of MD simulations have been conducted to explore the quenching behavior of these refractory Ta-based metals by Kbirou et al.,¹⁷ as well as nanoscale heterogeneity in Al₇₀Ni₁₅Co₁₅ metallic glass,¹⁸ atomic packing and fractal behavior in Al–Co metallic glasses,¹⁹ and the dynamics of Cu–Zr metallic glass.²⁰ Additionally, in 2022, Li et al. reported a data-driven discovery of a universal indicator for metallic glass-forming ability.²¹ These seminal studies have provided invaluable insights into the atomic-level processes governing rapid cooling and metallic glass formation, establishing a solid foundation for ongoing and future work in this field.

Therefore, we investigate and highlight the evolution of the first peak of the RDF during the rapid quenching of tantalum at five different cooling rates (γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹, and γ₅ = 1 × 10¹⁴ K s⁻¹). Compared to the study of Al, Ni, Cu, Ag, Au, and In conducted a decade ago,¹¹ we not only provide computational evidence in this paper that the anomalous expansion of the first peak in the RDF is an objective physical phenomenon, but also propose that this anomaly arises from the combined effects of the geometric structure of newly formed atomic clusters and the electronic structure of their internal nearest-neighbor atoms. Our study offers solid evidence for a deeper understanding of the causes of the anomalous shift in the first maximum peak of the RDF with decreasing temperature in metallic melts and solids, and helps clarify the long-standing controversy surrounding this phenomenon.

2. Computational methods

2.1. MD simulation

In our molecular dynamics (MD) simulations, we employed the large-scale atomic/molecular massively parallel simulator (LAMMPS)²² to investigate the energy relaxation and cooling behavior of tantalum (Ta). A cubic simulation box containing 11 664 Ta atoms was constructed with three-dimensional periodic boundary conditions to eliminate finite-size effects. The simulation followed a melt–quench protocol comprising two main stages.²³ First, the monatomic Ta system was melted at 4000 K—well above its experimental melting temperature (∼3296 K)¹⁶—and equilibrated until the energy no longer varied with simulation time, ensuring that the system had reached a stable high-temperature liquid state. Subsequently, the identical initial atomic configuration was cooled from 4000 K to 300 K at five different cooling rates (γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹, γ₅ = 1 × 10¹⁴ K s⁻¹), respectively, to investigate the influence of cooling rates on the resulting atomic structure. For simulation control, the system was maintained in the NPT ensemble,²⁴ with temperature and pressure regulated via a Nosé–Hoover thermostat and barostat.²⁴ A time step of 1 fs was used, and neighbor lists were updated every 5 steps to optimize computational efficiency. Atomic interactions were modeled using the embedded atom method (EAM) potential developed by Sheng et al.,¹⁶ which was parameterized based on density functional theory (DFT) calculations^25,26 to accurately reproduce the potential energy surface of Ta. Thermodynamic properties, including temperature, total energy, pressure, and volume, were recorded at regular intervals, while atomic positions and coordinates were continuously output for subsequent structural analysis. Due to the different cooling rates, the total simulation time varied accordingly. For γ₁ = 1 × 10¹⁰ K s⁻¹, the total simulation consisted of 370 000 000 steps, corresponding to 370 ns. For γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹ and γ₅ = 1 × 10¹⁴ K s⁻¹, the total number of simulation steps and times were 37 000 000 (37 ns), 3 700 000 (3.7 ns), 370 000 (0.37 ns), and 37 000 (0.037 ns), respectively. The collected data provide valuable insights into the structural evolution and thermodynamic stability of Ta under different cooling conditions. For the RDF and Voronoi analyses, due to computational and workload constraints, only representative temperature points were selected for plotting, corresponding to data from approximately 700 configurations. The relative error in the peak position of the RDF is approximately 0.0001 Å, while the error in the numerical integration of the area between the RDF curve and the x-axis is about 0.0001 Ų.

2.2. DFT calculation

We applied the DMOL³ software^27,28 to optimize the binding energy of the ten most frequent Voronoi polyhedron. The total energy calculation utilized the Perdew–Burke–Ernzerhof (PBE)²⁹ functional based on the generalized gradient approximation (GGA)³⁰ to represent electronic exchange–correlation functions. The electronic wave functions of Ta atoms were characterized using a double-numerical basis set with d-polarization functions (DN).²⁷ In the process of atomic position relaxation in Ta clusters, the convergence thresholds for geometry optimization were set as follows: energy change of 1 × 10⁻⁴ Ha, maximum force of 0.02 Ha Å⁻¹, and maximum displacement of 0.05 Å. The total energy and electronic structure calculations were performed after geometry optimization, with a self-consistent field (SCF) tolerance of 1 × 10⁻⁴ Ha. Fig. 8 shows the results of these calculations obtained using DMOL³ (ref. 27 and 28) for various sizes of Ta clusters. The electron localization functions (ELFs)³¹ for the Ta atomic structure (comprising approximately 60–80 atoms), extracted from MD simulations for DFT calculations at 3500 K and 300 K, are shown in Fig. 12. These selected atoms were then placed in a vacuum box of size 15 Å × 15 Å × 15 Å. Subsequently, different planes were randomly sliced in various directions and labeled as Slice-1 and Slice-2. In order to calculate the Mulliken bond overlap population and the average bond length of the nearest atoms in Ta solid as a function of temperature, as shown in Fig. 13, we applied Cambridge serial total energy package (CASTEP) software.³² Ultrasoft pseudopotentials were represented in reciprocal space with an exchange–correction function of PBE²⁹ under the GGA framework³⁰ for all Ta atoms involved in the calculation, as shown in Fig. 13. The cutoff energy for atomic wave functions (PWs) was set to 260 eV. The SCF convergence threshold was set to 2.0 × 10⁻⁶ eV per atom, with a maximum of 1000 SCF iterations. The total energy and electronic structure calculations followed a single-point energy calculation with an SCF tolerance of 2.0 × 10⁻⁶ eV per atom.

2.3. Visualization and analysis tool

In this work, the post-processing of MD simulations,^33–36 such as the visualization of atomic structures, was partially carried out using open visualization tool (OVITO) software.³³ OVITO is a 3D visualization software designed for post-processing atomistic data obtained from MD simulations.^33–36 It was developed by Alexander Stukowski at Darmstadt University of Technology, Germany.³³

2.4. Parallel tempering (PT)

Parallel tempering (PT)²² is a powerful technique for improving ergodic sampling, particularly in systems with high energy barriers, such as glass-forming liquids at low temperatures. Employing PT would fundamentally alter the nature of the rapid solidification process by allowing the system to cross energy barriers that it physically would not overcome in a real quenching experiment. To verify this, we further performed PT simulations using the LAMMPS software²² under a multi-replica setup. To enhance ergodicity, we implemented a parallel temperature exchange scheme with the following parameters: Eight temperature replicas were set at 4000 K (S0), 3500 K (S1), 3000 K (S2), 2500 K (S3), 2000 K (S4), 1500 K (S5), 1000 K (S6), and 300 K (S7), with a time step of 0.001 ps and a total simulation duration of 37 ns (37 000 000 steps in total). Temperature swaps were attempted every 5000 timesteps to ensure frequent energy exchanges between different replicas, thereby improving global potential energy surface (PES) exploration. The trajectories of each replica were also tracked to prevent the system from being trapped in specific energy wells for extended periods, further enhancing global ergodicity. As shown in Fig. 1(a) and (b), temperature and energy exchanges successfully occurred between different replicas as the simulation progressed. By monitoring the structural evolution across all eight temperature replicas and performing Voronoi polyhedral (VP) analysis,³⁷ we found that the final geometric configuration converged to the perfect global ground-state BCC (0,6,0,8) structure. For instance, Fig. 1(c) and (d) depict the VP distributions for the S7 sample at 4000 K and 300 K, respectively, highlighting the ten most frequent structures. At 4000 K (Fig. 1(c)), the characteristic BCC index (0,6,0,8) was absent. However, at 300 K (Fig. 1(d)), this index (0,6,0,8) accounted for 92.9% of the total VP population, clearly demonstrating that the system formed a global BCC ground-state structure at ambient temperature. A detailed comparison of the RDF³⁴ curves for sample S7 at 4000 K and 300 K (Fig. 1(e) and (f)) shows that the position of the first RDF peak remains unchanged at both temperatures. Notably, no anomalous peak broadening was observed at 300 K in Fig. 1(e) and (f). The consistent position of the first RDF peaks in both BCC Ta and high-temperature liquid Ta further supports our findings, suggesting that the formation of the global BCC ground state does not lead to peak broadening upon cooling. Instead, the anomalous shift in the first maximum peak of RDF with decreasing temperature are attributed to metastable geometric structures formed by interlocking VP (0,0,12,2) and (0,0,12,0) in Fig. 11. These findings demonstrate that PT is effective in enhancing ergodic sampling and locating global minimum-energy configurations by enabling the system to overcome high energy barriers. However, this also leads to the neglect of many competitive metastable structures. As a result, while PT is valuable for equilibrium studies, it is not suitable for investigating non-equilibrium solidification processes and metastable metallic glasses, as addressed in our study.

Fig. 1 Evolution of (a) temperature and (b) potential energy for eight different replicas during the parallel tempering process. Population of the ten most frequent Voronoi polyhedron for S7 at (c) 4000 K and (d) 300 K. (e) Total RDF curves for the simulated BCC and liquid Ta at 300 K and 4000 K, respectively. (f) Zoomed-in views of the first RDF peak for each corresponding curve in (e).

2.5. Physical rationale for metastable states in rapid solidification

Given that our work employs a rapid solidification technique far from thermodynamic equilibrium, it is essential to elaborate on the physical reasoning underlying our methodology, particularly for readers who may be unfamiliar with the non-equilibrium aspects of glass formation. From the standpoint of metallic glass formation, these materials form when high-temperature metallic melts are cooled so rapidly that crystallization is effectively suppressed, locking atoms into a disordered, amorphous structure. What distinguishes a glass from a crystal is not composition, but time: the system simply does not have enough time to organize itself. As the liquid cools below its liquidus temperature, atomic rearrangement toward order begins. If cooling is slow, atoms can diffuse, nucleate, and grow into crystals.^16,34 However, when the cooling rate is sufficiently high—exceeding the characteristic timescales of atomic diffusion and structural relaxation, the system becomes dynamically arrested before crystalline nuclei can form and grow. In this regime, vitrification becomes a competition between the external cooling rate and the intrinsic kinetics of atomic rearrangement, rendering glass formation a fundamentally non-equilibrium, time-dependent process. Moreover, this non-equilibrium transition can be more rigorously understood within the framework of the potential energy landscape (PEL). At elevated temperatures, the system has sufficient thermal energy to surmount energy barriers between basins, thereby exploring the configuration space extensively and exhibiting ergodic behavior.³⁴ As the temperature decreases rapidly, the system becomes confined to increasingly localized regions of the PEL, namely, deep metastable basins or megabasins separated by high barriers—leading to the breakdown of ergodicity.³⁴ In our simulations, we explicitly exploit this kinetic trapping mechanism by employing ultrafast quenching protocols that emulate the physical conditions of metallic glass formation. The aim is not to reach the global energy minimum—as sought in ergodic ensemble techniques like PT, but to faithfully reproduce the metastable states that emerge under realistic rapid solidification conditions. Non-ergodic sampling is not a methodological artifact but a physically grounded necessity for modeling the dynamic pathways of glass formation.¹⁶

Despite their strengths, current simulations are limited by spatial and temporal constraints, often employing system sizes (i.e., the number of atoms modeled) far smaller than those relevant to real-world solidification. In particular, to reduce computational cost, simulations often use cooling rates far exceeding those in experimental conditions. This, combined with the nonergodic nature of the process, means that simulations typically explore only a narrow region of phase space. Therefore, these methods should be viewed not as definitive predictors, but as evolving frameworks grounded in nonequilibrium statistical mechanics. Continued advancements in high-performance computing and artificial intelligence (e.g., machine learning) are expected to substantially enhance their predictive accuracy and applicability.

3. Results and discussions

3.1. Radial distribution functions (RDFs)

For RDFs, defined as , ρ is the number density of atoms in the system of N atoms. r⃑_ij is the interatomic distance between i and j atoms.³⁴ For crystalline materials, the RDF shows a series of sharp peaks as r increases. For amorphous materials, however, the intensity of the peaks in the RDF gradually decreases and weakens as r increases, and the second peak also shows signs of splitting at the onset of the glass transition during rapid solidification. This behavior is used as an indicator of the formation of amorphous materials.^34–36 The atomic structure information is embedded in the peak position, peak width, and relative intensity of the RDF. Thus, we can observe that not only do energy, volume, and pressure development strongly depend on the cooling rate (Fig. S1, ESI†), but also the variation in the RDF curves of Ta_11 664 is affected by different cooling rates (Fig. S2, ESI†). To examine how temperature affects the shift of the first maximum peak in the RDFs, Fig. 2(a–e) show the total RDF curves of the simulated Ta_11 664 system at selected temperatures under five different cooling rates (γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ =1 × 10¹³ K s⁻¹, and γ₅ = 1 × 10¹⁵ K s⁻¹), respectively. The local zoomed-in views of the first peak of each corresponding curve are shown in Fig. 2(a′–e′). For γ₁ = 1 × 10¹⁰ K s⁻¹, one can find that the position of the first peak, r, gradually increases from 2.8875 Å at 4000 K to 2.9625 Å at 300 K, as shown in Fig. 2(a′). The first peak in the RDFs shifts to a higher r and intensifies as the temperature decreases. This behavior contrasts with the typical expansion of metals when heated and contraction when cooled. In other words, this indicates that the average bond length between nearest-neighbor atoms in molten Ta at high temperatures is shorter than in the crystalline material at low temperatures. Similar anomalous behavior, where the first peak maximum shifts to higher r values for γ₂ = 1 × 10¹¹ K s⁻¹, is also observed in Fig. 2(b′), but this is not observed in Fig. 2(c′–e′). The vertical dashed lines in Fig. 2(c′–e′) further reinforce the observation that the position of the first peak in the RDF is independent of temperature for the cooling rates γ₃, γ₄ and γ₅. Additionally, the RDF at 4000 K from our study closely aligns with the data in ref. 16, supporting the reliability of our results.

Fig. 2 (a–e) The total RDF curves of the simulated Ta_11 664 system at selected temperatures under five different cooling rates (γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹ and γ₅ = 1 × 10¹⁵ K s⁻¹), respectively, and the figures in (a′–e′) are local zoomed-in views of the first peak of each corresponding curve in (a–e). The RDF data, represented by a bright green circle, correspond to a temperature of 4000 K, as reported in ref. 16.

Because high-resolution peaks of the RDF usually provide key information for unraveling the structures of liquids, to further demonstrate that the anomalous broadening of the first peak observed in Fig. 2 is indeed a physical phenomenon rather than a data statistical error, we further present in Fig. 3 the temperature dependence of the positions of the first and second peaks of the RDF at five different cooling rates. For cooling rates of γ₁ = 1 × 10¹⁰ K s⁻¹ and γ₂ = 1 × 10¹¹ K s⁻¹, as shown in Fig. 3(a) and (b), when the temperature decreases continuously from 4000 K to 300 K, the first RDF peak position shifts abruptly to a higher r, from 2.8875 Å to 2.9625 Å at 2010 K and 1850 K, respectively. However, at the same cooling rates, the second peak position undergoes a sudden shift from a higher r region at high temperature to a lower r region at the crystallization temperature. In fact, not only does the second peak follow the expected thermal expansion and contraction, but the third, fourth, and even higher-order peaks also exhibit this behavior (for clarity, the changes in the peaks at larger r are not shown in Fig. 3). This is in stark contrast to the anomalous expansion behavior observed for the first peak at temperatures below the crystallization point. Moreover, as shown in Fig. 3(b), the first peak remains at 2.9625 Å for most temperatures below 1850 K. However, at 700 K, 900 K, and 1100 K, it shifts to around 2.8875 Å, exhibiting significant fluctuations and anomalies. This is primarily attributed to temperature variations and system fluctuations during non-equilibrium solidification. For higher cooling rates of γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹, and γ₅ = 1 × 10¹⁴ K s⁻¹, the first peak position remains largely unchanged with decreasing temperature, while the second peak gradually shifts from a higher r region to a lower r region as the temperature decreases. The blue dashed lines in Fig. 3(c)–(e) highlight the clear trend, as indicated by the pink arrows. Furthermore, the position of the first peak (r_max) of the RDF as a function of five different cooling rates at 300 K and 4000 K is depicted in Fig. 3(f). From Fig. 3(f), it is observed that r_max for cooling rates of γ₁ = 1 × 10¹⁰ K s⁻¹ and γ₂ = 1 × 10¹¹ K s⁻¹ at 300 K is higher compared to that at 4000 K. Conversely, for cooling rates of γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹ and γ₅ = 1 × 10¹⁵ K s⁻¹, there is no discernible relationship between r_max and temperature. Clearly, both Fig. 2 and 3 demonstrate that at cooling rates of γ₁ = 1 × 10¹⁰ K s⁻¹ and γ₂ = 1 × 10¹¹ K s⁻¹, the first RDF peak shows anomalous expansion behavior with decreasing temperature, which contradicts established thermodynamic principles. There must be a deeper underlying cause for this phenomenon.

Fig. 3 Temperature dependence of the position of the first and second peak of RDF at five different cooling rates. (a) γ₁ = 1 × 10¹⁰ K s⁻¹; (b) γ₂ = 1 × 10¹¹ K s⁻¹; (c) γ₃ = 1 × 10¹² K s⁻¹; (d) γ₄ = 1 × 10¹³ K s⁻¹; (e) γ₅ = 1 × 10¹⁴ K s⁻¹; (f) the first peak position r_max of RDF as a function of five different cooling rates at 300 K and 4000 K, respectively.

To understand the underlying mechanisms behind this phenomenon, we conduct a more refined analysis of the enclosed area under the first peak in Fig. 4. Fig. 4(a) shows the first peak and subsequent peaks up to the fifth peak of the RDF, with different filling patterns for distinction, as a schematic diagram. The area (integral) under the RDF curve represents the atomic number density or probability in that region. At the same time, the evolution of the integration of the first peak of the RDF as a function of temperature at five different cooling rates (γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹, and γ₅ = 1 × 10¹⁵ K s⁻¹) is plotted in Fig. 4(b–f), respectively. From Fig. 4(b) and (c), we can see that as the temperature decreases, for γ₁ = 1 × 10¹⁰ K s⁻¹ and γ₂ = 1 × 10¹¹ K s⁻¹, the integrated value of the first RDF peak sharply decreases at the crystallization temperature. At the same time, their respective crystallization (or phase transition) temperatures are nearly identical to those shown in Fig. S1 (ESI†), based on energy scales. For γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹, and γ₅ = 1 × 10¹⁵ K s⁻¹, although the integrated values of the first peak in the RDF show a monotonically decreasing trend with decreasing temperature, the values at 300 K are still relatively large compared to γ₁ and γ₂. This indicates that, during solidification at these three cooling rates, the slight changes in the number of atoms in the first coordination shell are not sufficient to cause a reverse shift in the position of the RDF first peak's maximum. Therefore, whether it is the process of vitrification or crystallization, if the cutoff distance r is the same, both processes will result in a constant shell volume. Under this premise, while the number of atoms in the nearest-neighbor region decreases sufficiently with decreasing temperature, this will lead to an increase in the interatomic distances within the nearest-neighbor shell. Naturally, this change in interatomic distances will cause the maximum peak position in the RDF to shift towards larger r values. Next, we will focus on the cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹ as an example and analyze the underlying physical mechanism behind the anomalous expansion behavior of the position of the maximum RDF peak as the temperature decreases.

Fig. 4 (a) Schematic diagram of a series of peaks of radial distribution function (RDF) of the simulated Ta_11 664 system at 4000 K under a cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹. (b–f) The temperature dependence of the integral value of the first peak area of RDF functions for γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹ and γ₅ = 1 × 10¹⁴ K s⁻¹, respectively.

3.2. Coordination number and Voronoi polyhedron

How to characterize the local atomic configurations in a rapidly solidified system remains a significant scientific challenge.³⁴ To identify the atomic clusters in a rapidly solidified system with a cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹, we decided to employ the Voronoi analysis method to study the distribution of coordination numbers and polyhedral atomic clusters in detail, with the aim of accurately capturing the key structural information that causes the anomalous expansion of the first peak in the RDF. The Voronoi polyhedron method (VPM) is a typical method introduced by Finney³⁷ that does not require predefined templates. In VPM, space is divided into numerous polyhedra, each centered on an atom, and their local structures can be characterized by various Voronoi indices (n₃, n₄, n₅, n₆, …, n_i, …), where n_i denotes the number of i-edge polygons.³⁸ This method is mathematically rigorous, but there is some ambiguity due to the high sensitivity of Voronoi polyhedra (VPs) to structural distortions caused by thermal fluctuations, atomic perturbations,³⁹ and even the influence of the atomic sizes themselves. As a result, it is difficult to establish a universal criterion for removing ultra-short edges and tiny surfaces.⁴⁰ To establish a reasonable critical threshold, we focused on a cooling process with γ₁ = 1 × 10¹⁰ K s⁻¹, as shown in Fig. 5. Initially, we analyzed the distribution of areas and edge lengths of all original VPs at temperatures of 4000 K and 300 K, without any parameter approximation. At 4000 K (Fig. 5(a)), there was a distinct valley around 0.65 Ų between tiny surfaces and actual polyhedral areas, which served as the critical value for removing tiny surfaces. At 300 K (Fig. 5(b)), there was a broader, flatter valley ranging approximately from 0.13 Ų to 1.0 Ų. Thus, selecting thresholds within this range was deemed acceptable for analyzing VPs at this temperature. Considering both temperature scenarios (4000 K and 300 K), we decided to exclude surfaces smaller than 0.65 Ų during VP analysis, and anything equal to or larger than this value was considered an actual surface during data collection. Regarding edge distributions in both scenarios (4000 K and 300 K), no clear valleys were observed in Fig. 5(c) and (d). To prevent disregarding genuine small edges while analyzing VPs, we referred to Ovito's user manual, where a default edge length threshold of 0.1 Å was suggested: including edges greater than or equal to this value while excluding those smaller than it.

Fig. 5 Distribution of (a) 4000 K and (b) 300 K areas, along with (c) 4000 K and (d) 300 K edge lengths of Voronoi polyhedra during the nonequilibrium solidification process with a cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹.

After determining the threshold for removing ultra-short edges and tiny surfaces, we proceeded to analyze the distribution of coordination numbers and the number of Voronoi polyhedron, as shown in Fig. 6 and 7. Fig. 6(a) reveals that at a high temperature (4000 K), the distribution of coordination numbers is broad, ranging from 9 to 16, while coordination numbers of 8 and 17 are relatively infrequent, and 18 is one of the few rare points. Conversely, at lower temperatures (300 K) as shown in Fig. 6(b), coordination numbers predominantly cluster around 12, 13, 14, and 15, with a significant reduction in clusters having coordination numbers of 16 or beyond. Similarly, the number of clusters with coordination numbers less than 11 also decreases significantly, even disappearing. The marked change in coordination numbers reflects a significant transition in the microstructure. To further explore the detailed structural information of the microclusters, Fig. 7(a) and (b) depict the distribution of the top ten Voronoi polyhedral types by frequency for liquid Ta at 4000 K and 300 K, respectively.

Fig. 6 Population of the coordination number of those particles at (a) 4000 K and (b) 300 K, under the cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹, respectively.

Fig. 7 Population of the ten most frequent Voronoi polyhedron clusters at (a) 4000 K and (b) 300 K, under the cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹, respectively.

Specifically, Fig. 7(a) shows that at 4000 K, the top ten Voronoi polyhedral are (0,3,6,4), (0,2,8,4), (0,1,10,2), (0,2,8,2), (0,3,6,5), (1,2,7,3), (0,3,6,3), (1,3,5,4), (0,4,4,6), and (0,2,8,3), with their respective percentages being 2.6%, 2.0%, 1.7%, 1.5%, 1.5%, 1.2%, 1.2%, 1.2%, and 1.0%. Although these ten types of Voronoi polyhedral are the most frequent, they still account for only a small proportion of the total number of atoms, making up just 13.9% of the system's total atoms. This indicates that, at high temperatures, the majority of atomic clusters remain disordered, which is a characteristic feature of liquid metals. In contrast, Fig. 7(b) shows a marked difference in the distribution of Voronoi polyhedral at 300 K. The top ten Voronoi polyhedral at this temperature are (0,0,12,2), (0,0,12,0), (0,0,12,3), (0,0,10,2), (0,3,6,4), (0,2,8,4), (0,1,10,4), (0,1,10,3), (0,2,8,5), and (0,3,6,5). Notably, the combined percentage of (0,0,12,2), (0,0,12,0), and (0,0,12,3) accounts for 79.4% of the total number of atoms. In particular, the cage-like Voronoi polyhedral (0,0,12,2), formed by the interpenetration of two hexagonal bipyramids sharing one atom, and the icosahedron (0,0,12,0) account for significant proportions of 41.8% and 27.3%, respectively. Our previous research²⁶ has revealed that during non-equilibrium solidification, the Voronoi polyhedral (0,0,12,2) and the icosahedron (0,0,12,0) can interpenetrate to form a new Ta₂₆-C_2v atomic cluster²⁶ with both five-fold and six-fold symmetry. Extensive studies have demonstrated that even in non-equilibrium solidification, the distribution of atomic clusters in liquid metals tends to adhere to the principle of minimum energy. In this study, we compare cluster energies through the calculation of binding energy, which is determined using the following formula,

E_b = [E_total(Ta_n) − nE(Ta)]/n (1)

where E_total(Ta_n) and E(Ta) are the total energies of the neutral Ta_n clusters and single gaseous Ta atom, respectively. n represents the number of Ta atoms. The binding energy can be utilized to characterize the relative stability of metal clusters. Namely, the lower the E_b is, the higher the stability of clusters is. Fig. 8(a) and (b) illustrate the average binding energy per atom²⁶ and the evolution of the Voronoi polyhedron numbers for the ten most frequent types at temperatures of 4000 K and 300 K, respectively. Clearly, regardless of temperature, Voronoi polyhedra with larger numbers generally exhibit lower binding energies, while the binding energy gradually increases as the number of polyhedra decreases from left to right. This suggests that during non-equilibrium solidification, Voronoi polyhedra with lower energies are more likely to form and thus dominate in number. Additionally, it implies that the anomalous expansion of the first peak in the RDF of the metal Ta, solidified at a cooling rate of γ₁, is closely related to the presence of (0,0,12,2) and (0,0,12,0). The subsequent analysis of both the total and local RDFs in Fig. 9 further supports this view. Fig. 9 illustrates the RDFs of Ta metal solidified at 300 K with a cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹, alongside the RDFs of the principal atomic clusters obtained through Voronoi analysis. Fig. 9(b) provides an enlarged view of Fig. 9(a). From the RDF peak intensities, it is evident that the primary contributors to the first RDF peak at 300 K are the (0,0,12,2) and (0,0,12,0). Although clusters with Voronoi indices (0,0,12,3), (0,0,10,2), (0,3,6,4), and (0,2,8,4) also contribute to the overall RDF, their contributions are relatively minor. In terms of the RDF's first peak position, the Voronoi polyhedral (0,0,12,2) has its peak prominently situated in the higher r region, with the main peak between 3.25 Å and 3.50 Å, and the secondary peak between 2.50 Å and 2.84 Å. In contrast, the Voronoi polyhedral (0,0,12,0), a typical representative of five-fold symmetry, has its main peak between 2.65 Å and 3.00 Å, and its position is slightly higher than the secondary peak of the (0,0,12,2) in Fig. 9(b). These observations clearly indicate that (0,0,12,2) and (0,0,12,0) are the primary factors affecting the intensity and position of the RDF's first peak during cooling at γ₁ = 1 × 10¹⁰ K s⁻¹.

Fig. 8 The binding energy of the ten most frequent Voronoi polyhedron at (a) 4000 K and (b) 300 K, under the cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹, respectively.

Fig. 9 (a) Radial distribution function (RDF) curves of simulated Ta_11 664 system and the six most frequent Voronoi polyhedron at 300 K under the cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹, respectively. (b) is an enlarged partial view of (a).

We are left with two key questions. First, as the temperature decreases, which specific atoms within the Voronoi polyhedron (0,0,12,2) and (0,0,12,0) contribute to the shift of the RDF's first peak to higher r values? Second, since the (0,0,12,2) and (0,0,12,0) in the system are interpenetrated rather than isolated, how do they interpenetrate? To address the first question, we focus on Fig. 10, which provides insights into the specific atomic contributions to the anomalous expansion of the RDF's first peak by examining the Voronoi polyhedral (0,0,12,2) and (0,0,12,0) extracted from MD simulations. Due to the high symmetry of the standard (0,0,12,2) and (0,0,12,0), which feature numerous equivalent atomic positions and equal bond lengths, Fig. 10 only illustrates selected central atoms and shell atoms, as well as some bond lengths and atom indices between shell atoms. For the (0,0,12,2), the theoretical bond length between the central atom O and the vertex atoms B and G is equivalent. The bond length L_OB is 2.745 Å. The bond lengths between the central atom and the 12 surrounding shell atoms are approximately 3.33 Å, which is about 0.60 Å longer than L_OB(L_OG). Additionally, the bond lengths between shell atoms within the (0,0,12,2) range from 3.30 Å to 3.50 Å. This bond length distribution is consistent with the RDF bond length distribution for (0,0,12,2) shown in Fig. 9. For the standard icosahedral (0,0,12,0), which possesses I_h point group symmetry, it represents the highest symmetry structure in nature besides the sphere. The bond lengths between the central atom and shell atoms are uniform, and the distance between the central atom and the shell is 5% shorter than the distance between shell atoms.²⁵ For example, in Fig. 10, the bond length L_OE between the central atom and shell atoms is 2.796 Å, while the bond length L_AB between shell atoms is 2.940 Å, with a difference of 0.144 Å, which is 5.1% of L_OA, aligning with theoretical expectations.²⁵ Thus, from a geometric perspective, the anomalous expansion of the RDF's first peak is primarily due to two factors: the longer bond lengths between shell atoms and between the central atom and most shell atoms in the (0,0,12,2), as well as the bond lengths of shell atoms in (0,0,12,0) icosahedron. These results are a consequence of the mathematical characteristics inherent to the high-symmetry cage-like structures of (0,0,12,2) and (0,0,12,0).

Fig. 10 Schematic diagram of the local atomic configuration of the (0,0,12,2) and (0,0,12,0) from MD simulation.

For the second question regarding how the (0,0,12,2) and (0,0,12,0) clusters interpenetrate in the actual system, we first examine the formation of (0,0,12,2) in Fig. 11. Fig. 11 illustrates that a tetrahedral structure (Ta₄-D_2d) composed of 4 atoms, forms the basis, which then combines with 6 tetrahedra to create a dodecahedron containing 8 atoms. Furthermore, two dodecahedrons share a vertex atom to form a trisoctahedron (0,0,12,2) with 15 atoms. On this basis, adding a 7-atom hexagonal pyramid (light green atoms) to the apex results in Ta₂₂-D_6h, which can also be viewed as two isolated (0,0,12,2) nested together by sharing the hexagonal bipyramids. If, following the direction of the arrows in Fig. 11, we add 4 more atoms (dark red atoms) to the middle outer region of the Ta₂₂-D_6h cluster, we obtain Ta₂₆-C_2v, which can also be considered as an isolated standard icosahedron (0,0,12,0) nesting with the Ta₁₃-C_2v cluster on its left.

Fig. 11 The evolutionary schematic diagram of the local atomic configuration of the (0,0,12,2) and Ta₂₆-C_2v structure from the Ta₄ to Ta₂₆ cluster with point group symmetry.

In the Ta₂₆-C_2v geometric configuration highlighted by the dashed circle in Fig. 11, both (0,0,12,2) and (0,0,12,0) are present. Thus, during non-equilibrium solidification, the Voronoi polyhedra (0,0,12,2) and (0,0,12,0) actually exist in the form of Ta₂₆-C_2v. The low binding energy of Ta₂₆-C_2v,²⁶ as shown in Fig. 8(b), provides a reasonable theoretical explanation for its substantial presence.

3.3. Correlation between electronic structure and bond lengths

In 1947, Linus Pauling⁴¹ reported that the shortening of bonds in metals is often due to an increase in bond order. In fact, our MD simulations have found that the average bond length of nearest-neighbor atoms at high temperatures does indeed shorten during rapid solidification at cooling rates of γ₁ = 1 × 10¹⁰ K s⁻¹ and γ₂ = 1 × 10¹¹ K s⁻¹, while it elongates at low temperatures. We believe the origin of this phenomenon is not solely due to a decrease in the coordination number of nearest neighbors at high temperatures, as suggested in the existing literature.^11–15 Based on Linus Pauling's theory,⁴¹ we propose that the electronic interactions between different atoms likely have a more significant effect on changes in bond lengths. Therefore, to illustrate the bonding features, we performed a qualitative and intuitive comparison by calculating the electron localization function (ELF)^31,42,43 for Ta, as shown in Fig. 12. The ELF is a measure of relative electron localization in a system, and higher ELF values are usually found in regions with a high likelihood of electron pairing, such as cores, bonds, and lone pairs.^31,42,43 However, due to the complex atomic structure of Ta metallic solids and the large number of atoms in MD simulations, it is not feasible to directly use rapidly solidified structures with tens of thousands of atoms for DFT calculations. Therefore, to facilitate DFT calculations, a random atomic structure is extracted from the MD simulation (Fig. 12(a)) for ELF analysis. This chosen structure is directly used for DFT computations without any modifications or geometric optimization, leading to a notable reduction in the computational cost of DFT (Fig. 12(b)). At 3500 K, Fig. 12(b) shows clear electron localization in randomly selected bonds of A–B, C–D, D–E, and E–F on Slice-1, as well as G–H, I–J, and L–K bonds exhibiting distinct electron localization on the Slice-2 plane. In contrast, compared to the ELF at 3500 K, the ELF of M–N, P–M, P–N, P–R, R–S, S–Q, Q–O in Slice-1 and W–T, T–U, U–V, V–W, U–X, U–Y, and X–Y in Slice-2 shows weaker electron localization at 300 K. Studying the ELF can reveal whether the bond between two atoms is covalent, metallic, or ionic, and how many electrons each atom contributes to each bond. Thus, the ELF values at high temperature (3500 K) are significantly higher than those at low temperature (300 K) in Fig. 12, which confirms our physical intuition that strong electronic interactions lead to bond contractions at high temperatures.

Fig. 12 (a) A schematic of the atomic structure extracted from MD simulations with DFT calculations, and (b) electron localization functions (ELFs) for different slices of the aforementioned (a) atomic structure at 3500 K and 300 K, respectively.

Moreover, in order to further analyze the correlation between electronic interactions and bond lengths of Ta metal from a quantitative perspective, Mulliken population is used to calculate the magnitude of covalent bonds.^44,45 The bond overlap population Q_i–j between i and j atoms is defined as follows:^44,45

where P_μv(κ) and S_μv(κ) are the density matrix and the overlap matrix,^44,45 respectively. ω_K is the weight associated with the calculated K-points in the Brillouin zone. Usually, the magnitude of Q_i–j can be used as an approximate measure of covalent bonding strength.^44,45 In this section, the Mulliken bond overlap population and bond length of the nearest atoms as a function of temperature by DFT calculations, are plotted in Fig. 13(a) and (b), respectively. In Fig. 13(a), as the temperature decreases, there is a notable decrease in Q_i–j. Meanwhile, there is a notable increase in bond lengths in Fig. 13(b). Although the trend of change between bond order and temperature is completely opposite to that of bond length and temperature, the physical essence reflected by their correlations is completely consistent and logically self-consistent. Specifically, the change of Mulliken bond overlap reveals a stronger electronic interaction between the nearest atoms at high temperature compared to lower temperature. This highlights that at high temperatures, the nearest atoms of liquid Ta have a stronger electronic interaction, in contrast, a weaker electronic interaction exists in the Ta solid. Obviously, temperature-induced bond contractions in the Ta melt indeed can be attributed to the stronger electronic interaction among the nearest atoms at high temperatures. We have provided ample evidence from first principles calculations in Fig. 13, which strongly supports this point.

Fig. 13 (a) The Mulliken bond overlap population and (b) bond length of the nearest atoms as a function of temperature by DFT calculations, respectively.

4. Conclusions

In this paper, we investigate the evolution of the first RDF peak of rapidly quenched tantalum using molecular dynamics (MD) simulations and density functional theory (DFT) calculations at five different cooling rates (γ₁ = 1 × 10¹⁰ K s⁻¹, γ₂ = 1 × 10¹¹ K s⁻¹, γ₃ = 1 × 10¹² K s⁻¹, γ₄ = 1 × 10¹³ K s⁻¹, and γ₅ = 1 × 10¹⁴ K s⁻¹). For γ₁ and γ₂, MD simulations show that the first peak in the RDF shifts to a higher r and intensifies as the temperature decreases. This behavior indicates that the average bond length between nearest-neighbor atoms in molten Ta at high temperatures is shorter than that in the crystalline material at low temperatures. We focused on the rapid solidification process with a cooling rate of γ₁ = 1 × 10¹⁰ K s⁻¹ as an example for a systematic analysis, and found that this anomalous behavior is primarily due to the system forming a significant number of cage-like icositetrahedral Voronoi polyhedra (0,0,12,2) and standard icosahedral Voronoi polyhedra (0,0,12,0) at 300 K. Moreover, these Voronoi polyhedra nest into each other, forming a stable Ta₂₆-C_2v atomic configuration with longer bond lengths and lower binding energies compared to their high-temperature counterparts. This unique geometric arrangement fundamentally contributes to the anomalous expansion of the RDF's first peak. On the other hand, at higher temperatures, the electronic localization functions (ELFs) and the Millikan bond overlap populations are significantly increased, leading to stronger electronic interactions and a denser arrangement of nearest-neighbor atoms with shorter bond lengths. Consequently, the temperature-induced anomalous shift in the first peak of the RDF appears to originate from the unique geometric arrangement of (0,0,12,2) and (0,0,12,0) at lower temperatures, alongside the enhanced electronic interactions among nearest-neighbor atoms at higher temperatures. Our findings provide strong evidence supporting this conclusion.

Data availability

The data supporting this article have been included in this article, and other data are available upon request from the authors.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The authors acknowledge funding support from the National Natural Science Foundation of China (Grants No. 52263025, 51871096, 91961204, and 11974134), the Scientific Research Project of Education Department of Jiangxi Province (Grants No. GJJ2202021 and GJJ2202011), the Jiangxi Provincial Natural Science Foundation (Grants No. 20202BAB204004 and 20171BAB216001), the Qinglan Scholars Program and 2023 Teaching Reform Research Project (Grants No. NSJG-23-08) of Nanchang Normal University, and the Jilin Province Outstanding Young Talents project (Grants No. 20190103040JH). Yuanqi Jiang also appreciates the support from the Young and Middle-aged Academic (Professional) Leaders Program of Higher Education Institutions in Jiangxi Province (2024).

References

1. Y. Song, N. Shi, S. Deng, X. Xing and J. Chen, Negative thermal expansion in magnetic materials, Prog. Mater. Sci., 2021, 121, 100835.
2. P. Mohn, A century of zero expansion, Nature, 1999, 400(6739), 18–19.
3. A. Sleight, Zero-expansion plan, Nature, 2003, 425(6959), 674–676.
4. C. E. Guillaume, Recherches sur les aciers au nickel. Dilatations aux temperatures elevees; resistance electrique, CR. Acad. Sci., 1897, 125(235), 18.
5. T. A. Mary, J. S. O. Evans, T. A. Vogt and W. Sleight, Negative Thermal Expansion from 0.3 to 1050 Kelvin in ZrW₂O₈, Science, 1996, 272, 90–92.
6. J. Lin, P. Tong, W. Tong, Y. Zou, C. Yang, F. Zhu, X. K. Zhang, L. F. Li, M. Wang and Y. Wu, et al., Large and constant coefficient of negative thermal expansion covering a wide temperature range in Zn_1−xMn_xNMn₃ (0 ≤ x ≤ 0.3), Scripta Mater., 2018, 152, 6–10.
7. Q. Gao, J. Wang, A. Sanson, Q. Sun, E. Liang, X. Xing and J. Chen, Discovering large isotropic negative thermal expansion in framework compound AgB(CN)₄ via the concept of average atomic volume, J. Am. Chem. Soc., 2020, 142(15), 6935–6939.
8. B. K. Greve, K. L. Martin, P. L. Lee, P. J. Chupas, K. W. Chapman and A. P. Wilkinson, Pronounced negative thermal expansion from a simple structure: cubic ScF₃, J. Am. Chem. Soc., 2010, 132(44), 15496–15498.
9. N. Lock, M. Christensen, C. J. Kepert and Bo. B. Iversen, Effect of gas pressure on negative thermal expansion in MOF-5, Chem. Commun., 2012, 49(8), 789–791.
10. A. L. Goodwin, M. Calleja, M. J. Conterio, M. T. Dove, J. S. Evans, D. A. Keen, L. Peters and M. G. Tucker, Colossal positive and negative thermal expansion in the framework material Ag₃[Co (CN)₆], Science, 2008, 319(5864), 794–797.
11. H. Lou, X. Wang, Q. Cao, D. Zhang, J. Zhang, T. Hu, H. Mao and J. Jiang, Negative expansions of interatomic distances in metallic melts, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 10068–10072.
12. A. K. Gangopadhyay, M. E. Blodgett, M. L. Johnson, J. McKnight, V. Wessels, A. J. Vogt, N. A. Mauro, J. C. Bendert, R. Soklaski and L. Yang, et al., Anomalous thermal contraction of the first coordination shell in metallic alloy liquids, J. Chem. Phys., 2014, 140, 044505.
13. S. V. Sukhomlinov and M. H. Müser, Determination of accurate, mean bond lengths from radial distribution functions, J. Chem. Phys., 2017, 146, 024506.
14. J. Ding, M. Xu, P. F. Guan, S. W. Deng, Y. Q. Cheng and E. Ma, Temperature effects on atomic pair distribution functions of melts, J. Chem. Phys., 2014, 140, 064501.
15. J. Ding and E. Ma, Computational modeling sheds light on structural evolution in metallic glasses and supercooled liquids, npj Comput. Mater., 2017, 9, 1–11.
16. L. Zhong, J. Wang, H. Sheng, Z. Zhang and S. X. Mao, Formation of monatomic metallic glasses through ultrafast liquid quenching, Nature, 2014, 512, 177–180.
17. M. Kbirou, A. Atila and A. Hasnaoui, On the structure and icosahedral interconnectivity in Tantalum monatomic glass produced under pressure, Phys. Scr., 2024, 99, 085946.
18. A. Atila, M. Kbirou1, S. Ouaskitb and A. Hasnaouic, On the presence of nanoscale heterogeneity in Al₇₀Ni₁₅Co₁₅ metallic glass under pressure, J. Non-Cryst. Solids, 2020, 550, 120381.
19. M. Kbirou, M. Mazroui and A. Hasnaoui, Atomic packing and fractal behavior of Al-Co metallic glasses, J. Alloys Compd., 2018, 735, 464–472.
20. D. Iabbaden, J. Amodeo, C. Fusco, F. Garrelie and J. Colombier, Dynamics of Cu–Zr metallic glass devitrification under ultrafast laser excitation revealed by atomistic modeling, Acta Mater., 2024, 263, 119487.
21. M. X. Li, Y. T. Sun, C. Wang, L. W. Hu, S. Sohn, J. Schroers, W. H. Wang and Y. H. Liu, Data-driven discovery of a universal indicator for metallic glass forming ability, Nat. Mater., 2022, 21, 165–172.
22. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys., 1995, 117, 1–19.
23. L. Gao, H. B. Yu, T. B. Schrøder and J. C. Dyre, Unified percolation scenario for the α and β processes in simple glass formers, Nat. Phys., 2025, 21, 471–479.
24. G. J. Martyna, D. J. Tobias and M. L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys., 1994, 101, 4177.
25. Y. Q. Jiang and P. Peng, Nearly golden-ratio order in Ta metallic glass, Chin. Phys. B, 2020, 29, 046105.
26. Y. Q. Jiang, J. Lv, W. X. He and P. Peng, Shed light on new growth pattern in rapid solidified Ta metal: MD simulation, DFT calculation and CALYPSO search, Vacuum, 2023, 216, 112423.
27. B. Delley, An all-electron numerical method for solving the local density functional for polyatomic molecules, J. Chem. Phys., 1990, 92(1), 508.
28. B. Delley, From molecules to solids with the DMol³ approach, J. Chem. Phys., 2000, 113(18), 7756.
29. M. Ernzerhof and G. E. Scuseria, Assessment of the Perdew–Burke–Ernzerhof exchange-correlation functional, J. Chem. Phys., 1999, 110(11), 5029.
30. P. K. John, B. Perdew and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865.
31. E. Räsänen, A. Castro and E. K. U. Gross, Electron localization function for two-dimensional systems, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 115108.
32. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson and M. C. Payne, First principles methods using CASTEP, Z. Kristallogr., 2005, 220, 567–570.
33. A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
34. Y. Q. Cheng and E. Ma, Atomic-level structure and structure–property relationship in metallic glasses, Prog. Mater. Sci., 2011, 56, 379–473.
35. Y. Q. Jiang and P. Peng, The influence of engineering strain rates on the atomic structure, crack propagation and cavitation formation in Al based metallic glass, J. Mater. Res. Technol., 2023, 26, 8997–9008.
36. L. Xu, Z. Wang, J. Chen, S. Chen, W. Yang, Y. Ren, X. Zuo, J. Zeng, Q. Wu and H. Sheng, Folded network and structural transition in molten tin, Nat. Commun., 2022, 13, 126.
37. J. L. Finney, Random packings and the structure of simple liquids. I. The geometry of random close packing, Proc. R. Soc. London, Ser. A, 1970, 319, 479.
38. W. Brostow, M. Chybicki, R. Laskowski and J. Rybicki, Voronoi polyhedra and Delaunay simplexes in the structural analysis of molecular-dynamics-simulated materials, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 13448.
39. D. Q. Yu, M. Chen and X. J. Han, Structure analysis methods for crystalline solids and supercooled liquids, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 051202.
40. Y. D. Wei, P. Peng, Z. Z. Yan, L. T. Kong, Z. A. Tian, K. J. Dong and R. S. Liu, A comparative study on local atomic configurations characterized by cluster-type-index method and Voronoi polyhedron method, Comput. Mater. Sci., 2016, 123, 214–223.
41. L. Pauling, Atomic radii and interatomic distances in metals, J. Am. Chem. Soc., 1947, 69, 542.
42. A. D. Becke and K. E. Edgecombe, A simple measure of electron localization in atomic and molecular systems, J. Chem. Phys., 1990, 92, 5397–5403.
43. F. Peng, Y. Yao, H. Liu and Y. Ma, Crystalline LiN₅ Predicted from First-Principles as a Possible High-Energy Material, J. Phys. Chem. Lett., 2015, 6, 2363–2366.
44. R. S. Mulliken, Criteria for the construction of good self-consistent-field molecular orbital wave functions, and the significance of LCAO-MO population analysis, J. Chem. Phys., 1962, 36, 3428–3439.
45. Y. Q. Jiang, P. Peng, D. D. Wen, S. C. Han and Z. Y. Hou, A DFT study on the heredity induced coalescence of icosahedral basic clusters in the rapid solidification, Comput. Mater. Sci., 2015, 99, 156–163.

---

Footnote

† Electronic supplementary information (ESI) available: Fig. S1 shows the systemic potential energy, volume and pressure of Ta_11 664 as a function of the temperature at five different cooling rates, respectively. Fig. S2 shows the radial distribution function (RDF) curves of the simulated Ta_11 664 system at 300 K under five different cooling rates, respectively. See DOI: https://doi.org/10.1039/d5cp00247h

---