Short- to Intermediate-Range Structure, Transport, and Thermophysical Properties of LiF–NaF–ZrF4 Molten Salts

LiF–NaF–ZrF4 multicomponent molten salts are identified as promising candidates for coolant salts in molten salt reactors and advanced high-temperature reactors. This study focused on low-melting point salt compositions of interest: 38LiF–51NaF–11ZrF4, 42LiF–29NaF–29ZrF4, and 26LiF–37NaF–37ZrF4. Ab-initio molecular dynamics (AIMD) calculations were performed and compared with available experimental data to assess the ability of rigid ion models (RIM) to reproduce short- to intermediate-range structure, transport, and thermophysical properties of the LiF–NaF–ZrF4 salt mixtures. It is found that as ZrF4 mol% increases, the average cation–anion coordination number (CN) of monovalent cations (Li+, Na+) obtained from RIM calculations decreases, while multivalent Zr4+ CN varied from 15% to 19% in comparison to corresponding AIMD values. In addition, RIM is found to predict the existence of 7, 8, and 9 coordinated fluorozirconate complexes, while AIMD and the available experimental data showed an occurrence of 6, 7, and 8 coordinated complexes in the melt. The intermediate-range structure analysis revealed that while the RIM parameters are able to reproduce a local structure for lower ZrF4 mol% salts such as in 38LiF–51NaF–11ZrF4, an extensive fluorozirconate network formation is observed in RIM simulations for higher ZrF4 mol% compositions. The network generated by RIM parameters is found to be mainly connected by “corner-sharing” fluorozirconate complexes as opposed to both “edge-sharing” and “corner-sharing” connectively portrayed by AIMD. It is found that a close agreement between AIMD and the RIM salt structure for the 11-mol% ZrF4 salt resulted in good agreement in the calculated Zr diffusivities and the viscosity values. However, due to the inaccurate short- to intermediate-range structure prediction by RIM for higher ZrF4 mol% compositions, thermophysical properties such as densities and heat capacity differ by up to 26% and 27%, respectively, upon comparison with AIMD and experimental values. Also, the network-dominated properties such as diffusion coefficients and viscosities differed by up to two and three orders of magnitude, respectively. This study signifies the importance of accurate salt structure generation for an accurate prediction of transport and thermophysical properties of multicomponent molten salts.


INTRODUCTION
Molten salts are a promising class of ionic liquids for applications in advanced clean energy systems such as next-generation nuclear reactors and solar-thermal storage plants. Fluoride salts have been previously identified as good candidates for primary coolant applications in advanced high-temperature reactors (AHTR) [1] and molten salt reactors (MSR) [2]. Among them, ternary systems containing BeF 2 and ZrF 4 were recommended as coolant salts. While there has been significant interest in using BeF 2 -based salts due to their low neutron absorption, there remain substantial issues in using beryllium salts due to their toxicity and required processing facilities. As such, Zr salts present a compelling alternative due to their acceptable neutron economy, vapor pressures, thermal hydraulics, and lower costs [3]. In order to achieve low vapor pressure at higher temperatures (<1 mm Hg at 700°C-900°C), the ZrF4 mole fraction in the salt mixture should be maintained withiñ 20%-45% [3].
Here, the eutectic compositions 26LiF-37NaF-37ZrF 4 (mol%) and 42LiF-29NaF-29ZrF 4 (mol %) with a freezing point around 436°C and 460°C, respectively, have been identified as promising candidates. The 38LiF-51NaF-11ZrF 4 (mol%) salt composition is also considered in this study due to low ZrF 4 mol%. Due to lack of data in thermophysical property databases [4], this salt composition has been recommended for further study [2]. The precise experimental interrogation of salt structure and properties is challenged by high-temperature conditions, cost, material handling, and difficulties in interpreting experimental data. Particularly, it is difficult to access the structure of multivalent cations using techniques such as Raman spectra and EXAFS alone due to their existence in multiple coordination states and intermediate-range ordering [5][6][7]. Therefore, ab-initio and classical molecular dynamics simulations can be used to interpret experimental data and predict temperaturedependent structures, transport, and thermophysical properties of LiF-NaF-ZrF 4 molten salts.
AIMD simulations have been shown to accurately predict the salt structure and transport properties of binary and ternary fluoride salts [8,9]. This provides further motivation for using AIMD to study the structure, transport, and thermophysical properties of LiF-NaF-ZrF 4 molten salt. Despite the accurate salt structures and property prediction by AIMD, the computational cost of AIMD simulations limits their use in terms of system size (a few hundred atoms) and timescale (picosecond). Therefore, the classical molecular dynamics calculations using a well-fitted interatomic potential can be quite useful. In classical molecular dynamics modeling of ionic liquids involving multivalent charged species, accurately capturing the charge-dipole and dipole-dipole polarizability effects is essential for generating an accurate salt structure [10]. As the properties such as diffusion coefficients, conductivity, and viscosity of the salt melt are strongly influenced by the formation of coordination complexes, their lifetime, and the degree of their connectivity (network formation) [11], it is crucial to regenerate the correct short-to intermediaterange structure using the interatomic potential in classical molecular dynamics simulations. However, correctly accounting for all pair interactions in multivalent cations can be a tedious and challenging process as it depends on the quality of both force fitting and dipole fitting. Additionally, the quality of experimental and first-principle data used for parameter optimization further affects the accuracy of the developed rigid ion model force field.
In the past, the rigid ion model potential has shown great promise in the simulations of multicomponent molten salts involving multivalent cations [11,12]. Efforts have been made to develop RIM parameters for the configurations LiF, NaF, KF, ZrF 4 , LiF-ZrF 4 , and NaF-ZrF 4 using a force-fitting procedure [13]. While the fitted RIM parameters were used by Salanne et al. [14] to study electrical conductivities of molten LiF-NaF-ZrF 4 mixtures, a detailed study on the salt structure and its effect on the properties such as heat capacity and viscosity is yet to be reported.
In this study, classical molecular dynamics using the rigid ion model [13] was used to first explore the salt structure of 38LiF-51NaF-11ZrF 4 , 42LiF-29NaF-29ZrF 4 , and 26LiF-37NaF-37ZrF 4 at 750°C, 727°C, and 700°C, respectively. These predictions were evaluated against Born-Oppenheimer ab-initio molecular dynamics simulations. The first part of this study focuses on determining and comparing the short-to intermediate-range structure/network formation in the LiF-NaF-ZrF 4 salt melt using both AIMD and RIM interatomic potential molecular dynamics (IPMD) simulations. In order to study the effect of structures on properties, diffusivity and thermophysical properties are evaluated. In doing so, the second part incorporates the evaluation of ionic diffusivity, diffusion coefficients, densities, viscosity, and heat capacity. The RIMcalculated values are compared with the calculated AIMD values and available experimental data. The effect of salt structure generated by RIM parameters on the calculated properties is discussed in each section.

Ab-Initio Molecular Dynamics Simulations
Born-Oppenheimer AIMD simulations were performed using the Vienna Ab-Initio Simulation Package (VASP) with the projector-augmented wave (PAW) method, a plane-wave basis set, and the Perdew-Burke-Ernzerhof (PBE) generalizedgradient-approximation (GGA) exchange correlation functional [15,16]. PAW-PBE potentials provided by VASP were used for Li_sv (1s 2 2s 1 ), Na_sv (2s 2 2p 6 3s 1 ), Zr_sv (4s 2 4p 6 4d 2 5s 2 ), and F_s (2s 2 2p 5 ). A plane-wave basis set with an energy cutoff of 650 eV was used. The convergence criterion of 1E-5 eV was set for electronic self-consistent steps. A gammacentered 1 × 1 × 1 k-point mesh was used for reciprocal space sampling. The parameters chosen yielded convergence within 2 meV/atom and is in agreement with previous studies [17]. Charges were calculated within the VASP code. The density functional theory (DFT)-D3 formulation proposed by Grimmes [18] was used to account for the effect of dispersion interactions. Three different compositions of LiF-NaF-ZrF 4 mixtures are considered in this study. The salt systems were generated by randomizing the atom positions in a simulation box with periodic boundary conditions [19]. Simulation cells contained 17 Li, 23 Na, 5 Zr, and 60 F atoms (Composition A: 38LiF-51NaF-11ZrF 4 ), 13 Li,9 Na,9 Zr, and 58 F atoms (Composition B: 42LiF-29NaF-29ZrF 4 ), and 7 Li, 10 Na, 10 Zr, and 57 F atoms (Composition C: 26LiF-37NaF-37ZrF 4 ). All calculations were performed allowing for spin polarization. The canonical ensemble (NVT) using a Nosé-Hoover thermostat [20] was employed. The integration time step of 2 fs was used for all the calculations. All simulations were run for 40-80 ps, and the system is allowed to equilibrate for at least 20 ps before the structure and property evaluation. Reported equilibrium densities were calculated based on the equation of state fit on 5-6 volumes, as described in [21]. The calculated equilibrated densities obtained from the AIMD are reported in Section 3.3.1.

Classical Molecular Dynamics Modeling Using Rigid Ion Model Potential
The rigid ion model potential is used to study the three compositions of interest and was performed using CP2K [22]. The simulation cell size for RIM-IPMD calculations for all three compositions was eight times that of the AIMD simulations. The interatomic interactions are defined using RIM potential developed in Ref. [13]. Formal charges were used for F (-1), Li (+1), Na (+1), and Zr (+4). For a direct comparison with the first-principle calculations, the simulations were conducted at densities estimated by both additive molar volume density and equilibrium AIMD calculations. These simulations were used for sampling from over 4-ns trajectories each. Additionally, for each composition, the equilibrium system densities were also obtained by equilibrating the salt in the NPT ensemble fixed at atmospheric pressure using a Nosé-Hoover thermostat and barostat [20].The equilibrated densities obtained from the RIM parameters are reported in Section 3.3.1. Ovito [23] and Vesta [24] were employed for visualization.

Salt Structure Analysis
The local coordination behavior in the salt was analyzed and compared using radial distribution function and cation-anion coordination numbers. The Zr-Zr RDFs were examined to observe the presence of fluorozirconate network formation. Further, the angle distribution function was introduced to discuss the fluorozirconate complex connectivity. The intermediate-range structure was examined by snapshots from both RIM and AIMD trajectories, and the fluorozirconate chains were quantified to evaluate the percentage of Zr involved in them. As the ionic charges were not estimated from AIMD, the specified charges while discussing AIMD results only refer to the formal charge on ions for convenience such as that used in [8,25,26].

Local Coordination Behavior
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 4 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 AIMDand 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 5 ] -[ZrF 6 ] 2-[ZrF 7 ] 3-, and [ZrF 8 ] 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 4 melt with ZrF 4 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 7 ] 3-, and [ZrF 8 ] 4occurring 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 4 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 4 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 4 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 AIMDcalculated Zr CN grows from 15% to 19% going from compositions A to C (increasing ZrF 4 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 4 mol% as in composition A, the differences in the CN of both cations are less than 2%. For 37 ZrF 4 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 4 . 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 shortrange 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 4 mol% for both AIMD and RIM calculations. Overall, the Zr-Zr pre-peak observed in RIM calculations is quite sensitive to the ZrF 4 mol % when compared to AIMD calculations. In the case of higher ZrF 4 mol% as in compositions B and C, RIM-IPMD Zr-Zr prepeak 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 4 mol% is higher than 10% in this ternary salt system.
The Zr-Zr pre-peak completely disappears for composition A in which the ZrF 4 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 4 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 4 in this ternary salt mix, the mixture tends to behave as a welldissociated ionic melt consisting of Li + , Na + [ZrF ] 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 4 melts [31].

Angular Distribution Function: Fluorozirconate Complexes Connectivity
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°b etween 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°c ompared 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 4 , the polarization of F − plays a key role in local structure by screening the repulsive electrostatic interactions between two Zr 4+ 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 4+ and Zr 4+ connected by bridging fluorine, which in turn results in larger values of θ Zr-F-Zr and d Zr-Zr (cornersharing zone). Consequently, a low-density region is observed at smaller values of d Zr-Zr and θ Zr-F-Zr (edge-sharing zone).

Intermediate-Range Fluorozirconate Chain Formation
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 4 (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  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 4 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 Zrbridging 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 4 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 4 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 4 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 4 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 AIMDtrained neural network potentials for salts containing divalent  Be 2+ cations [32,33] also provides motivation for the use of machine-learned force-field for ZrF 4 molten salts.

Self-Diffusion Coefficients
In order to study the effect of salt structure on the transport properties, the diffusion coefficients of Zr 4+ are evaluated from AIMD and RIM simulations. In order to obtain diffusion coefficients, mean squared displacement (MSD) calculated using the block-averaging method (Eq. 1) [34] was used in the Einstein equation (Eq. 2). In the block-averaging scheme, the equilibrated simulation trajectory was divided into nearly 10 equal blocks (n t ). The MSD vs. time curve was fitted for all but initial steps to avoid the initial quadratic region in the MSD plot. The remaining data are fitted to obtain diffusion coefficients for each block. The final diffusion coefficient is calculated as the average of values from each block. The error on D Zr4+ for each block is calculated based on the 95% confidence interval and is finally averaged over the block count. It should be noted that the relative uncertainties are smaller for RIM-IPMD calculations due to better sampling in the larger simulation cell. Table 4 compares the diffusion coefficients for Zr 4+ as calculated from AIMD and RIM-IPMD simulations. For direct comparison, the diffusion coefficients reported in Table 4 were evaluated at the same densities (ρ AIMD ) for both AIMD and IPMD calculations, while the values calculated at RIM equilibrated densities (ρ RIM ) are reported in Supplementary  Table S1. For composition A, the RIM-calculated diffusion coefficient is nearly 56% lower than the AIMD value, which is caused due to higher Zr coordination complexes ( Table 1) and short-lived chain formation involving~27% of Zr on average ( Table 3). Zr diffusion coefficient values for compositions B and C are lower by one and two orders of magnitude, respectively, compared to the corresponding AIMD values. A previous study involving a smaller system size comprising 574 atoms for composition B [14] reported the Zr 4+ diffusion coefficient as nearly 5 × 10 -6 cm 2 /s, which differs slightly from the value of 0.43 ± 0.14×10 -6 cm 2 /s calculated in this study. The observed discrepancies may be caused due to the smaller system size and comparatively shorter simulation length (~2 ns) in [14], which may have resulted in higher sampling uncertainty (not reported in Ref. [14]). This can have a significant effect where diffusivities are relatively low due to the formation of intermediate-ranged structures.
Overall, the disagreement in diffusivities as calculated from our RIM-IPMD and AIMD simulations is caused due to the unrealistic and long-lived continuous network formation observed from the RIM-IPMD-generated salt structure ( Figures 8B,C). Such large discrepancies in the diffusion coefficient emphasize the importance of accurate intermediaterange structure prediction.

Thermophysical Properties
In this section, the densities from the three salt compositions as obtained from the AIMD and RIM calculations are reported and discussed. Both viscosity and heat capacity values are evaluated using the RIM parameters, and the values will be compared with the available experimental data.

Density
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 AIMDcalculated 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.

Viscosity
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: 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 2 [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 4 ), 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 4 [38], viscosity values tend to decrease with the decrease in the mol% of networkforming 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).

Specific Heat Capacity
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°i ntervals. 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][40][41][42][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: The specific heat capacities values of 0.256 ± 0.002 and 0.24 ± 0.006 cal/gm-°C were obtained for compositions A and

CONCLUSION
The ab-initio and classical molecular dynamics simulations were performed for 38-51-11 (composition A), 42-29-29 (composition B), and 26-37-37 (composition C) mol% of the LiF-NaF-ZrF 4 salt system in the molten state. The short-to intermediate-range salt structure generated by AIMD and RIM-IPMD were analyzed and compared using cation-anion coordination number, cation-anion and Zr-Zr RDF, angular distribution, and fluorozirconate chain length quantification. Densities obtained from RIM-IPMD were compared with AIMD, semi-empirical expression a.k.a. additive molar volume method, and experimental values. It was found that the AIMDobtained density for composition C is in reasonable agreement (within 11%) with the experimental values, while RIM-IPMDcalculated densities are off by 26%.
For all compositions considered in this study, RIM suggested the formation of 7, 8, and 9 coordinated fluorine zirconium complexes, while the AIMD calculations and previous Raman spectra indicate the occurrence of 6, 7, and 8 coordinated fluorozirconate complexes in the salt mix. Both RIM and AIMD simulations showed a decrease in the average Zr coordination number with the increase in ZrF 4 mol% going from compositions A to C. For all compositions, RIM reported up to 19% higher values for average Zr coordination numbers. The monovalent cation-anion (Li, Na) coordination numbers showed good agreement with those obtained from AIMD at lower ZrF 4 mol%, while their values deviated by 7% upon increase in ZrF 4 mol%. At higher ZrF 4 mol% (compositions B and C), there is a deviation between AIMD and RIM-calculated Zr-F coordination numbers. This indicates poor representation of the short-range interactions in RIM when a higher mol% of multivalent cation is present. At the intermediate-range, isolated zirconium fluorine complexes ([ZrF x 4−x ]) were dominant in both RIM and first-principle calculations when ZrF 4 mol% was limited (composition A). However, as ZrF 4 mol % increased such as in compositions B and C, both AIMD and RIM showed the presence of fluorozirconate chains ([ZrF x 4−x ] n ). Evidently, RIM parameters depicted an extensive network formation in the salt structure involving up to 99% of Zr in the melt compared to only 67% of Zr involved in chain formation in AIMD. Moreover, angle distribution functions revealed that the AIMD-generated network is connected by both edge-sharing and cornersharing fluorine-zirconium complexes, while the RIMgenerated network is mainly connected by corner-sharing complexes.
Such differences in coordination complexes ([ZrF x 4−x ]) resulted in~56% smaller Zr diffusion coefficients in composition A. An additional unrealistic network formation ([ZrF x 4−x ] n ) in higher ZrF 4 mol% compositions contributed to differences in AIMD and RIM Zr self-diffusion coefficients up to one order and two orders of magnitude for compositions B and C, respectively. Moreover, for compositions B and C, the overprediction of the network by RIM resulted in the viscosity values to be two and three orders higher than the reported experimental data. Where AIMD and RIM salt structures showed good agreement (composition A), reasonable Zr diffusivities and viscosity values calculated from RIM are predicted. It was also found that the specific heat capacities calculated from the RIM structure showed a decrease in values with the increase in ZrF 4 mol%, which agrees with the previously available experimental and simulation data. Upon comparison, the specific heat capacity values from RIM deviated by 27% from the available experimental data for composition B.
Overall, the discrepancies in coordination behavior can be attributed to the increase in inaccuracy of short-range interactions calculated by RIM potential as ZrF 4 mol% increases. Additionally, the observed discrepancies in the intermediate-range structure for LiF-NaF-ZrF 4 salt are likely caused by several factors including noninclusion of dispersion interaction in the parameter fitting process and omission of higher dipole moments while fitting for the polarization terms. Fluorine polarizability significantly affects the Zr-F-Zr angle and hence can influence the connectivity of the adjacent flurozirconates (corner-or edge-sharing).
While it can be argued that a re-parameterization of RIM parameters for the LiF-NaF-ZrF 4 salt system containing higher ZrF 4 mol% may lead to obtaining of an accurate short-to intermediate-range salt structure, the transferability of the refitted force-field parameters would be reduced. It is known that the parameter fitting for multicomponent molten salts containing multivalent charged species is a tedious process which involves multiple challenges. To this end, AIMD calculations are particularly required in interpolating and extrapolating across compositional space. In future studies, DFT-trained neural network interatomic potentials can also provide an alternative pathway for generating the accurate salt structure and evaluating structure-dependent transport as well as thermophysical properties.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be available by the authors, without undue reservation upon reasonable request.