The local coordination behavior of cation–anion is depicted using the partial radial distribution function (RDF) and cation–anion coordination numbers. Figure 1 compares the cation–anion partial RDFs as obtained from AIMD and IPMD simulations for each composition. In all cases, the RDF peak for the F–Zr pair is more intense and has a smaller width, which suggests a strong association between F and Zr in the melt. Much broader peaks for F–Na and F–Li suggest comparatively weak interactions among the respective cation–anion pairs. A very low minimum in the F-Zr RDF plot indicates strong solvation with limited exchange of fluorine from the Zr first solvation shell. However, a much higher minimum in the F–Na and F–Li plot suggests weak solvation allowing a rapid exchange of “free” and coordinated fluoride ions. Also, a smaller value of the first peak radius, r _avg,F-Zr = 2.04 Å, further supports a stronger association between F and Zr in the molten salt mixture. Similar observations have been made by Salanne et al. for a different LiF–NaF–ZrF₄ salt composition studied using RIM potential parameters [14].

Overall, the peak positions in the cation–anion RDFs are similar for both AIMD and IPMD calculations. However, the F–Zr peak heights for IPMD calculations are higher than those of AIMD, resulting in higher F–Zr coordination numbers as reported in Table 1, which is discussed in the upcoming section.

Table 1 reports and compares the average F–Zr coordination numbers (CN) and the % of individual complexes in the AIMD- and RIM-IPMD-simulated structure. CN calculation is based on r _cut,F-Zr = 2.84 Å, which is the first minimum in the F-Zr RDF curve. The coordination number was averaged over all Zr atoms over equilibrated trajectories. Supplementary Figure S1 can be referred to for the %CN complexes plotted over the equilibrated trajectory from AIMD for composition C. In AIMD calculations, the Zr coordinations of 5, 6, 7, and 8 were observed as [ZrF₅]^- [ZrF₆]^2- [ZrF₇]^3-, and [ZrF₈]^4-, respectively, during the 60-ps simulation trajectory. The 5-coordinated species are short-lived (<200 fs) and may have occurred due to thermal vibrations considering the higher simulation temperatures. It is in agreement with the Zr coordination values of 6, 7, and 8 obtained from the Raman spectrum measurements of LiF–NaF–ZrF₄ melt with ZrF₄ mol% varying from 14% to 40% [27]. On the other hand, Zr coordination numbers of 7, 8, and 9 were observed in RIM-IPMD calculations with [ZrF₇]^3-, and [ZrF₈]^4- occurring as dominating species. Overall, the AIMD melt consists of 6 (octahedral), 7, and 8 coordination zirconium fluoride complexes (Figure 2A), while RIM simulated trajectories comprised 7, 8, and 9 coordination zirconium fluoride complexes (Figure 2B).

The average Zr–F coordination number is decreasing as the ZrF₄ mol% increases, which is observed from both AIMD and RIM IPMD calculations. This is caused by the decreasing F/Zr ratios with the increase in Zr mol%. As ZrF₄ is more fluoroacidic, it binds more readily to the “free” fluorines, resulting in higher values of average Zr coordination number [28]. With the increase in number of Zr in the melt, fewer F per Zr are available to form a higher coordinated complex, while when the Zr concentration in the mix is less, more F atoms are available per Zr atom which leads to formation of a higher coordinated complex as observed in composition A in both AIMD and RIM IPMD calculation. This is also evidenced by an increase in the percentage of higher coordinated complexes as ZrF₄ mol% goes down (Table 1). However, this trend is similar in AIMD and RIM calculations, the average Zr–F coordination number predicted by RIM is consistently higher than that calculated in AIMD for all three compositions. It can be noted that deviation in RIM and AIMD-calculated Zr CN grows from 15% to 19% going from compositions A to C (increasing ZrF₄ mol%). The occurrence of higher Zr–F coordination numbers in RIM simulations was also previously pointed out based on the F–Zr cation–anion RDF discussion.

Some inconsistencies in the coordination number of Li and Na are also observed between AIMD and RIM-generated structures, which are reported in Table 2. The Li–F and Na–F cutoff distances were calculated based on the first minimum in their respective RDF curves, which are also reported in the same table. For both AIMD and RIM, the cutoff distances agree well and indicate a shift to the right going from compositions C to A. It should be noted that the shift could also be caused due to relatively higher temperatures considered for compositions A and B. The reported CN for both monovalent cations is calculated by averaging over equilibrated trajectories. For lower ZrF₄ mol% as in composition A, the differences in the CN of both cations are less than 2%. For 37 ZrF₄ mol% (composition C), the differences between RIM and AIMD values of Li and Na CN grow to 7%. This trend is similar to that observed for the Zr–F CN obtained from the AIMD and RIM study. Nevertheless, the deviation in CN of Zr (15%–19%) is higher compared to that observed for monovalent cations (up to 7%).

The differences in coordination behavior can be explained by examining the repulsion and columbic terms in the RIM potential. In general, the short-range structural properties are dominated by the competition between the overlapping repulsion and the Columbic interactions. Specifically, the electrostatic term has a role in inducing the ordering effects around an ion, which decides the cation-ion coordination number. As the RIM-calculated CNs disagree with the reference AIMD values, a refitting of the mentioned short-range terms could lead to better prediction of CNs [29]. However, such refitting would likely result in loss of generality of the RIM toward other compositions and mixtures of LiF–NaF–ZrF₄. In our study, the higher deviation (between AIMD and RIM) in CN of Zr (15%–19%) compared to monovalent cations (up to 7%) indicates a comparatively poorly fitted short-range interaction terms for F–Zr interactions.

In addition to cation–anion RDFs, the Zr–Zr RDFs are also examined in order observe the occurrence of fluorozirconate network formation (Figure 3). The pre-peak between 3 and 5 Å in the Zr–Zr RDF plot occurred due to the bridging fluorine connecting the two Zr atoms from two adjacent [ZrF_x ^4−x] polyhedrons. It should be noted that the pre-peak in Zr–Zr RDF grows stronger with the increase in ZrF₄ mol% for both AIMD and RIM calculations. Overall, the Zr–Zr pre-peak observed in RIM calculations is quite sensitive to the ZrF₄ mol% when compared to AIMD calculations. In the case of higher ZrF₄ mol% as in compositions B and C, RIM-IPMD Zr–Zr pre-peak is quite prolific compared to that observed in the AIMD calculation, indicating an overprediction of network formation in the RIM simulations. To our knowledge, the intermediate-range salt structure and Zr–Zr RDF have not been reported in the literature when ZrF₄ mol% is higher than 10% in this ternary salt system.

The Zr–Zr pre-peak completely disappears for composition A in which the ZrF₄ concentration is 11%, which suggests no chain formation in AIMD simulation. A similar behavior is observed in the RIM-IPMD Zr–Zr RDF plot. A much smaller pre-peak observed in our RIM simulations for lower ZrF₄ is also in agreement with a study from Rollet et al. [30]. A close agreement between RIM and AIMD structure for composition A indicates that the RIM parameters are acceptably reproducing the salt structure in case of low Zr mol% in the salt melt, i.e., when polarizability effects are relatively mild. Both RIM-IPMD and AIMD simulations suggest that at lower concentrations of ZrF₄ in this ternary salt mix, the mixture tends to behave as a well-dissociated ionic melt consisting of Li⁺, Na⁺ [ZrF_x ^4−x], and F⁻ ions.

Based on RDF and coordination analysis, the inclusion of ZrF₄ in the salt melt forms the [ZrF_x ^4−x] complexes in the melt. Thereafter, the increase in ZrF₄ mol% leads to an increase in the connection formed by [ZrF_x ^4−x] leading to the formation of larger [ZrF_x ^4−x]_n chains, as is realized in Figure 3. The occurrence of fluorozirconate chains in this study is in agreement with the chain formation reported from NMR spectra of LiF–ZrF₄ melts [31].

The angular distribution functions for F–Zr–F and Zr–F–Zr provide information on the individual [ZrF_x ^4−x] complexes and their connectivity in the [ZrF_x ^4−x]_n chain, respectively. Figures 4A,B shows and compares the F–Zr–F angle distribution between three compositions as obtained from AIMD and RIM-IPMD calculations. The first peak in the F–Zr–F angle distributions provides the angle formed by adjacent fluorines with the central Zr in a fluorozirconate complex. With the decrease in Zr concentration (increase in F/Zr ratio), the first peak in the F–Zr–F angle distribution plot shifts to the left from the peak value of 81.9°–76.5° between composition C and composition A as observed for AIMD calculations (Figure 4A). This shift points to a decrease in the mean value of the F–Zr–F angle due to increased average Zr coordination number. This observation is consistent with the Zr coordination values reported in Table 1. A similar behavior is observed in the plot obtained from RIM-IPMD calculations (Figure 4B). However, due to the higher average Zr CN reported by RIM, the first peak is located at angles lowered by up to 8° compared to the corresponding values obtained from AIMD for each composition. For RIM calculations, the first peak shift between composition C and composition A is negligible, which is also corroborated by a small change in average Zr coordination number between compositions C and A as obtained from RIM calculations (Table 1).

Figure 5 displays the 2D density plots for θ _Zr-F-Zr plotted against Zr–Zr distances (d _Z-Zr ) in an AIMD simulation for compositions B and C. Here, θ _Zr-F-Zr angles are formed between two fluorozirconate complexes connected by a bridging fluorine atom, while d_Z-Zr is the distance between the Zr atoms from the connected fluorozirconates. The histograms along the abscissa and ordinate are the addition of all the bin values along the corresponding directions. The two distinct peaks in θ _Zr-F-Zr angle distribution corresponding to the high-density region in the corresponding density plot indicate the presence of edge-sharing and corner-sharing complexes as observed in AIMD simulation. The number of face-sharing complexes is almost negligible for both the compositions as indicated by a sparse region in the density plot and supported by no distinctive third peak in the angle distribution. The corner-sharing complex average connectivity angle is ∼145° while the angle is nearly 110° for edge-sharing complexes. Both compositions B and C show a network connected via corner sharing and edge sharing. No such analysis was done for composition A as no chain formation (hence, no bridging fluorine) was observed in this case resulting in a negligible number of short-lived Zr-bridging F–Zr triplets. Compared to composition C, a lower number density is observed for composition B due to a low number of Zr-bridging F–Zr triplets (Figure 4).

The θ _Zr-F-Zr and d _Zr-Zr plot for RIM IPMD calculations (Figure 6) looks significantly different from that obtained from AIMD calculations. The presence of only one major peak in θ _Zr-F-Zr and d _Zr-Zr histograms corresponds to one dense region in the density plot. This implies the absence of edge-sharing connectivity between the two adjacent fluorozironates in the network for compositions B and C. Even though Figures 3B,C suggested an overprediction of the fluorozirconate network in RIM-IPMD simulations, the whole network is largely connected by corner-sharing complexes. It could be caused due to the exclusion of higher-magnitude dipole moment values during the dipole-fitting process when optimizing the parameters for the polarization term by Salanne et al. [13]. In multivalent compounds such as ZrF₄, the polarization of F⁻ plays a key role in local structure by screening the repulsive electrostatic interactions between two Zr⁴⁺ cations in the anion coordination shell. The exclusion of higher-magnitude dipole moment values during potential fitting may result in less F⁻ polarization effect leading to less shielding between Zr⁴⁺ and Zr⁴⁺ connected by bridging fluorine, which in turn results in larger values of θ _Zr-F-Zr and d _Zr-Zr (corner-sharing zone). Consequently, a low-density region is observed at smaller values of d _Zr-Zr and θ _Zr-F-Zr (edge-sharing zone).

Qualitatively, intermediate-range fluorozirconate chain formation can be visualized from equilibrated AIMD and RIM simulation trajectories. Figure 7 and Figure 8 show the snapshots from AIMD and RIM simulations, respectively. For both compositions B and C, AIMD shows the formation of fluorozirconate chains ([ZrF_x ^4−x]_n) constructed by edge- and corner-sharing polyhedrons as observed in Figures 7B,C. In comparison, the RIM-generated network is mainly connected by the corner-sharing polyhedrons in both compositions B and C (Figures 8B,C). Such differences in network connectivity were also observed from Figure 5 and Figure 6 based on the θ _Zr-F-Zr vs. d _Zr-Zr density plots. For compositions B and C, Figures 8B,C confirm an overprediction of network by RIM when compared to AIMD snapshots (Figures 7B,C). For lower 11 mol% ZrF₄ (composition A), both AIMD and RIM snapshots (Figure 7A and Figure 8A) show no such long-lived chain formation.

In order to quantify the fluorozirconate chain lengths, the number of Zr atoms (representing count of fluorozirconate complexes) in the two largest chains (Chain 1: longest chain and Chain 2: second longest chain) was counted. Supplementary Figure S2 plots and compares the values of number of Zr involved in the longest chain (a.k.a., Chain 1) formation for each composition during equilibrated simulation trajectories from the respective AIMD and RIM-IPMD calculations. The values averaged over AIMD and RIM trajectories for both chains (Chain 1: longest chain and Chain 2: second longest chain) are reported in Table 3. For composition C, AIMD predicts ∼67% of Zr atoms (i.e., 6–7 Zr atoms out of a total of 10) involved in the longest chain formation, while RIM-IPMD calculations indicate more than 98% Zr (i.e., 79 Zr atoms out of a total of 80) involved in the longest chain formation. Similarly for 29 ZrF₄ mol% as in composition B, AIMD predicts ∼26% of Zr atoms (i.e., two to three Zr atoms out of a total of 9) involved in the chain formation, while RIM-IPMD calculations indicate more than 92% Zr (i.e., 66–67 Zr atoms out of a total of 72) involved in the chain formation. The decrease in % of Zr in the network formation in composition B compared to composition C can also be noticed by the comparatively lower number density of Zr-bridging F–Zr triplets in Figure 5A and Figure 6A compared to Figure 5B and Figure 6B, respectively.

Both AIMD and RIM suggest no long-lived chain formation for lower ZrF₄ mol% as can be noticed in the snapshots for composition A in Figure 7A and Figure 8A. For instance, Table 3 reports that AIMD showed the longest chain (Chain 1) to be consisting of nearly 1 Zr, indicating the presence of isolated fluorozirconate complexes in composition A. On the other hand, a value of 1.36 for average Zr atoms involved in the longest chain indicates a short-lived chain formation including between one and two fluorozirconate complexes in the RIM simulated trajectories.

Overall, the network structure predicted by RIM for salts with lower ZrF₄ mol% (composition A) is in close agreement with that obtained from AIMD, whereas there is an overprediction of network formation when using RIM parameters for higher ZrF₄ mol% (compositions B and C). In addition, for both compositions B and C, RIM shows similar motifs containing a continuous fluorozirconate network spread across the periodic boundaries of a larger simulation cell (Figures 8B,C), which may affect the transport and thermophysical properties of the salt melt as discussed in Section 3.2 and Section 3.3. It should be noted that during the RIM force-field fitting, the dispersion parameters were not fitted directly from the first-principle calculations and instead were adjusted later to reproduce the experimental density of the pure compounds. It is known that dispersion parameters along with the polarization term also play a crucial role in the intermediate salt structure of multicomponent salts involving multiply charged species [11]. Recently, the discrepancies in the cation coordination numbers and structures due to inaccurate rigid ion model (RIM) parameters for NaCl and FLiNaK molten salt were pointed out by Lee et al. [29], which were fixed upon reparametrizing the RIM force-field parameters. While it suggests that a similar re-parameterization of RIM parameters for the LiF–NaF–ZrF₄ salt system while implicitly including the dispersion effects may lead to obtaining an accurate short to the intermediate-range salt structure, the transferability of the refitted force-field parameters would be compromised. In quest to obtain an AIMD-accurate salt structure in comparatively efficient molecular dynamics simulations, recent success of AIMD-trained neural network potentials for salts containing divalent Be²⁺ cations [32, 33] also provides motivation for the use of machine-learned force-field for ZrF₄ molten salts.

The equilibrated densities for compositions A, B, and C at 750°C, 727°C, and 700°C, respectively, are reported in Figure 9. For each case, RIM-IPMD-calculated densities are higher than the AIMD-calculated values. Both AIMD and the additive molar volume calculation of density are higher than the available experimental data for composition C at 700°C. For this composition, AIMD predicts the density nearly 11% higher than the experimental value while RIM overpredicts the density by ∼26%. The density prediction from AIMD is previously reported to be greatly influenced by the choice of dispersion interactions included in the DFT calculation [35]. A comparative study involving different dispersion methods such as vdW-DF, DFT-D3, and DFT-D2 can be nontrivial. As the current AIMD density is in reasonable range corresponding to experimental values, no such comparison is included. Regardless, the RIM parameters are further overpredicting the density by ∼26% compared to the experimental value. No uncertainty values were reported for the experimental measurement, while the densities estimated from the additive molar volume are considered to be within 5% of the experimental values [36]. The uncertainty values for both RIM and AIMD are based on the 95% confidence interval of the mean estimate. No simulation as well as experimental data for densities are available for comparison for compositions A and B.

In RIM-IPMD simulations, the Green–Kubo method is used to evaluate the viscosity by calculating the time integral of the shear stress autocorrelation function (ACF) as shown in Eq. 3: η = V k B T ∫ 0 ∞ σ α β ( t ) . σ α β ( 0 ) d t (3) where V is the volume of the simulation cell, k _B is the Boltzmann constant, T is the temperature, and σ _αβ is one of the components of the stress tensor. For better statistics, the ACF is averaged over five stress components being σ _xy , σ _xz , σ _yz , σ _xx-yy , and σ _2zz-xx-yy . The stress tensor data were obtained from RIM-IPMD run from nearly 10-, 8-, and 24-ns simulations for compositions A, B, and C, respectively, at their equilibrium densities. In each case, ACF is calculated by dividing the simulation trajectory in two halves. It can be noted that even after significantly longer simulation times, the stress ACF is not converging to zero, hence delaying the plateau behavior in the viscosity values. It is also observed by other researchers while simulating a highly polymerized system such as LiF-BeF₂ [37]. Figure 10 reports the ACF and viscosity values evaluated from Eq. 3 for each composition using RIM parameters.

It can be noted that the viscosity values for composition C drastically surpass its experimental value of 6.9 cP reported at 700°C [38] by nearly two orders of magnitude. It can be explained from the perspective of long-lived extensive network formation (Figure 8) and significantly lower Zr self-diffusivities (Table 4) observed from RIM-IPMD calculations. Clearly, an unrealistic overprediction of network formation significantly affects the viscosity of the melt. Similarly, as composition B also showed extensive network formation by involving nearly 93% of Zr in the longest chain formation (Table 3), the predicted viscosity values are still higher by an order of magnitude. Due to a better agreement between AIMD and RIM-IPMD regarding network formation in composition A (11 mol % of ZrF₄), the viscosity value is nearly 10 cP. As such, no experimental data for viscosities are available for comparison for compositions A and B. As previously observed for Flibe [37] and NaF-ZrF₄ [38], viscosity values tend to decrease with the decrease in the mol% of network-forming multivalent cations. Following such trend, the viscosity values for composition A would be expected to be lower than the experimental value of 6.9 cP for composition C. A comparatively higher value of viscosity calculated from the RIM simulation of composition A can be due to the higher coordination number of Zr (Table 1) and short-lived chain formation involving ∼27% of Zr in the melt (Table 3).

The heat capacity values for all three compositions are evaluated from RIM-IPMD calculations from nearly 2-ns NPT calculations performed at five different temperatures. For composition A, NPT simulations were performed at 750°C–950°C with 50° intervals. For composition B, specific enthalpies were calculated from NPT simulations performed at 727°C, 777°C, 800°C, 850°C, and 900°C, while 700°C, 750°C, 800°C, 850°C, and 900°C temperatures were considered for composition C. For each calculation, the pressure was targeted to 1 atm. The specific enthalpies were calculated using ACF, and the error bars were estimated based on the 95% confidence interval from the mean estimated values. For compositions A and B, the specific enthalpies follow a linear trend when calculated between the given temperature range as shown in Figures 11A,B, indicating a constant value of C _p within the temperature range. Weighted least square regression was used to obtain a linear trend between specific enthalpies and temperature for each composition. A similar linear trend in specific enthalpies and temperature for short temperature ranges (200°C–300°C) at higher temperatures (>T_melt) was noticed in previous literature [13, 39–43], while a nonlinear relationship was obtained when fitting to a wider temperature range [37, 38, 44, 45]. For composition C, a second-order polynomial fit better described the specific enthalpies and temperature trend in temperature range of 700°C–900°C (Figure 11C). The slope of the fitted specific enthalpy vs. temperature correlation gives specific heat capacity as described in Eq. 4: C p = ( ∂ H ∂ T ) P       (4) C p ,   c o m p o s i t i o n   C = ( 8.285 ± 6   ) × 10 − 4   T ( ° C ) − 0.4823 ± 0.477           ( 700 − 900 ° C ) (5)

The specific heat capacities values of 0.256 ± 0.002 and 0.24 ± 0.006 cal/gm-°C were obtained for compositions A and B, respectively. It should be noted that the constant values of specific heat capacity only represent the approximate C _p vs. T trend in the considered temperature ranges. For composition C, Eq. 5 describes the linear correlation between C _p and T in the temperature range of 700°C–900°C. The reported uncertainty values are evaluated while fitting the equation to the specific enthalpy values using weighted least square regression. Previous experimental measurements using differential scanning calorimetry for 42 LiF–29 NaF–29 ZrF₄ (composition B) reported a C _p value nearly 0.334 cal/gm-°C with the 5% uncertainty [46]. In this case, the C _p value evaluated from RIM-IPMD is ∼28% smaller compared to the reported experimental values. To our knowledge, no previous experimental data are reported for the other two compositions studied in this paper. Overall, the trend of decrease in the C _p values with the increased ZrF₄ mol% is in agreement with the observations from experimental measurements for the LiF-KF-ZrF₄ system [47].