Investigation of Microscopic Structure and Ion Dynamics in Liquid Li(Na, K) Eutectic Cl Systems by Molecular Dynamics Simulation

: Molten chloride salts are the main components in liquid metal batteries, high-temperature heat storage materials, heat transfer mediums, and metal electrolytes. In this paper, interest is centered on the inﬂuence of the LiCl component and temperature on the local structure and transport properties of the molten LiCl-NaCl-KCl system over the temperature range of 900 K to 1200 K. The liquid structure and properties have been studied across the full composition range by molecular dynamics (MD) simulation of a sufﬁcient length to collect reliable values, such as the partial radial distribution function, angular distribution functions, coordination numbers distribution, density, self-diffusion coefﬁcient, ionic conductivity, and shear viscosity. Densities obtained from simulations were underestimated by an average 5.7% of the experimental values. Shear viscosities and ionic conductivity were in good agreement with the experimental data. The association of all ion pairs (except for Li-Li and Cl-Cl) was weakened by an increasing LiCl concentration. Ion clusters were formed in liquids with increasing temperatures. The self-diffusion coefﬁcients and ionic conductivity showed positive dependences on both LiCl concentration and temperature, however, the shear viscosity was the opposite. By analyzing the hydrodynamic radii of each ion and the coordination stability of cation-anion pairs, it was speculated that ion clusters could be the cation-anion coordinated structure and affected the macro properties.


Introduction
Molten chloride salts have been widely used in liquid metal batteries [1][2][3], high-temperature heat storage materials, heat transfer mediums [4], and metal electrolytes, owing to their favorable characteristics, such as wide operating temperature range, low vapor pressure, moderate heat capacity, low cost, and high thermal stability [5]. Both the structure and properties of molten chloride salts are hard to measure due to the extreme high-temperature, corrosive, and volatile conditions, needless to consider the accuracy of testing. Fortunately, computer simulation methods have been developed rapidly in recent years, which could give deep insight into the microstructure and the properties of molten salts.
The mixtures of molten salts have been investigated intensively in the past years, and the composition dependence of properties has often been the central topic. Endoh and Okada [6] firstly studied the internal mobility of molten (Li, Na, K)Cl of its nearly eutectic composition at 773 K, 873 K, and 973 K by molecular dynamics (MD) simulation. Molecular dynamics simulation had been carried x units. The molecular dynamics simulations on mixtures composed of MCl 3 (where M 3+ is Y 3+ , Tb 3+ , La 3+ , and U 3+ ) in (LiCl-KCl) eut were conducted, and the correlation between the viscosity and structure was discussed in depth [11]. Manga et al. [12] investigated a molten ZnCl 2 -NaCl-KCl system by MD, and the interplay between the extent of the network structure, composition, and transport properties (viscosity, thermal conductivity, and diffusion) was characterized. The network modification due to the addition of K ions leads to the formation of non-bridging terminal Cl ions, and leads to a positive temperature dependence of thermal conductivity in these melts, a valuable property for heat transfer fluid design.
Most of the work concentrated on complex multi-valent molten salts, neglecting the essential study of molten alkali chlorides, which are the solvent or main components in many applications. The alkali chlorides ternary system has a much lower melting point compared to the corresponding pure ones. There are few systematic reports on the study of ternary alkali metal chlorides systems in a temperature range of 900-1200 K and no comprehensive analysis on the correlation between the local structure and properties. In this context, our group determined to systematically investigate the ternary systems. This paper aimed to study the LiCl effect on eutectic NaCl-KCl systems (molar ration NaCl:KCl = 0.506:0.494).
In this paper, molecular dynamics simulation was carried out to investigate the influence of LiCl concentration and temperature (900 K to 1200 K) on the local structure and properties of a molten LiCl-NaCl-KCl system. Correlations between the structure and properties were also analyzed. The paper is organized as follows: Section 1 is the Introduction. In Section 2 the simulation details and theoretical methods are introduced, as well as experimental details. In Section 3, we analyze and discuss the simulation results, with the available experimental results. Partial radial distribution function, coordination numbers distribution, and angular distribution functions are displayed to explore the local structure. Then, the density, shear viscosity, and ionic conductivity of ternary salts are listed and compared with experiments. Density and viscosity experimental data are used to evaluate the simulated results. The effects of the temperature and concentration on the micro-structure and properties are also discussed. In the end of this section, the correlation between the local structure and transport properties is intensively discussed. Lastly, conclusions are made in Section 4.

Simulations
Molecular dynamics simulations have been carried out on ternary systems composed of LiCl (where LiCl concentration ranged from 0% to 100% with a 20% interval) and (NaCl-KCl) eutectic (mole fraction ratio of NaCl:KCl is 0.506:0.494) in the temperature range from 900 to 1200 K. The simulation studies were based on the Born-Mayer-Huggins-Tosi-Fumi (BMHFT) rigid-ion (RI) interionic potential for the molten salts. The BMHFT force field reproduces the properties of molten salts well, and the formula is listed below: where the detailed implications of these symbols are presented in nomenclature and full exposition has been provided in a previous work [13]. It should be noted that a correction factor of 0.92 was multiplied by the original calculated ρ values in the chloride melts because the parameters given by Tosi and Fumi yield shorter distances by 10-30 pm for the first peak position of the cation-anion pair correlation functions in molten ternary chlorides [6,7]. The potential parameters used in this paper are listed in Table 1. Table 1. Potential parameters for different pairs in the ternary systems.

Ion Pairs Li-Li
Li-Na Li-K Li-Cl Na-Na Na-K Na-Cl K-K K-Cl Cl-Cl The open source software, Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [14], was chosen to conduct all the simulation. A cubic box with 2000 ions (1000 cation and 1000 anion) was adopted to perform the simulation. The ion numbers of simulated compositions with a 20% interval are listed in Table 2, and were mainly discussed throughout the paper to explore the composition and temperature effect. The periodic boundary conditions method (PBC) was used to eliminate the boundary effect and keep the particle numbers constant. The cutoff radius was r c = L/2 where L indexed the box length. For faster calculation, the particle-particle particle-mesh solver method was employed and the precision was set as 1 × 10 −4 . The Newton's equation of motion was solved by the Verlet algorithm, and the timestep during the simulation was ∆t = 1 fs. More details of the simulated set up can be found in our previous paper [13].   For initial isothermal-isobaric ensemble (NPT) simulations, the liquid densities could be determined with the following equation: where N is the number of different alkali chlorides, M is the molar mass of each salt, N A denotes the Avogadro constant, and v E is the equilibrated volume of the simulation cell at the given temperature in the NPT ensemble simulations. For the structure and transport properties' collection, a random velocity according to the Gaussian distribution was assigned to the ions at the desired temperature. Then, an equilibration process for the system was conducted in NPT ensemble under a pressure fixed at 0 GPa at the desired temperature for 3 ns. The production runs were performed in a canonical ensemble (NVT) with the Nose-Hoover thermostat method at the equilibrated cell volume. The damping parameters, which determined how rapidly the temperature and pressure were relaxed, were set to be 0.1 ps and 0.5 ps, respectively. The NVT stages lasted more than 10 ns for each system to ensure good statistics for the computation of the transport coefficients. During the production stages, the pair distribution functions, coordination numbers distribution, angular distribution, and self-diffusion coefficients were calculated.

Evaluated Properties
Theoretical methods of calculating the partial radial distribution function (RDF), angular distribution function, coordination distribution function, density, diffusion coefficient, and ionic conductivity were introduced in our previous reports [15][16][17]. The angular distribution function and coordination numbers distribution function were processed based upon the trajectory files by python.
The shear viscosity, η, is a measure of the resistance of a fluid to be deformed by shear stress. The reverse non-equilibrium molecular dynamics (RNEMD) method has been shown to be more accurate and efficient than equilibrium molecular dynamics simulation (EMD) and non-equilibrium molecular dynamics simulation (NEMD) [18][19][20]. Here, the shear viscosity was calculated through the RNEMD method based on linear response theory: where j z (p x ) is the momentum flux, which is the x component of the momentum transported in z direction per given time and per unit area. The term ∂v x ∂z , is the velocity gradient built up in z direction, which is also denoted as the shear rate. p x is the total momentum transferred in a simulation. In this method, momentum swaps were conducted between the middle and bottom bin of the simulation box, and the generated velocity gradient was measured. The viscosity was then calculated as the ratio of the total flux transferred and the velocity gradient as a follower:

Experiments
The density and viscosity of molten LiCl-NaCl-KCl were measured to test the simulation results. Three different materials were used in the experiments: Lithium chloride, sodium chloride, and potassium chloride. All the materials were analytical reagents and the main contents were ≥99.5% (Sinopharm Chemical reagent Co. Ltd., Shanghai, China). All the mixture samples were prepared in the glove box (Mikrouna (China) Co., Ltd., Shanghai, China), then pre-melted in a muffle furnace (Y-feng furnace Co., Ltd., Shanghai, China) at the desired temperature, which is higher than the melting point. The pre-melting process was for at least one hour to ensure the sample mixing was uniform. After cooling in the drier, the samples were ground into a powder for testing. A self-designed device for molten salts based on the Archimedes method was used for density experiments. The density self-designed device was constructed following Cheng et al.'s work [21], and the set-up is illustrated in Figure 1. The Anton Paar Rheology EC-twist502 (Anton Paar GmbH, Graz, Austria) and furnace CTD1000 (Anton Paar GmbH, Graz, Austria) were used to measure the viscosity by the coaxial cylinders method. The shear rate was selected to be 20 s −1 after shear rate scanning.

Results
In this section, the addition of LiCl to a eutectic proportion (Na, K) Cl system has been investigated across the full composition range from 900 to 1200 K with a 25 K temperature interval. Furthermore, a higher temperature of 1500 K for the extreme condition was also investigated.

Partial Radial Distribution Function
(a) The LiCl concentration dependence Partial radial distribution function (PDF) is calculated in histogram form by binning pairwise distances into 200 bins from 0.0 to the maximum force cutoff (defined as rc = L/2). The bins are of uniform size in a radial distance. Thus, a single bin encompasses a thin shell of distances in 3d and a thin ring of distances in 2d. PDFs of A-2, A-4, A-6, and A-8 at 1100 K are shown in Figure 2 to illustrate the influence of the LiCl concentration on the structure of ternary systems. The quantitatively structural information containing the first peak position, rmax, the first peak height, hmax, the first valley position, rmin, and running integral coordination number, N(rmin) are summarized in Table 3, which were determined by fitting PDFs (R 2 = 1.0). The running integral coordination number was obtained by a spherical integration of g(r) from 0 to rmin. The structural parameters of the ternary systems were in good agreement with results of a similar composition [22].
For cation-anion pairs, though the variations of the first peak positions of potassium are a little greater than lithium sodium ions across the concentration, the changes (maximum change was 0.05 Å) were subtle. However, the first peak heights all decreased distinctively, while the heights of the first minimum increased with an increasing LiCl. This indicated the degree of unlike charged ions' association decreased with LiCl addition. The coordination numbers were roughly 4, 5, and 6 for Li + , Na + , and K + , respectively. With an increasing LiCl molar fraction, the coordination numbers of the cation-anion pairs all increased, which suggested that cations tended to accommodate more anions.
For cation-cation pairs, as shown in Figure 2, the PDF curves of the cation-cation pairs are not as smooth as the ones of cation-anion pairs for the relatively weak association between like charge ions. The first peak positions shifted right and the corresponding height decreased with LiCl addition (except for Li-Li). For Li-Li, rmax showed an overall increasing tendency, but hmax increased up to 60% of the LiCl concentration then decreased. The LiCl concentration had more influence on the structure of cation-cation pairs than the cation-anion pairs. The changing rates of rmax and hmax for Li-Na and Li-K pairs were smaller than those for Na-Na, Na-K, and K-K pairs, which indicated a larger cation radius difference was a more stable structure even in cation-cation pairs.

Results
In this section, the addition of LiCl to a eutectic proportion (Na, K) Cl system has been investigated across the full composition range from 900 to 1200 K with a 25 K temperature interval. Furthermore, a higher temperature of 1500 K for the extreme condition was also investigated.  Figure 2 to illustrate the influence of the LiCl concentration on the structure of ternary systems. The quantitatively structural information containing the first peak position, r max , the first peak height, h max , the first valley position, r min , and running integral coordination number, N(r min ) are summarized in Table 3, which were determined by fitting PDFs (R 2 = 1.0). The running integral coordination number was obtained by a spherical integration of g(r) from 0 to r min . The structural parameters of the ternary systems were in good agreement with results of a similar composition [22].
For cation-anion pairs, though the variations of the first peak positions of potassium are a little greater than lithium sodium ions across the concentration, the changes (maximum change was 0.05 Å) were subtle. However, the first peak heights all decreased distinctively, while the heights of the first minimum increased with an increasing LiCl. This indicated the degree of unlike charged ions' association decreased with LiCl addition. The coordination numbers were roughly 4, 5, and 6 for Li + , Na + , and K + , respectively. With an increasing LiCl molar fraction, the coordination numbers of the cation-anion pairs all increased, which suggested that cations tended to accommodate more anions.
For cation-cation pairs, as shown in Figure 2, the PDF curves of the cation-cation pairs are not as smooth as the ones of cation-anion pairs for the relatively weak association between like charge ions. The first peak positions shifted right and the corresponding height decreased with LiCl addition (except for Li-Li). For Li-Li, r max showed an overall increasing tendency, but h max increased up to 60% of the LiCl concentration then decreased. The LiCl concentration had more influence on the structure of cation-cation pairs than the cation-anion pairs. The changing rates of r max and h max for Li-Na and Li-K pairs were smaller than those for Na-Na, Na-K, and K-K pairs, which indicated a larger cation radius difference was a more stable structure even in cation-cation pairs. For anion-anion pairs, the first peak positions shifted from 4.22 Å to 3.75 Å across the concentration, while the first peak heights rose from 1.73 to 2.04. The equilibrium distance of Cl-Cl was shortened and the association between Cl-Cl was strengthened. The coordination number of Cl-Cl decreased as LiCl increased, which indicated the closed pack between anions was inhibited. This could be attributed to the significant anion-anion overlap in LiCl salts relative to the other alkali chlorides [23]. For anion-anion pairs, the first peak positions shifted from 4.22 Å to 3.75 Å across the concentration, while the first peak heights rose from 1.73 to 2.04. The equilibrium distance of Cl-Cl was shortened and the association between Cl-Cl was strengthened. The coordination number of Cl-Cl decreased as LiCl increased, which indicated the closed pack between anions was inhibited. This could be attributed to the significant anion-anion overlap in LiCl salts relative to the other alkali chlorides [23].   Table 3. Structural parameters for six typical ternary systems at 1100 K.   temperature range of 900 K to 1200 K. The increasing temperature led to a lower h max and an increase of the minima offset, which indicated the looser coordination structure.

K xLiCl Li-Cl Na-Cl K-Cl Li-Li Li-Na Li-K Na-Na Na-K K-K Cl-Cl rmax/Å
However, for 1500 K, the r max of cation-cation pairs obviously shifted to the right, while the second peak became weaker, which revealed the structure of ions became more liquid-like and disordered. However, the r max of cation-anion pairs unexpectedly shifted left. The equilibrium distances of cation-cation pairs increased, however, the cation-anion was the opposite, which means ion clusters were formed. The increasing temperature caused free space in melts to become larger, thus cations drifted apart. Meanwhile, the distances of cations and anions became shorter, thus the ion clusters were likely to be cation-anion clusters. Another possibility for this phenomenon could be that the free space became larger in melts with an ascending temperature, so less correlated, similarly charged ion pairs drifted apart, and only the more tightly coordinated anions around the cations were left.
The temperature influence on Cl-Cl pairs was similar with the cation-cation pairs. At 1500 K, the r max of Cl-Cl pairs shifted right slightly and the corresponding h max became weaker. However, it is interesting to note that with the LiCl concentration increasing, the Cl-Cl pairs were less affected by the temperature. We assumed that this was caused by the stronger interactions between Li-Cl pairs.
In conclusion, though the ternary system is chemically simple, there are complicated changes in the association of each ion with its first coordination shell with the composition changing. The addition of LiCl to molten ternary systems decreases the association between ion pairs, except for Li-Li and Cl-Cl pairs. It is worthy to note that the addition of LiCl strengthened the Cl-Cl association, and the closed pack between anions was inhibited. Ion clusters were formed in liquid and the distance of ion clusters became larger with an increasing temperature, which has great influence on the macro-properties. To further study the local structure, PDFs must be combined with the coordination numbers distribution and angular distribution function.

Coordination Numbers Distribution
The simulation enabled deep insight into the numbers of anions in the coordination shell. We determined the number of anions within a spherical shell by collection of every step ions' position and the shell radius was set to the position of the minimum (r min ) in the PDFs. This constitutes an instantaneous coordination number for that cation. Figure 3 shows the coordination numbers distribution around each cation at 1100 K across the concentration. In the ternary melts, there coexisted 3-to 8-fold coordination, with different proportions. The coordination numbers distribution of Li + are mainly located in 4 and coexisted with 3-fold and 5-fold coordination. With an increasing LiCl composition, the mainly located 4-fold coordination of Li + slightly decreased. Meanwhile, the proportion of 3-fold distribution decreased, and the proportion of 5-fold coordination number increased. For Na-Cl, the coordination numbers, mainly located in 5, coexisted with 4-fold and 6-fold coordination. The 5-fold coordination almost remained constant, with no regard to the composition. However, with an increasing LiCl concentration, the 4-fold coordination largely decreased and 6-fold coordination increased. For K-Cl, the coordination numbers of K + were mainly 4, 5, 6, 7, and 8. When x LiCl ≤ 60%, the coordination numbers were mainly located in 6. However, when x LiCl > 60%, the 6-fold and 7-fold dominated. With an increasing LiCl, the 4-fold and 5-fold coordination of K + all decreased, while the higher coordination 7-fold and 8-fold increased.  The results showed that the proportion of higher coordination numbers for each cation increased as the LiCl molar fraction increased. Combined with our previous study [13,15], the local structure of Li-Cl in melts is an anomalous, distorted tetrahedron, while Na-Cl and K-Cl are a distorted octahedron. It can be clearly seen that the tendency to form an octahedron became stronger for Na-Cl with an increasing LiCl addition. The anions coordination numbers distribution of K + were the most complex. Though sexadentate dominated, the tendency of accommodating more anions for K + became stronger. As for Li-Cl, though the portion of 5-fold coordination increased, the tetrahedron structure still dominated. The sharp peak in the PDF indicates a tight coordination of Li-Cl and its height decreases with LiCl addition. It can be found from Table 3 that when the equilibrium distance between the cation and anion gets longer, the coordination numbers of the cation become greater. Combined with PDFs analysis, the LiCl addition increased the mobility of Cl -, so the tighter coordination structure (the relatively lower coordinated structure) reduced. With the temperature increasing, the coordination became amorphous (the temperature influence on the coordination numbers distribution is provided in Appendix A Figures A5-A7). In Figure 4, there is a snapshot of the A-2 system at 1100 K, and the coordinates around different cations were in a dynamic equilibration. The results showed that the proportion of higher coordination numbers for each cation increased as the LiCl molar fraction increased. Combined with our previous study [13,15], the local structure of Li-Cl in melts is an anomalous, distorted tetrahedron, while Na-Cl and K-Cl are a distorted octahedron. It can be clearly seen that the tendency to form an octahedron became stronger for Na-Cl with an increasing LiCl addition. The anions coordination numbers distribution of K + were the most complex. Though sexadentate dominated, the tendency of accommodating more anions for K + became stronger. As for Li-Cl, though the portion of 5-fold coordination increased, the tetrahedron structure still dominated. The sharp peak in the PDF indicates a tight coordination of Li-Cl and its height decreases with LiCl addition. It can be found from Table 3 that when the equilibrium distance between the cation and anion gets longer, the coordination numbers of the cation become greater. Combined with PDFs analysis, the LiCl addition increased the mobility of Cl -, so the tighter coordination structure (the relatively lower coordinated structure) reduced. With the temperature increasing, the coordination became amorphous (the temperature influence on the coordination numbers distribution is provided in Appendix A Figures A5-A7). In Figure 4, there is a snapshot of the A-2 system at 1100 K, and the coordinates around different cations were in a dynamic equilibration.
Combined with PDFs analysis, the LiCl addition increased the mobility of Cl -, so the tighter coordination structure (the relatively lower coordinated structure) reduced. With the temperature increasing, the coordination became amorphous (the temperature influence on the coordination numbers distribution is provided in Appendix A Figures A5-A7). In Figure 4, there is a snapshot of the A-2 system at 1100 K, and the coordinates around different cations were in a dynamic equilibration.

Angular Distribution Function
The angular distribution function (ADF) of Cl-A-Cl (A = Li, Na, K) from A-2 to A-8 at 1100 K is shown in Figure 5. Here, "bond" does not mean a real chemical bond, but just an atom pair within some distance. In the ternary melts, the Cl-A-Cl bond angles were distributed between 60 • and 150 • , all showing a broad peak and a smearing phenomenon, which indicates that the coordination bonds are orientated randomly.

Angular Distribution Function
The angular distribution function (ADF) of Cl-A-Cl (A = Li, Na, K) from A-2 to A-8 at 1100 K is shown in Figure 5. Here, "bond" does not mean a real chemical bond, but just an atom pair within some distance. In the ternary melts, the Cl-A-Cl bond angles were distributed between 60° and 150°, all showing a broad peak and a smearing phenomenon, which indicates that the coordination bonds are orientated randomly.
With an increasing LiCl addition, the peak values of the bond angle for Cl-Li-Cl shifted from 100º to 93°, for Cl-Na-Cl, shifted from 88° to 79°, and for Cl-K-Cl, shifted from 78° to 68°. The angles become smaller with an increasing LiCl concentration. The corresponding bond angles for pure melts of LiCl, NaCl, and KCl were 98.18°, 87.27°, and 83.64°, respectively [15]. The bond angles of the corresponding pure alkali chloride crystal were mainly distributed in 90° and 180° [24] because the three salts in the crystal type were a face-centered cubic structure and sexadentate. Though the longrange order disappears in melts, at the micro-level, there was still some tendency to keep the crystal structure. Considering the coordination numbers of each cation become larger and the bond angles become smaller, we assumed that the deforming in the micro-level structure was larger than in that of the pure [15] and binary melt systems [25]. The ADFs were broader and no peak shifting was observed as the temperature increased (the temperature influence on the ADFs is provided in Appendix A Figure A8).

Density
The comparison of simulated and experimental densities for A-2, A-4, A-6, and A-8 systems are plotted in Figure 6. The experimental results were 6.3%, 8.1%, 3.8%, and 4.5% larger than the simulated results from A-2 to A-8. The discrepancies can be contributed to the selected force field and experimental error. Despite many efforts to keep the samples pure before testing, the LiCl should be hydrated and this affected the results. The BMHFT potential underestimated the density by about 5.3%, 6.3%, and 9.0% for pure molten LiCl, NaCl, and KCl, respectively [26]. Thus, it seems that the BMH pair potential, Larsen combination rule, and correction factors do a reasonable job of capturing molten phase properties. With an increasing LiCl addition, the peak values of the bond angle for Cl-Li-Cl shifted from 100º to 93 • , for Cl-Na-Cl, shifted from 88 • to 79 • , and for Cl-K-Cl, shifted from 78 • to 68 • . The angles become smaller with an increasing LiCl concentration. The corresponding bond angles for pure melts of LiCl, NaCl, and KCl were 98.18 • , 87.27 • , and 83.64 • , respectively [15]. The bond angles of the corresponding pure alkali chloride crystal were mainly distributed in 90 • and 180 • [24] because the three salts in the crystal type were a face-centered cubic structure and sexadentate. Though the long-range order disappears in melts, at the micro-level, there was still some tendency to keep the crystal structure. Considering the coordination numbers of each cation become larger and the bond angles become smaller, we assumed that the deforming in the micro-level structure was larger than in that of the pure [15] and binary melt systems [25]. The ADFs were broader and no peak shifting was observed as the temperature increased (the temperature influence on the ADFs is provided in Appendix A Figure A8).

Density
The comparison of simulated and experimental densities for A-2, A-4, A-6, and A-8 systems are plotted in Figure 6. The experimental results were 6.3%, 8.1%, 3.8%, and 4.5% larger than the simulated results from A-2 to A-8. The discrepancies can be contributed to the selected force field and experimental error. Despite many efforts to keep the samples pure before testing, the LiCl should be hydrated and this affected the results. The BMHFT potential underestimated the density by about 5.3%, 6.3%, and 9.0% for pure molten LiCl, NaCl, and KCl, respectively [26]. Thus, it seems that the BMH pair potential, Larsen combination rule, and correction factors do a reasonable job of capturing molten phase properties. Densities across the whole composition and temperature range are shown in Figure 7. It decreased linearly with the ascending temperature and Li + addition. A database has been built for the ternary system across the full concentration, and the fitting expression is shown below: Densities across the whole composition and temperature range are shown in Figure 7. It decreased linearly with the ascending temperature and Li + addition. A database has been built for the ternary system across the full concentration, and the fitting expression is shown below: where x is the concentration of LiCl. Neglecting the minor terms and fitting again: ρ density = −0.00172x − 6.58567 × 10 −4 T − 1.34192 × 10 −6 x 2 + 4.29788 × 10 −8 T 2 + 1.33476 × 10 −6 xT + 2.08316 (6) = -0.00164x -8.14236 × 10 -4 T -1.51351 × 10 -6 x 2 + 1.94672 × 10 -7 T 2 -1.94938 × 10 -9 x 3 -4.86205 × 10 -11 T 3 +1.24708 × 10 -6 xT + 4.17543×10 -13 x 2 T 2 + 2.13573 (5) where x is the concentration of LiCl. Neglecting the minor terms and fitting again: = -0.00172x -6.58567 × 10 -4 T -1.34192 × 10 -6 x 2 + 4.29788 × 10 -8 T 2 + 1.33476 × 10 -6 xT + 2.08316 (6) Figure 7. Surface map of simulated density results.
As listed in Table 4, the absolute values of the slopes of fitted lines become smaller with the addition of lithium ions, which indicated that the slope decreased slower with ascending temperature as the LiCl addition increased. Derived from Wang et al.'s [15] work, pure LiCl has the smallest slope in LiCl, NaCl, and KCl. Therefore, the higher LiCl concentration was, the more slowly the density changed along the temperature. As listed in Table 4, the absolute values of the slopes of fitted lines become smaller with the addition of lithium ions, which indicated that the slope decreased slower with ascending temperature as the LiCl addition increased. Derived from Wang et al.'s [15] work, pure LiCl has the smallest slope in LiCl, NaCl, and KCl. Therefore, the higher LiCl concentration was, the more slowly the density changed along the temperature.

Self-Diffusion Coefficient
The self-diffusion coefficient is an important parameter in electrokinetics, metal extraction, and pyrochemical processing of nuclear waste; it is important to understand how the diffusion coefficient value is affected by changes in the solvent, and how it depends on the ionic charge and size.
The self-diffusion coefficients were calculated from the slope of the mean square displacement in the linear region for each ion in the melt. The simulated data are presented in Table 5. The results were in good agreement with corresponding experiments' data of pure salts [27,28]. The values fluctuated as the component changed, and showed an overall positive dependence on the LiCl concentration and temperature. The fluctuation may be attributed for two reasons: First, according to Fürth hole theory [29], holes were formed by the action of the irregular thermal movement. Larger holes increased as the temperature increased. Second, ion clusters may form and the very short existence of ion clusters influences the mobility of ions. It could be observed from Table 5 that ions in melts were less diffused in ternary systems than in the corresponding pure melt [28,30] and the binary systems (LiCl-KCl [30] and LiCl-NaCl [16]). The Li ions' diffusivity was smaller than the Na ions' and K ions' at a low LiCl concentration and exceeded them at a higher concentration, which indicated the chemla effect. The diffusion coefficient of each ion exhibited an overall increase tendency as the LiCl concentration increased, which was consistent with the h max decrease in PDFs. The increasing rate of D(Li + ) become greater when the molar fraction of LiCl was beyond 60%. The anions' mobility incensement was smaller than the cations'.

Ionic Conductivity
The ionic conductivity was computed based on the Green-Kubo correlation [31] via a programme compiled by language C in LAMMPS [13,15]. The ionic conductivity, λ z , can be calculated from the time integral of the charge flux autocorrelation function through the EMD method: where the charge flux vector, J Z (t), is defined as: where z i e and v i are the charge and velocity of atom i, respectively. Each of the three independent components (i.e., J xz , J yz , and J zz ) of the charge flux vector provides an independent estimate of the ionic conductivity, and the averaged value is taken as the final ionic conductivity value. Figure 8a shows the time evolution of the normalized charge flux auto-correlation function (ACF) and the running integral (integration of charge flux ACF) for a molten A-6 system at 900 K. The normalized ACFs all decayed to zero quickly. The running integrals reached a plateau after about 0.5 ps for all the systems. During ionic conductivity calculations, a time of 5 ps was adopted for the charge flux ACFs calculation. Figure 8b shows the simulated ionic conductivity results of ternary melts and the available experimental data. The ionic conductivities obtained here were in good agreement with the experiment, especially for molten LiCl. Matiasovsky and coworkers [32] investigated the specific electrical conductivity of LiCl-NaCl-KCl experimentally, and the results are plotted in Figure 8b (marked in line and symbols). In the experimental temperature range, the results of the calculation agreed with the experimental data within the error of 8-12%. Further research, including the polarization effect, should improve the accuracy of the calculated results. Besides, take A-6 system as an example, the Nernst-Einstein fitting [33,34], which neglects ionic correlations, was used to validate the simulation results here. The Nernst-Einstein relation, unlike the Einstein-Helfand relation, which is equivalent to the Green-Kubo relation at a long-time limit, is missing the cross term: D + D - [35]. Most of the ionic conductivity fitting for molten salts used the Nernst-Einstein formula, and the formula used here [33] is listed below: Appl

Shear Viscosity
Viscosities were calculated by the RNEMD method. The systems were divided into 20 slabs in the z axis, and the velocities in the "x" direction were swapped between the 1st and 11th bins [19,20,25]. Velocity profiles were equilibrated for 1 ns. Next, the momentum transferred and velocity gradient were accumulated over 10 ns to ensure the average of the velocity equilibrium. Once the system relaxed, the velocity gradient in the "x" direction was almost linear.
Different swap rates' effect on the linearity of velocity profiles, z direction lengths' effect on the viscosities, and the temperature dependence of velocity profiles were investigated, and the simulation conducted with optimized RNEMD parameters. It was found that momentum swap rates between 20 and 150 time steps showed good linearity velocity gradients, as illustrated in Figure 9a, and 100 was chosen. The larger ones caused a large deviation from the experimental values, as shown in Figure 9b. The z direction box length effect on the results was tested at 900 K. In Figure 10a, all the velocity profile linearity is acceptable (a slightly nonlinear of the velocity gradients will not affect the accuracy of the calculated viscosity values [19]). In Figure 10b, the results showed no obvious change beyond the 2z length in the z direction, and, taking the computation efficiency into account, 2z was picked. The temperature dependence of the velocity gradient was testified, and there is a minor effect of temperature on the linearity of the velocity gradient, as illustrated in Figure 11. In Figure 8b, when the LiCl addition was less than 60%, the ionic conductivity increased slightly. When LiCl addition was more than 60%, the ionic conductivity increased rapidly. Therefore, in practical industrial applications, it should be noted that the addition of LiCl largely increased ionic conductivity beyond a certain concentration. Furthermore, it also showed that the ionic conductivities of molten ACl increased with the ascending temperature.

Shear Viscosity
Viscosities were calculated by the RNEMD method. The systems were divided into 20 slabs in the z axis, and the velocities in the "x" direction were swapped between the 1st and 11th bins [19,20,25]. Velocity profiles were equilibrated for 1 ns. Next, the momentum transferred and velocity gradient were accumulated over 10 ns to ensure the average of the velocity equilibrium. Once the system relaxed, the velocity gradient in the "x" direction was almost linear.
Different swap rates' effect on the linearity of velocity profiles, z direction lengths' effect on the viscosities, and the temperature dependence of velocity profiles were investigated, and the simulation conducted with optimized RNEMD parameters. It was found that momentum swap rates between 20 and 150 time steps showed good linearity velocity gradients, as illustrated in Figure 9a, and 100 was chosen. The larger ones caused a large deviation from the experimental values, as shown in Figure 9b. The z direction box length effect on the results was tested at 900 K. In Figure 10a, all the velocity profile linearity is acceptable (a slightly nonlinear of the velocity gradients will not affect the accuracy of the calculated viscosity values [19]). In Figure 10b, the results showed no obvious change beyond the 2z length in the z direction, and, taking the computation efficiency into account, 2z was picked. The temperature dependence of the velocity gradient was testified, and there is a minor effect of temperature on the linearity of the velocity gradient, as illustrated in Figure 11.

Shear Viscosity
Viscosities were calculated by the RNEMD method. The systems were divided into 20 slabs in the z axis, and the velocities in the "x" direction were swapped between the 1st and 11th bins [19,20,25]. Velocity profiles were equilibrated for 1 ns. Next, the momentum transferred and velocity gradient were accumulated over 10 ns to ensure the average of the velocity equilibrium. Once the system relaxed, the velocity gradient in the "x" direction was almost linear.
Different swap rates' effect on the linearity of velocity profiles, z direction lengths' effect on the viscosities, and the temperature dependence of velocity profiles were investigated, and the simulation conducted with optimized RNEMD parameters. It was found that momentum swap rates between 20 and 150 time steps showed good linearity velocity gradients, as illustrated in Figure 9a, and 100 was chosen. The larger ones caused a large deviation from the experimental values, as shown in Figure 9b. The z direction box length effect on the results was tested at 900 K. In Figure 10a, all the velocity profile linearity is acceptable (a slightly nonlinear of the velocity gradients will not affect the accuracy of the calculated viscosity values [19]). In Figure 10b, the results showed no obvious change beyond the 2z length in the z direction, and, taking the computation efficiency into account, 2z was picked. The temperature dependence of the velocity gradient was testified, and there is a minor effect of temperature on the linearity of the velocity gradient, as illustrated in Figure 11.   The viscosities obtained from experiments and simulations are plotted in Figure 12. Viscosities calculated by the RNEMD method are in good agreement with the experimental values. The uncertainties of simulated viscosities for A-0, A-2, A-6, and A-8 were under 7%. Since the RNEMD method underestimated viscosities by about 5.8%, 5.2%, and 5.7% for pure molten LiCl, NaCl, and KCl, respectively [26], the simulation results were satisfying. The relatively large discrepancies of A-4 and A-10 (about 20% discrepancy) could be caused by experimental errors, since the A-4 sample corroded the crucible badly and there were slags in the salts during testing. The viscosity monotonously decreased with the increasing temperature. The uncertainties of simulated viscosities for A-0, A-2, A-6, and A-8 were under 7%. Since the RNEMD method underestimated viscosities by about 5.8%, 5.2%, and 5.7% for pure molten LiCl, NaCl, and KCl, respectively [26], the simulation results were satisfying. The relatively large discrepancies of A-4 and A-10 (about 20% discrepancy) could be caused by experimental errors, since the A-4 sample corroded the crucible badly and there were slags in the salts during testing. The viscosity monotonously decreased with the increasing temperature.

The Correlation between Local Structures and Properties
The hydrodynamic radius was calculated from the correlation between the self-diffusion coefficient and shear-viscosity of ionic liquids, which was an effective radius that takes solvation and coordination into account. It can be calculated according to the Stokes-Einstein correlation [39]: The results are listed in Table 6 (the crystal radii is from Ref. [40]). The hydrodynamic radius of

The Correlation between Local Structures and Properties
The hydrodynamic radius was calculated from the correlation between the self-diffusion coefficient and shear-viscosity of ionic liquids, which was an effective radius that takes solvation and coordination into account. It can be calculated according to the Stokes-Einstein correlation [39]: The results are listed in Table 6 (the crystal radii is from Ref. [40]). The hydrodynamic radius of each cation decreased significantly with LiCl addition, which indicated the hydrodynamic radius was influenced by the concentration. It seems the hydrodynamic radii of cations appeared to decrease when the imagined strength of the correlation of the cation with its coordination shell of anions weakened. With an ascending temperature, the variation in the hydrodynamic radius of each ion showed an overall increase. It was found that for the larger, singly charged K + and Clions, which could be anticipated to be most weakly coordinated to their surroundings, the calculated radii were similar to their crystal radii. For the remaining ions, the calculated radii were significantly larger than the bare ions under a 60% LiCl concentration, and became similar beyond 60%. At an 80% LiCl concentration, the hydrodynamic radius for Li + was still larger than for bare ions, but for Na + , it was smaller than the crystal radius under 1150 K. For the analysis, we speculated that the ions were diffusing with an intact coordination shell of the central cation, which was in agreement with the hypothesis of Morgan [34] that Li + moves with the surrounding Cl -. Additionally, some of the coordinates were broken apart by LiCl increasing.
Next, we attempted to identify the strength of the A-Cl ionic "bonds". Many microscopic events can contribute to the local structure, and the breaking apart of the A + coordination shells should be involved. The stability of the 'bond' could be estimated using the potential of the mean force that determines the exchange of Cl − with the first solvation shell of A + , and can be evaluated by [11]: where β = 1/k B T, and r min is the minimum of RDF as a constant. As shown in Figure 13, deep first minimums appeared at the position of the corresponding first peak position of the RDF. The height between the minimum and maximum was called the energy barrier (E barrier ), which was the energy that a coordinated Cl − anion must overcome to escape. The barrier to dissociation is higher for the smaller cation, which indicates the stability as follows: Li-Cl > Na-Cl > K-Cl. The heights almost all decreased with an increasing concentration of LiCl, which is reflected in the higher self-diffusion coefficients and ionic conductivity, and smaller viscosity. As shown in Figure 13, deep first minimums appeared at the position of the corresponding first peak position of the RDF. The height between the minimum and maximum was called the energy barrier (Ebarrier), which was the energy that a coordinated Cl − anion must overcome to escape. The barrier to dissociation is higher for the smaller cation, which indicates the stability as follows: Li-Cl > Na-Cl > K-Cl. The heights almost all decreased with an increasing concentration of LiCl, which is reflected in the higher self-diffusion coefficients and ionic conductivity, and smaller viscosity.

Conclusions
In this paper, comprehensive molecular dynamics simulations have been carried out to calculate local structure and transport properties of molten alkali chlorides ternary systems. Densities of the liquid phases obtained from simulations were underestimated by 5.7% of the experimental values. Shear viscosities calculated by the RNEMD method, as well as the ionic conductivity, were in good agreement with experimental data. The BMHFT force field combined with the Larsen combination rule and the correction factor used in this study did a reasonable job of capturing the molten phase properties of the ternary salts.
The liquid structures were fully calculated to analyze the influence of the temperature and component on the local structures of molten salts. The result suggested that ion clusters were formed in liquid as temperature increased, which has an influence on the macro-properties. The association of all ion pairs (except for Li-Li and Cl-Cl) was weakened by an increasing LiCl concentration. Therefore, the self-diffusion coefficients and ionic conductivity showed positive dependences on both the LiCl concentration and temperature, however, the shear viscosity was the opposite. By analyzing the hydrodynamic radius of each ion and the coordination stability of cations, it was speculated that ion clusters could be the cation-anion coordinated structure and affected the macro properties. In the system, the properties changed significantly beyond a 60% LiCl concentration, which should be noted in applications.  Figure A2. PDFs of unlike cation-cation pairs in the four ternary systems (A-2, A-4, A-6, and A-8) at 900 K, 1000 K, 1100 K, 1200 K, and 1500 K.