First-principles AIMD simulations were performed using the Vienna Ab-Initio Simulation Package (VASP) based on density functional theory (DFT) and the Born–Oppenheimer approximation with periodic boundary conditions (Vosko et al., 1980; Kresse and Hafner, 1993; Kresse and Hafner, 1994; Kresse and Furthmüller, 1996a; Kresse and Furthmüller, 1996b). The interactions between electrons and nucleus are defined by the projector augmented wave (PAW) method, and the revised Perdew–Burke–Ernzerhof (rPBE) DFT of the generalized gradient approximation (GGA) is used for the electron exchange-correlation. The kinetic energy cutoff is 400 eV, and the system has a 1 × 1 × 1 k-point mesh. The timestep chosen is two femtoseconds (fs) to avoid energy drift. Fermi smearing with a smearing parameter of 0.2 eV was used for the partial occupancies of the wave function (Grimme, 2006a; Grimme et al., 2011; Hacene, et al., 2012; Hutchinson and Widom, 2012).

In this research, the molten salt system under investigation is a ternary chloride salt MgCl₂–NaCl–KCl, 44.8 mol% MgCl₂–29.4 mol% NaCl–25.8 mol% KCl, provided by NREL. Based on the experimental compositions, two systems were generated and used in this study: 142-atom (26 Mg₂+, 17 Na+, 15 K+, 84 Cl-, 58 cations and 84 anions) (Grimme et al., 2010) and 83-atom (15 Mg₂+, 10 Na+, 9 K+, 49 Cl-, 34 cations and 49 anions). The following were regarded as the valence electrons: Mg₂+ 3s2, Na+ 2s22p63s1, K+ 3s23p64s1, and Cl 3s23p5. The 83-atom system was used for the initial evaluation of the dispersion forces. Both systems were used to predict the structure and properties of the ternary salt.

The starting configuration file of the two systems was generated using PACKMOL, which calculated the initial geometries of the system by randomly packing atoms into a given volume based on the experimental composition and density at a specified temperature (Martínez et al., 2009). The cell sizes for the 83-atom system ranged 13.75–14.20 Å; the cell sizes for the 142-atom system ranged 16–17 Å. Before these initial configuration files could be used for AIMD simulations, the system needed to be pre-equilibrated with classical interatomic potential molecular dynamic (IPMD) simulations. These were performed with a Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) and utilized the Born–Mayer–Huggins potentials for each component of our system (Mayer, 1933; Fumi and Tosi, 1964; Tosi and Fumi, 1964; Thompson, et al., 2022). Even though there were limitations with the force field potentials for complex ternary salts, the accuracy of the pre-equilibration stage was not essential. The pre-equilibration IPMD simulation used an NVT ensemble at the specified temperature for 5 ps to lose the memory of the initial configuration from PACKMOL (Martínez et al., 2009). This method of pre-equilibration using IPMD followed Bengtson et al. (2014) and Nam, et al. (2014) The final configuration of the IPMD simulation was used for the following AIMD simulations. With the final configuration of the system from the IPMD simulation, the configuration underwent another step of equilibration using an isobaric–isothermal (NPT) ensemble for total simulation time (Steinmann and Corminboeuf, 2011). A generalized-gradient approximation exchange hole model was used for dispersion coefficients of approximately 100 ps at the respective temperature with a timestep of 1 fs and pressure of 1 atm. The NPT ensemble used the Langevin thermostat with a temperature coefficient set to 10 ps⁻¹. The goal of the NPT equilibration was to evaluate the density and energy of the system. The configuration used for the next simulation was selected from the trajectory of the NPT equilibration simulation. This configuration was selected from a timestep where the density of the system was approximately its average density from the overall simulation. With this selected configuration file, the system underwent an NVT ensemble using a Nosè thermostat for 100 ps with a timestep of 2fs, referred to as the “production run”. The production runs from each temperature were used to estimate the thermodynamic properties at that temperature. The first 5 ps from the production run were neglected from the analysis and served as further equilibration. The simulated temperatures for both the 83-atom and 142-atom systems are 723 K, 773 K, 823 K, 873 K, 923 K, and 973 K.

Properties were predicted using the trajectories from the production runs at the respective temperature. Each trajectory had a total simulation time of approximately 100 ps, used an NVT ensemble, and then was analyzed using MDANSE (Goret, 2017). The results given in Section 3.1 were simulated with an NPT ensemble while those in Section 3.2 were simulated with an NVT ensemble—the production runs.

The density of the salt in the molten phase was measured with a density meter employing Archimedes’ principle of buoyancy. This meter utilizes the weight of an object, such as high purity nickel cylinder, and a quartz container, which holds the molten salt. The measurement relies on measuring the object’s mass before and after submerging in the molten salt. The scale used had a full capacity of 50 g and an uncertainty of ±0.1%. Utilizing the disparity in the mass between the object in air and submerged in molten salt, the density of the molten salt across various temperatures is calculated using Eq. 7: ρ f = ∆ M M ρ M (7) where ∆ M is the difference in the weight of the object before and after submerging caused by the buoyancy force of the molten salt, M is the weight of the object measured in air, ρ M is the density of the object, and ρ f is the density of the molten salt. This method has been utilized in the literature with good reliability. The systematic error of the experimental density can be calculated using the partial derivative error propagation method based on Eq. 7 (Wang et al., 2021; Xu, et al., 2018).

For comparison and validation of our AIMD simulation results, X-ray total scattering techniques are utilized to investigate the atomic structure of a material system by collecting Bragg, diffuse, and inelastic scattering. The diffraction pattern of a system is obtained from the scattering pattern. In total scattering studies, the pair distribution function (PDF) is obtained by performing a Fourier transformation on the system’s scattering pattern. The PDF illustrates the probability of finding interatomic distances between pairs of atoms in a system (Egami and Billinge, 2012). The experimental PDF will be compared to the PDF of the simulated system.

A scattering experiment measures the probability of X-rays scattered at a certain angle and with a certain energy. X-ray scattering experiments collect the static scattering function of a sample, I Q → . The angle is translated to the wavevector transfer Q → , which is related the incident and scattered wave-vector of the neutron or X-ray that hit the sample. In principle, the static scattering function is only determined by the structure of a sample and does not depend on the energy of the incident particle. The intensity of the static scattering function is related to the system’s intra-particle structure factor P Q → and the inter-particle structure factor S Q → . In liquid systems, the general equation for the static structure factor is: I ( Q  → ) = n P Q  → S Q  → (8)

Here, n is the number density of particles in a sample. The inter-particle structure factor is extracted from the experimental static scattering function data after subtracting the background. This inter-particle structure factor, S Q → , is converted to the pair distribution function G r via a sine Fourier transform (8) (Egami and Billinge, 2012): G r = 2 π ∫ 0 ∞ Q → S Q → − 1 sin Q r → d Q → (9)

Total X-ray scattering data were collected on the 11-ID-B beamline at the Advanced Photon Source at Argonne National Laboratory. (Borkiewicz et al., 2019). The sample container was a 1 mm diameter quartz tube 75% filled with our salt sample. The incident X-ray wavelength was 0.2115 Å. The X-ray scattering measurements were converted to X-ray diffraction patterns with GSAS-II (Toby and Von Dreele, 2013). The diffraction patterns were then converted into experimental PDFs with PDFgetx2 (Qiu, 2004).

The 83-atom system was used to observe the effect of the total simulation time of the production runs in the calculated properties. Figure 2 show the thermo-kinetic properties of the 83-atom system with a shorter and longer total simulation time of production run. The trajectory simulation time for the shorter simulation was approximately 25 ps, and the longer simulation was approximately 100 ps. For both simulations, the first 5 ps of the production run is neglected in the calculation as the original trajectory length was 5 ps longer. The methods for calculation of the properties are given in Section 2.1.2.

Figure 2 compare the shorter vs. longer simulation time for the 83-atom system. The properties calculated were the diffusion coefficient of each species (total, Cl, K, Mg, Na) (upper left quadrant), ionic conductivity of each species (Cl, K, Mg, Na) (upper right quadrant), and the viscosity and (4) heat capacity of the system. The results of the properties were predicted for the 83-atom system with a total simulation time of 25 ps (shorter) and 100 ps (longer). The experimental values of these properties for our salt system with the same composition are limited and are not available for comparison. In Figure 2, the calculated diffusion coefficient at each temperature for the 83-atom system with shorter and longer simulation times shows that the simulation time does not have a significant effect on the results. At temperatures below 825 K, there is very little difference between the diffusion coefficients from the shorter and longer simulations. This trend can also be seen in the comparison of the calculated ionic conductivity values for the 83-atom system shorter and longer simulation times (Figure 2). The accuracy of this calculation is questionable, considering that the ionic conductivity should theoretically increase with increasing temperature and there is little difference between values at the lowest and highest temperature. In Figure 2, the calculated viscosity values are shown for the 83-atom system with shorter and longer simulation times. This shows that there is very little difference between the values from the shorter and longer simulation times. This trend can also be seen in Figure 2 for the comparison of the calculated heat capacity values.

Figure show the thermo-dynamic properties for the 142- and 83-atom systems. The methods for calculation were discussed in section 2.1.2. The trajectory used for analysis had a simulation time of approximately 100 ps and a timestep of 1fs and 2fs for the 83-atom and 142-atom systems, respectively. In Figure, the properties calculated were (a) diffusion coefficient of each species (total, Cl, K, Mg, Na), (b) ionic conductivity of each species (Cl, K, Mg, Na), and the (c) viscosity and (d) heat capacity of the 83- and 142-atom systems simulated at 723 K–973 K, respectively. The experimental values of these properties for our salt system with the same composition are limited and are not available for comparison.

Figure 3 shows the comparison of the calculate diffusion coefficient from the 83- and 142-atoms systems. The calculated diffusion coefficient for the 142-atom system is almost consistent with theoretical predictions, with temperature increasing as the diffusion coefficient increases. Figure 3 shows the ionic conductivity of the 83- and 142-atom systems. Since the ionic conductivity is a scalar of the diffusion coefficient, it is a safe assumption that it also follows the same temperature dependency trend. Figure 3 shows the viscosity values for the 83- and 142-atom systems. Theoretically, viscosity should decrease as the temperature increases, and the viscosity values for the 142-atom system follow that trend. This is one example where it is evident that the 142-atom system is more accurate than the 83-atom system. The calculations for the 142-atom simulation shows near consistency with theoretical predictions— as temperature increases, viscosity decreases. Figure 3 shows the calculated heat capacity values for the 83- and 142-atom systems. In theory, the heat capacity of the molten salt should increase with increasing temperatures. It can be seen in this comparison how the unit cell size affects the predicted properties. The values for the 142-atom system are consistent with the theory that heat capacity increases with increasing temperature. The line for the 142-atom system is quite straight without much scatter, unlike the rest of the calculations; this is concerning but the calculations have been reviewed.

This method of calculating viscosity and heat capacity has questionable accuracy because the Stokes–Einstein relationship assumes that the system has spherical particle random motion, which is not the case for fluids with intermediate range structures, such as ionic liquids. In future calculations, a more sophisticated method of statistical mechanics, such as green-kubo analysis, will be used for viscosity and heat capacity calculation. This requires the calculation of the autocorrelation function of the system’s stress tensor.

In Figure 5, the experimental PDF of the liquid system was compared to the theoretical PDF from the simulated system at the specified temperature. The trajectory was from the 142-atom production run with a simulation time of approximately 100 ps. The theoretical PDF was obtained from MDANSE. It can be seen at both 723 K and 973 K that the positions of the first and second major peak in the PDF are at the same interatomic distance, approximately 2.4 and 3.6 Å, respectively. The literature values for the bond length between Mg–Cl, Na–Cl, and K–Cl in the solid phase are 2.8, 2.36, and 2.7 Å (Bickelhaupt, Sola and Fonesca Guerra, 2007). The slight difference between the first peak and the Na–Cl bond length is due to the system being in a liquid phase. The difference between them is the magnitude of the intensity, for which it can be assumed that the simulated system is in a fully diffusive regime. The discrepancy between the two systems is caused by the calculation of the simulated PDF not considering the periodicity of the simulation cell. Non-consideration of the periodicity also makes the theoretical PDF go to 0 instead of 1. These results show that the AIMD simulation estimates both the short- and medium-range atomic structure.

To contextualize our findings and attempt to fill the gaps in the fundamental knowledge of molten salts, the theoretical viscosity values of our system are compared to the experimental viscosity values of a similar system (Figure 6). Both systems have the same main components, but the difference is the concentration of each. Wang et al. (2021) measured MgCl₂–NaCl–KCl at the concentration of (wt%) 45.98%–15.11%–38.91%. These experimental values were reported in centi-Poise, which is equal to our units of milli-pascal seconds.

The theoretical viscosity value was calculated using the Einstein–Stokes approximation. The viscosity values of this system can be calculated using this approximation because the eutectic salt exhibits Brownian motion in the molten phase. However, the disparity observed between the theoretical and experimental viscosity values may stem from the possibility that the calculation relies on the assumption of the minimal interactions between ions, a factor critical for this approximation. This results in an underestimation of viscosity values due to the interactions of ions not being negligible. It is also possible that the size of the 142-atom system is not sufficient to capture the atomic behavior. While theoretical viscosity values do not display the same magnitudes as the experimental values, they do exhibit the same trend: viscosity decreases with increasing temperature.