Deep neural network based quantum simulations and quasichemical theory for accurate modeling of molten salt thermodynamics

With dual goals of efficient and accurate modeling of solvation thermodynamics in molten salt liquids, we employ ab initio molecular dynamics (AIMD) simulations, deep neural network interatomic potentials (NNIP), and quasichemical theory (QCT) to calculate the excess chemical potentials for the solute ions Na+ and Cl− in the molten NaCl liquid. NNIP-based molecular dynamics simulations accelerate the calculations by 3 orders of magnitude and reduce the uncertainty to 1 kcal mol−1. Using the Density Functional Theory (DFT) level of theory, the predicted excess chemical potential for the solute ion pair is −178.5 ± 1.1 kcal mol−1. A quantum correction of 13.7 ± 1.9 kcal mol−1 is estimated via higher-level quantum chemistry calculations, leading to a final predicted ion pair excess chemical potential of −164.8 ± 2.2 kcal mol−1. The result is in good agreement with a value of −163.5 kcal mol−1 obtained from thermo-chemical tables. This study validates the application of QCT and NNIP simulations to the molten salt liquids, allowing for significant insights into the solvation thermodynamics crucial for numerous molten salt applications.


Introduction
The study of molten salt properties has seen a major resurgence in recent decades due to promising applications in clean energy technologies such as molten salt reactors 1-3 and concentrated solar power storage. 4,5 The diverse application of molten salts is primarily due to their usage in high-temperature heat-transfer media. Molten salts are excellent candidates in these systems due to their favorable physicochemical properties (e.g., thermal conductivity, heat capacity, viscosity, etc.) and relatively low cost of production. 6 When designing and optimizing salt mixture compositions for various applications, however, exploration of a high-dimensional material space is required to select for candidates with optimal properties. In addition, it is important to understand corrosion at alloy/molten-salt interfaces and property-evolution during reactor operation. This makes the experiments (such as X-ray and neutron diffraction and electrochemical measurements 7-11 ) expensive and challenging under extreme conditions. Therefore molecular dynamics (MD) simulations have been exploited to calculate molten salt properties. Classical simulations are efficient for modelling systems over timescales on the order of nanoseconds in order to make viable predictions.
Empirical force-eld molecular dynamics (FFMD) simulations with classical potentials such as the Born-Mayer-Huggins-Tosi-Fumi (BMHTF) rigid-ion potential [12][13][14][15] have been demonstrated to lack full predictive capabilities due to the lack of many-body interactions that inuence local structure. Polarizable ion models (PIM) developed in conjunction with quantum mechanical calculations have led to signicant improvements in modeling structure, thermodynamics, and dynamic properties. 16,17 Ab initio molecular dynamics (AIMD) simulations have been shown to be capable of accurately capturing local structure and solute chemistry in comparison with experiments. [18][19][20][21][22][23][24][25][26] AIMD is computationally expensive, however, and therefore it is difficult to access long time and large length scales to predict the dynamic and thermal properties. Recently, state-of-the-art arti-cial neural network-based interatomic potential molecular dynamics simulations (NNIP-MD) have been demonstrated as a promising computational tool to explore the underlying physics of the high-dimensional molten salt compositions by simulating 10 4 atoms on the timescale of nanoseconds with an accuracy at the level of density functional theory (DFT). It has been demonstrated that the NNIP-MD studies [27][28][29][30] are capable of accurately predicting molten salt structure, heat capacity, selfdiffusion coefficients, thermal conductivity, electrical conductivity, viscosity and the melting/freezing point in comparison with experimental measurements.
Understanding phase behavior remains a great challenge at the forefront of molten salt research. In addition, during reactor operation, there is continuous evolution due to processes such as transmutations, ssion/corrosion product generation, gas bubble formation, precipitation of insoluble species and other reactions. 31,32 A challenge facing the molten salt research community is to directly model these properties starting from phase diagrams, but many of the phase diagrams have never been constructed.
The CALPHAD method is a principal method for phase diagram and molten salt database development. [32][33][34][35][36][37] Experimental thermodynamic data (from phase equilibria measurements and activities of solution species) and/or AIMD simulations are used to empirically t models and then predict stable phases at different temperatures or compositions. As mentioned above it would be difficult, expensive, and timeconsuming to obtain data for high-dimensional salt systems through either experiments or AIMD simulations. 19,38 The chemical potential is central for understanding molten salt thermodynamic properties and predicting phase behavior. Studies of the chemical potential are limited due to the abovementioned difficulties, however. The classical PIM potential force eld has been used to calculate the activity coefficient ratios for lanthanide cations (e.g. U 3+ to Y 3+ in the LiCl/KCl eutectic mixture 39 ) and the free energy change for the reaction of the Eu 3+ /Eu 2+ redox couple in molten KCl. 40,41 The Widom particle insertion method has been employed with AIMD simulations to calculate the chemical potential and the solubility of the sodium atom in molten NaCl. 42 The cutting-edge NNIP-MD methods hold signicant promise to play a vital role in exploring the thermodynamic properties since they provide quantum-level accuracy with efficiency similar to classical simulations. In the present work, molten NaCl is chosen as a prototype system to demonstrate the validity of the NNIP-based quasichemical theory (QCT) calculations, which is shown to be an efficient approach for directly calculating the solute ion solvation free energy via molecular dynamics simulations. [43][44][45] To our knowledge, this is the rst application of deep learning methods to the calculation of excess thermodynamic properties of molten salts.

Theoretical methods
The chemical potential of the solute X (Na + or Cl À ion) solvated in the molten salt liquid phase can be expressed as where m id X is the ideal contribution, and the second term m ex X is the excess chemical potential due to interactions of the solute X with the solvent ions.
The solute(X)-solvent interaction energy is dened as 3 X ¼ U solution À U solvent À U solute , where U is the system total potential energy. Then the excess chemical potential is given by where the angled brackets denote the ensemble average and the subscript "0" indicates that the solute X is absent during the simulation. In the QCT, the excess free energy is partitioned into three physical parts by manipulations involving insertion and withdrawal of a hard particle that carves out a cavity in the liquid: [43][44][45] m ex X ¼ ÀkT lnhe ÀMl/kT i 0 +kT lnhe ÀMl/kT i 3 X ÀkT lnhe À3X/kT i M l (3) where M l is a repulsive potential (half-harmonic potential in the current work) that pushes solvent ions away to the distance l.
The rst term (packing, PK) is the free energy change to grow a cavity of radius l in the liquid. The second term (inner-shell, IS) is minus the free energy change to grow the same cavity in the liquid around the solute X and the subscript 3 X indicates that the solute X is present during the simulation. The third term (long-ranged, LR) is the free energy change for inserting the solute X into the cavity center and the subscript M l indicates that the solute X is absent and there is a cavity generated by the repulsive potential M l during the simulation (referred to as "uncoupled" below).
The LR contribution can also be written as 43,44 where the subscript M l + 3 X indicates the solute X is present at the cavity center fully interacting with solvent ions during the simulation (referred to as "coupled" below). The PK and IS contributions are computed via thermodynamic integration (TI) 46 by slowly growing in the repulsive potential M l using the functional form where g is the coupling parameter that evolves from 0 to 1 (in the current work the step-size is 0.1). The repulsive potential M l (r) is a half-harmonic potential where k ¼ 100 kcal mol À1ÅÀ2 and l ¼ 4.0Å as exploited in previous AIMD calculations 45 for Na + ion solvation in water. The PK contribution is then expressed as and the solute X interacts with solvent ions through the potential function f l (g) during the simulation, while the innershell contribution is given by where the solute X fully interacts with solvent ions including the potential function f l (g) during the simulation. For the long-ranged contribution, a cumulant expansion to second order yields the (uncoupled) expression where the d3 X is the uctuation term for the solute-solvent ions interaction energy. The corresponding coupled expression is The regularization 47 due to the repulsive potential M l (that pushes the solvent molecules/ions away from the solute) produces near-Gaussian statistics for 3 X , and thus the two uctuation terms of the uncoupled (eqn (9)) and the coupled (eqn (10)) samplings 44 are of close magnitude. Consequently, the average of the two mean terms gives a relatively accurate approximation for the LR contribution.
Alternatively, at a value for which the two interaction energy distribution functions are equal, 43 the LR free energy contribution m ex LR is exactly equal to the interaction energy 3 at the intersection point, 46,48,49 and we utilize this equality below. Above P 0 (3) is the uncoupled distribution and P(3) is the coupled distribution.

Computational methods
In the following section, we discuss the implementation of the AIMD simulations with CP2K 2.6.1, 50 the training process for the NNIP model with DeePMD-Kit, 51,52 and the protocol for the NNIP-MD simulations with LAMMPS. 53 2.2.1 AIMD simulation setup. Using the QuickStep module of the CP2K package, we performed all the DFT-based simulations of molten NaCl with 64 solvent ion pairs and 1 solute ion (Na + or Cl À ) xed at the center of a periodic cubic box. The box size is determined by L ¼ 64 , where the ion number density is r n ¼ 15.61 (nm) À3 , 42 and r c ¼ 4.0 gÅ. The coupling parameter g is varied from 0 to 1 during the thermodynamic integration for the PK and IS contributions. The initial congurations were generated during classical molecular dynamics simulations using the GROMACS 4.5.5 package, 54 and the ions were modeled with the OPLS-AA force-eld. 55 With 1 fs as the time step, all classical force-eld based simulations were run for 1.5 ns in the NVT ensemble aer 1 ns of equilibration in the NPT ensemble. The temperature was set at 1150 K for the classical simulations. For the AIMD simulations, parameters were chosen based on previous studies that accurately predicted local structure, standard reduction potential and sodium solubility in molten NaCl. 42 This includes the dual basis sets of Gaussian-type orbitals (double zeta bases, DZVP-GTH) and plane waves with a 600 Ry cutoff. 56 Atomic cores were modeled with the Goedecker-Teter-Hutter pseudopotentials (GTH). 57 The Perdew-Burke-Ernzerhof (PBE) 58 functionals were used for all atoms in the system, and the D2 dispersion corrections 59,60 were utilized. D2 dispersion corrections were also employed in a recent NNIP study on molten NaCl liquid. 28 While more advanced dispersion methods such as D3 could be used, 61 systematic improvement in predictions across different properties is difficult to achieve, due to the semi-empirical nature of most computationally efficient DFT-based dispersion methods. In any case, from density predictions in molten NaCl, 42 we nd that additional corrections based on higher-level theory are likely necessary to achieve chemical accuracy.
The Ewald potential 45,62,63 was applied for the electrostatic interactions under periodic boundary conditions (PBC), and a Nosé-Hoover thermostat chain 64 of length 3 was coupled to each ion to maintain a temperature of 1150 K for all the NVT ensemble simulations. Due to the requirement of many simulations along the thermodynamic integration paths to grow the nano-scale cavities, we were not able to employ the path integral formalism for incorporating quantum effects. 65 These corrections are expected to be small at such a high temperature, however.
The time step was taken as 0.5 fs. For the LR contribution, two 25 ps simulations (50 000 congurations, coupled and uncoupled) were implemented. This is found to be sufficient for calculating ensemble-average energies based on previous AIMD studies with molten salts. [26][27][28]42 We performed 9-step integrations for the PK and IS calculations, and we ran simulations for 10 ps (20 000 congurations) for each step.
2.2.2 DNN setup for potential energy surface and forces. In the DeePMD-kit framework, a local coordinate frame should be constructed to preserve translational, rotational, and permutational (same species exchange) symmetry. 51 We set up the axes with the rst axis along the direction to the nearest atom of the same ion type and the second axis along the direction to the nearest atom of the other ion type. Within this local coordinate system, the (x ij ,y ij ,z ij ) are the Cartesian components of the distance vector R ij , where atom j is a generic neighbor of atom i. The full radial and angular information around atom i is included as D ij ¼ {1/R ij ,x ij /R 2 ij ,y ij /R 2 ij ,z ij /R 2 ij } for the closest 20 atoms and only radial information is included as D ij ¼ {1/R ij } for up to 40 more atoms within the cutoff distance of R c ¼ 6.8Å. Full details of the NNIP can be found in previous studies. 51,52 The descriptors of each atom are used as inputs into a feedforward DNN of 5 hidden layers with decreasing numbers of neurons, (512,256,128,32,8). Each neuron takes data D in l from the previous layers and outputs D out k to the next layer, implementing the linear transformationD k ¼ P l w kl D in l + b k , followed by the non-linear transformation D out k ¼ 4(D k ), where the hyperbolic tangent function 4 is used as the activation function. In the last layer, only the linear transformation is applied to produce atomic energy E i . The sum of all atomic energies then yields the total energy E. The forces on each atom were computed as negative derivatives with respect to position. The loss function was taken as where DE and DF i are the root mean square (RMS) errors of the energy of the system and the force F on atom i, N is the number of atoms, and p 3 and p f are the adjustable pre-factors. As the training proceeded, p 3 began at 0.02 and ended at 8, and p f changed from 1000 to 1. The initial learning rate was 0.0001 with a decay rate of 0.95 for 5000 total decay steps. The Adam stochastic gradient descent method 66 was used to minimize the loss function and nd the parameters w kl and b k of each hidden layers. The batch-size was 5 and the training process proceeded for 1 000 000 steps. The training data consisted of energies and forces from AIMD simulations. From the AIMD simulations, we generated 20 000 congurations for each of the rst 9 steps of the thermodynamic integration. For the last step, where g ¼ 1, 50 000 congurations were collected from the two coupled simulations (IS of Na + and Cl À ) and one uncoupled simulation (PK), respectively. We trained a model for each step of the TI process (IS and PK, 10 models in total) and another 2 models for the interaction energy calculations over the uncoupled and coupled congurations of the LR simulations, respectively. We also tried to obtain a single NNIP model by using all of the collected AIMD data. This all-inone model produced an error of magnitude about 3 kcal mol À1 for the solute-solvent interaction energy compared with the AIMD calculation, and a similar error in the free energy. Due to this relatively large error, we did not utilize the all-in-one method further.

NNIP-MD simulation setup.
DeePMD-kit provides LAMMPS support through a third-party package in order to produce classical MD simulations using the NNIP to compute the atomic interactions. In this way, large time-scale simulations are accessible with quantum accuracy. 51,52 The box size is determined in the same way as shown in the AIMD simulation setup discussion. We ran NVT simulations using the LAMMPS code for systems of each ion in the molten salt liquid. Following the calculation for the Na + ion hydration free energy, 45 a cavity of radius 4.0Å was included at the periodic box center. The variations in each contribution due to the change of the cavity size are discussed in detail in our previous calculations. 44 We applied a Nosé-Hoover thermostat chain of length 3 to maintain a temperature of 1150 K. The system sizes were calculated following the above method. The NNIP-MD simulations were run for 1000 ps with the rst 250 ps for equilibration and the subsequent 750 ps for data production. A time step of 0.5 fs was utilized and the congurations were recorded every 0.01 ps (every 20 steps).

Results and discussion
We rst present results related to the NNIP model validation. For the LR contributions, the solute-solvent ion interaction energies calculated via AIMD and NNIP-MD simulations are shown in Panel (a) of Fig. 1. Averaging over 11000 congurations drawn from NNIP-MD simulations for the systems of one solute ion and 64 solvent ion pairs, the interaction energy of the solute Na + calculated via NNIP exhibits a 1.4 kcal mol À1 deviation from AIMD for the uncoupled simulation and 1.3 kcal mol À1 deviation for the coupled case. The deviations for the solute Cl À are À0.4 kcal mol À1 (uncoupled) and 0.4 kcal mol À1 (coupled). These values are close to chemical accuracy ($1 kcal mol À1 ) and indicate sufficient accuracy of the NNIP model for the free energy calculations.
Next, as shown in Panel (b) and Panel (c) of Fig. 1, the overlapping of radial distribution functions (RDF or g(r)) indicates that NNIP sufficiently reproduces the local structure predicted by the AIMD simulations and uncovers the oscillating structures at larger distances as well. Additionally, the RDF of the Na + -Cl À pair in Panel (b) exhibits a rst maximum position at 2.77Å and a rst minimum position at 4.27Å, which are close to those reported recently 27,28 using NNIP training employing a different protocol. This indicates that the presence of a cavity with a radius of 4.0Å does not signicantly affect the average liquid structure in the nearby bulk.
Just outside of the cavity, it is observed that there is a slightly higher density of the Cl À ions than the Na + cations. Integrating to a distance of 5.5Å the coordination number of Cl À is 0.7 higher than that of Na + . This leads to a dipole layer in the vicinity of a cavity of 4Å size. When the solute ion is present at the center of the cavity as shown in Panel (c), however, both solvent ions exhibit charge-symmetrical behavior, as evidenced by the overlaps in the IS RDFs.
The LR contribution is estimated from the two distributions of solute-solvent ion interaction energies, which are calculated over congurations from uncoupled and coupled simulations. As shown in Panel (a) of Fig. 2 for the solute Na + , the mean uncoupled interaction energy is À125.62 kcal mol À1 with a uctuation contribution (eqn (9)) of 28.2 kcal mol À1 , while the mean coupled interaction energy is À180.81 kcal mol À1 with a uctuation contribution of 21.3 kcal mol À1 . The difference between the two uctuation terms leads to an error of 3.5 kcal mol À1 for the LR free energy contribution, indicating that the Gaussian approximation is inappropriate and a larger cavity is required to observe Gaussian behaviour as was previously seen in the calculation of the hydration free energy of the Na + ion in water. 44,45 As discussed in the theoretical methods section above, however, the long-range contribution is equal to the interaction energy at the intersection point of the two distributions. Thus we estimate the long-ranged contribution for Na + as À156.0 kcal mol À1 . In Panel (b) of Fig. 2, the distributions for the Cl À ion exhibit more closely Gaussian behavior since the uctuation terms are closer: 25.2 kcal mol À1 for uncoupled simulation with a mean value of 84.05 kcal mol À1 and 27.6 kcal mol À1 for coupled simulation with mean value of 21.40 kcal mol À1 .
Based on the estimation of the intersection point of the two energy distributions, we estimate the LR contribution for Cl À ion as 54.0 kcal mol À1 . The dramatic increase of the LR contribution (including the sign) for the anion is attributed to the electrostatic potential at the center of the cavity. 45 Assuming that the electrostatic interaction dominates the interactions between the center solute ion and solvent ions outside the cavity of 4Å, the average energies of both Na + and Cl À ions over the uncoupled congurations gives an estimate of the electrostatic potential at the center of the cavity as À4. 55 Volt, relative to the average potential in the bulk region of molten salt liquids.
Multipole electrostatic moment analysis of the cavity-water interfacial potential 69 reveals that there is a potential shi of À3.96 Volt from the liquid water bulk phase to the cavity center, which is primarily attributed to the water molecular Bethe potential (quadrupole) contribution. Since the current work focuses on the excess free energy for the solute ion pair (in which case the interfacial potential contributions cancel exactly), the Bethe potential for liquid NaCl is not discussed further here.
The numerical results for the IS and PK contributions are shown in Panel (c) of Fig. 2 by presenting the cumulative work computed during the NNIP-MD simulations. Both Na + and Cl À Fig. 2 The process of excess chemical potential calculation for the systems of solute Na + and Cl À ions with 256 solvent ion pairs. Panel (a) and Panel (b) are for the long-ranged (LR) contributions to the solvation free energy of the solute ions Na + and Cl À , respectively. The logarithms of the distribution of interaction energies from uncoupled (right, blue) and coupled (left, black) configurations are presented with the mean value 3 and standard deviation s. The right dashed curves and the left dash-dotted curves are the Gaussian fit to both distributions. Panel (c) is for the packing (PK) and minus the Inner-Shell (IS) cumulative contributions for both solute ions. Along with the expansion of the cavity with(IS)/without(PK) the solute ions at the center, the coupling parameter g increases from 0 to 1. The IS and PK contributions are listed in Table 1 ions share the same PK contribution (24.2 kcal mol À1 ), while the IS contribution to the Na + solvation free energy is À44.7 kcal mol À1 and that for the Cl À ion is À43.6 kcal mol À1 .
The numerical results for both the Na + and Cl À ions are listed in Table 1. The nite-size correction term is due to the articial effect of the Ewald potential on the free energy. It is shown to be 1 4p3 0 where Q is the net charge in unit of elementary charge e, L is the box size, x ¼ À2.837 297 for a cubic lattice and 1 4pe 0 ¼ 332.063301 kcal mol À1ÅÀ1 Â 10 2 . 43 The summation of the excess free energy for both solute ions with 256 solvent ion pairs is À178.5 kcal mol À1 , which is converged to within the standard error of the mean (SEM) of AE1.1 kcal mol À1 as the system size increases to 512 ion pairs. The SEM in parentheses is calculated using the block-average method (10 blocks) over the production congurations. To the best of our knowledge, this is the rst prediction of the absolute solvation free energy for ions in a molten salt using deep learning techniques. The computational cost by AIMD and NNIP-MD are shown in the last columns of Table 1, where it is shown that the efficiency of the NNIP-MD simulation is about 3 orders greater than that of AIMD for all the molten salt systems and the SEMs are reduced by about 1 order of magnitude compared to the calculation reported. 42 The AIMD cost for 128, 256 and 512 systems are estimated over 100 MD-steps.
The experimental data are analyzed in terms of a Born-Haber cycle as illustrated in Fig. 3. Assuming that all the gas states are well-approximated by the ideal gas, the sum of the solvation free energies of both solute ions is calculated from thermochemical tables 67,68 to be À163.5 kcal mol À1 . Consequently our QCT calculation at the DFT level over-estimates the excess Gibbs binding free energy by roughly À15 kcal mol À1 .
In a previous study by Gray-Weale et al., 42 it was noted that the DFT calculation over-estimates the Na atom solvation free energy in liquid Na by À19 kcal mol À1 and the total Gibbs free energy change for the redox reaction Na(l)+ 1 2 Cl 2 (g) # NaCl(l) by À15 kcal mol À1 . The molten salt density is over-estimated by up to 10% when using the PBE DFT functional with the D2 dispersion correction, while the system appears to be unstable without the D2 correction. These results indicate both the importance and the subtlety of dispersion forces in the molten salt liquids. It is not surprising that the over-estimation of the density correlates with the over-estimation of the magnitude of the excess free energy. Some over-binding is also observed in ion solvation for the aqueous solution. 45 Fig. 3 The Born-Haber cycle scheme for Na + and Cl À ion solvation. The ions are assumed to be solvated from the ideal gas state (g) to liquid state (l) with number density r o ¼ 15.61(nm) À3 and T 1 ¼ 1150 K (NVT ensemble). The ideal gas of ions are changed into the NPT ensemble with pressure r o ¼ 1 bar and temperature of 1150 K. Then the temperature of the ideal gases is reduced to T 0 ¼ 298.15 K isobarically. At room temperature the sodium ion is changed to the sodium atom and chloride ion to the chlorine radical. Both atoms are then changed into elements, respectively. The solid state (s) NaCl is formed from component elements at room temperature and 1 bar pressure. Then the solid salt are heated isobarically into the liquid phase (l) at T 1 ¼ 1150 K. The change of free energy from the NVT ensemble (r o ,T 1 ) to the NPT ensemble (p o ,T 1 ) of the liquid molten NaCl(l) is neglected. 43 All of the Gibbs free energy changes are from thermodynamic tables, 67,68 and the total change of the solvation free energy is À163.5 kcal mol À1 . Table 1 The excess chemical potentials (kcal mol À1 ) for the Na + ion and the Cl À ion in molten NaCl at 1150 K. The second column is the packing contribution. The third column is the inner-shell contribution. The fourth column is the long-ranged contribution, and the fifth column is the finite-size correction (FS). The sixth column is the excess chemical potential of each solute ion and solute ion pair. The number of solvent ion pairs is given after "NaCl". In the parentheses is the standard error of the mean (SEM) over 10 blocks. The last two columns give an estimate of the computational cost (core-hour/atom/MD-step) for the PK contribution where g ¼ 0. Herein to estimate the correction for the free energy based on the underlying DFT simulations, we implement higher-level density-tted MP2 theory 70 calculations over 400 cluster congurations with the basis set aug-cc-pvdz 71,72 in the Psi4 package. 73 The cluster congurations with 7 and 17 solvent ion pairs are carved from a trajectory generated by the DFT simulation.
The correction (in the direction DFT / MP2) to the solvation free energy of each ion Dm ex is given by where the rst term is the interaction energy correction and the second term is the uctuation correction. The subscript 3 DFT indicates that the sampling is over congurations produced by the DFT simulation. As presented in Table 2, the results of DFT calculations (with the same basis set, pseudo-potential and functional as those in the above simulation) are close to the results of the MP2 calculations obtained with the Psi4 code for the sampled clusters. (Note that in Table 2 the results are presented in reference to the MP2 data; thus, the sign of the energetic correction should be ipped when inserting the data into eqn (13).) Inclusion of the D2 correction over-estimates the interaction energy by À6.4 kcal mol À1 on average. As mentioned above, the over-binding between ions is also indicated by the 7% overestimation of the molten NaCl density at 1150 K. 42 The PBC interactions (related to the Ewald potential 45 To estimate the total free energy correction, including contributions from outside the clusters, we extrapolate the cluster interaction energy correction to 256 ion pairs. Considering the cancellation of Ewald potential contributions for the solute ion pair, we assume that it is convergent to 1.6 (1.4) kcal mol À1 . Assuming the dispersion interaction energy f 1 r 6 is proportional to the solvent ion density (which is approximately a constant outside the rst solvation shell), the spherical integration leads to the expression of dispersion energy correction as a + b/N, where N is the solvent ion pair number. Using the dispersion energy corrections presented in row number 2 of Table 2, we calculate the parameters a and b as (À11.59, 20.23) for Na + and (À8.88, 16.66) for Cl À , which leads to the extrapolation for the dispersion energy correction as À11.5(1.1) kcal mol À1 for Na + and À8.8(0.8) kcal mol À1 for Cl À (in the MP2 / DFT direction). The total dispersion energy correction for the solute ion pair is then À20.3(1.4) kcal mol À1 or +20.3 kcal mol À1 in the desired DFT / MP2 direction. The Ewald correction in this direction is À1.6 kcal mol À1 . Summing up the total correction (PBC/Ewald and dispersion) of the solute ions, we see that the correction from the cluster calculations with 17 solvent ion pairs dominates 89.3% of that with 256 solvent ion pairs. As a result, the uctuation term in eqn (13) is assumed to be primarily captured by the interaction energy correction of 17 solvent ion pair cluster. This leads to a À2.5(0.1) kcal mol À1 correction for each of the Na + and Cl À solute ions. Inserting these results into eqn (13), we obtain the total correction (in the DFT / MP2 direction) as +13.7(1.9) kcal mol À1 for the solute ions. Compared to the experimental reference value of À163.5 kcal mol À1 , the calculated total excess chemical potential of À164.8(2.2) kcal mol À1 provides strong initial validation of the methodology.

Conclusions
Through investigation of the solvation thermodynamics of the Na + and Cl À ions in the molten NaCl liquid, we have presented and validated an efficient and general methodology to calculate the solvation free energy of ionic species in the molten salts. The methodology incorporates ab initio molecular dynamics simulations, interatomic potentials based on deep neural network models, and quasichemical theory. Efficient molecular dynamics simulations with the NNIP model reduces the calculation uncertainty signicantly to roughly 1 kcal mol À1 . Due to the overestimation of the magnitude of the attractive interaction energy Table 2 Corrections for the interaction energy of the solute ion (Na + or Cl À ) with solvent ion pairs (7 and 17). The correction is referenced to the MP2 calculation using the Psi4 quantum chemistry package. In parenthesis is the SEM. In the square bracket is the standard deviation s of the sampling. The configurations are carved out from a CP2K simulation trajectory without a cavity around the solute ion. The DFT calculations are implemented via the CP2K package. The MP2 and the first two DFT calculations are for the isolated systems. The subsequent two DFT calculations are under Periodic Boundary Condition (PBC) in the simulation cell of size 16.0Å and 25.4Å, respectively. The row number 5 shows the dispersion correction of D2 contributions. The PBC effects on the corrections are presented in the last two rows Na + 07 Cl À 07 Na + 17 Cl À 17 0) MP2 0.0 0.0 0.0 0.0 1) DFT of the solute ion with the solvent ions at the DFT level (with D2 dispersion corrections), a high-level quantum chemical correction of roughly 15 kcal mol À1 to the free energy is required to make a quantiable prediction of the excess free energy. These results highlight the importance of dispersion and polarization interactions in the molten salt liquids.
The methodology sets the stage for larger-scale simulations of molten salt mixtures at a quantum mechanical level of accuracy, allowing for quantitative investigations of the activities, solubilities, and redox potentials of ionic species (including the corrosion and ssion products) in the liquid phase. [40][41][42] The methodology holds the potential to provide essential data for molten salt applications in both concentrated solar energy storage materials and molten salt reactors.

Data availability
Data related to the calculations presented in this paper will be made freely available upon request to the authors.

Conflicts of interest
There are no conicts to declare.

Notes
This manuscript has been authored in part by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The publisher acknowledges the US government license to provide public access under the DOE Public Access Plan (https://energy.gov/downloads/doe-public-access-plan).