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

Abstract

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⁻¹. 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⁻¹. A quantum correction of 13.7 ± 1.9 kcal mol⁻¹ 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⁻¹. The result is in good agreement with a value of −163.5 kcal mol⁻¹ 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.

---

1 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.⁶ 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-field molecular dynamics (FFMD) simulations with classical potentials such as the Born–Mayer–Huggins–Tosi–Fumi (BMHTF) rigid-ion potential^12–15 have been demonstrated to lack full predictive capabilities due to the lack of many-body interactions that influence local structure. Polarizable ion models (PIM) developed in conjunction with quantum mechanical calculations have led to significant 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–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 artificial 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⁴ 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–30 are capable of accurately predicting molten salt structure, heat capacity, self-diffusion 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, fission/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–37 Experimental thermodynamic data (from phase equilibria measurements and activities of solution species) and/or AIMD simulations are used to empirically fit models and then predict stable phases at different temperatures or compositions. As mentioned above it would be difficult, expensive, and time-consuming 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 above-mentioned difficulties, however. The classical PIM potential force field has been used to calculate the activity coefficient ratios for lanthanide cations (e.g. U³⁺ to Y³⁺ in the LiCl/KCl eutectic mixture³⁹) and the free energy change for the reaction of the Eu³⁺/Eu²⁺ 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.⁴²

The cutting-edge NNIP-MD methods hold significant 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–45 To our knowledge, this is the first application of deep learning methods to the calculation of excess thermodynamic properties of molten salts.

2 Methods

2.1 Theoretical methods

The chemical potential of the solute X (Na⁺ or Cl⁻ ion) solvated in the molten salt liquid phase can be expressed as

μ_X = μ^id_X + μ^ex_X (1)

where μ^id_X is the ideal contribution, and the second term μ^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 defined as ε_X = U_solution − U_solvent − U_solute, where U is the system total potential energy. Then the excess chemical potential is given by

μ^ex_X = −kT ln〈e^−ε_X/kT〉₀ (2)

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–45

μ^ex_X = −kT ln〈e^−M_λ/kT〉₀+kT ln〈e^−M_λ/kT〉_εX−kT ln〈e^−ε_X/kT〉_Mλ (3)

where M_λ is a repulsive potential (half-harmonic potential in the current work) that pushes solvent ions away to the distance λ. The first term (packing, PK) is the free energy change to grow a cavity of radius λ 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 ε_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_λ indicates that the solute X is absent and there is a cavity generated by the repulsive potential M_λ during the simulation (referred to as “uncoupled” below).

The LR contribution can also be written as^43,44

μ^ex_LR = kT ln〈e^ε_X/kT〉_Mλ+εX (4)

where the subscript M_λ + ε_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)⁴⁶ by slowly growing in the repulsive potential M_λ using the functional form

f_λ(γ)=γ³M_λ, (5)

where γ is the coupling parameter that evolves from 0 to 1 (in the current work the step-size is 0.1).

The repulsive potential M_λ(r) is a half-harmonic potential

where k = 100 kcal mol⁻¹ Å⁻² and λ = 4.0 Å as exploited in previous AIMD calculations⁴⁵ 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_λ(γ) during the simulation, while the inner-shell contribution is given bywhere the solute X fully interacts with solvent ions including the potential function f_λ(γ) during the simulation.

For the long-ranged contribution, a cumulant expansion to second order yields the (uncoupled) expression

where the δε_X is the fluctuation term for the solute–solvent ions interaction energy. The corresponding coupled expression is

The regularization⁴⁷ due to the repulsive potential M_λ (that pushes the solvent molecules/ions away from the solute) produces near-Gaussian statistics for ε_X, and thus the two fluctuation terms of the uncoupled (eqn (9)) and the coupled (eqn (10)) samplings⁴⁴ 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,⁴³

the LR free energy contribution μ^ex_LR is exactly equal to the interaction energy ε at the intersection point,^46,48,49 and we utilize this equality below. Above P₀(ε) is the uncoupled distribution and P(ε) is the coupled distribution.

2.2 Computational methods

In the following section, we discuss the implementation of the AIMD simulations with CP2K 2.6.1,⁵⁰ the training process for the NNIP model with DeePMD-Kit,^51,52 and the protocol for the NNIP-MD simulations with LAMMPS.⁵³

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⁻) fixed at the center of a periodic cubic box. The box size is determined by , where the ion number density is ρ_n = 15.61 (nm)⁻³,⁴² and r_c = 4.0 γ Å. The coupling parameter γ is varied from 0 to 1 during the thermodynamic integration for the PK and IS contributions.

The initial configurations were generated during classical molecular dynamics simulations using the GROMACS 4.5.5 package,⁵⁴ and the ions were modeled with the OPLS-AA force-field.⁵⁵ With 1 fs as the time step, all classical force-field based simulations were run for 1.5 ns in the NVT ensemble after 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.⁴² This includes the dual basis sets of Gaussian-type orbitals (double zeta bases, DZVP-GTH) and plane waves with a 600 Ry cutoff.⁵⁶ Atomic cores were modeled with the Goedecker-Teter-Hutter pseudopotentials (GTH).⁵⁷ The Perdew–Burke–Ernzerhof (PBE)⁵⁸ 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.²⁸ While more advanced dispersion methods such as D3 could be used,⁶¹ 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,⁴² we find 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⁶⁴ 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.⁶⁵ 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 configurations, 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–28,42 We performed 9-step integrations for the PK and IS calculations, and we ran simulations for 10 ps (20 000 configurations) 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.⁵¹ We set up the axes with the first 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²_ij,y_ij/R²_ij,z_ij/R²_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 feed-forward DNN of 5 hidden layers with decreasing numbers of neurons, (512, 256, 128, 32, 8). Each neuron takes data Dⁱⁿ_l from the previous layers and outputs D^out_k to the next layer, implementing the linear transformation D̃_k = ∑_lw_klDⁱⁿ_l + b_k, followed by the non-linear transformation D^out_k = φ(D̃_k), where the hyperbolic tangent function φ 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 aswhere ΔE and ΔF_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_ε and p_f are the adjustable pre-factors. As the training proceeded, p_ε 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⁶⁶ was used to minimize the loss function and find 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 configurations for each of the first 9 steps of the thermodynamic integration. For the last step, where γ = 1, 50 000 configurations 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 configurations of the LR simulations, respectively. We also tried to obtain a single NNIP model by using all of the collected AIMD data. This all-in-one model produced an error of magnitude about 3 kcal mol⁻¹ 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.

2.2.3 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,⁴⁵ 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.⁴⁴ 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 first 250 ps for equilibration and the subsequent 750 ps for data production. A time step of 0.5 fs was utilized and the configurations were recorded every 0.01 ps (every 20 steps).

3 Results and discussion

We first 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 configurations 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⁻¹ deviation from AIMD for the uncoupled simulation and 1.3 kcal mol⁻¹ deviation for the coupled case. The deviations for the solute Cl⁻ are −0.4 kcal mol⁻¹ (uncoupled) and 0.4 kcal mol⁻¹ (coupled). These values are close to chemical accuracy (∼1 kcal mol⁻¹) and indicate sufficient accuracy of the NNIP model for the free energy calculations.

Fig. 1 Validation of NNIP-MD simulations in comparison with AIMD simulations. Panel (a) is the NNIP-MD interaction energy extrapolation. The energies are calculated for 1000 configurations sampled with/without a (coupled/uncoupled) solute ion (Na⁺/Cl⁻) at the center of a cavity of radius 4 Å in the molten NaCl liquid of 64 solvent ion pairs. Panel (b) and (c) are radial distribution functions (g(r)) calculated over the configurations for LR contributions, where curves are from NNIP-MD simulations. Symbols are from AIMD simulations. Panel (b) is for the uncoupled sampling (without the solute ion centered in the cavity of 4 Å). Panel (c) is for the coupled sampling (with the solute ion at the center of the cavity).

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 first maximum position at 2.77 Å and a first 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 significantly 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 configurations 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⁻¹ with a fluctuation contribution (eqn (9)) of 28.2 kcal mol⁻¹, while the mean coupled interaction energy is −180.81 kcal mol⁻¹ with a fluctuation contribution of 21.3 kcal mol⁻¹. The difference between the two fluctuation terms leads to an error of 3.5 kcal mol⁻¹ 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

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 ε and standard deviation σ. 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 γ increases from 0 to 1. The IS and PK contributions are listed in Table 1. Panel (d) is the distributions of interaction energy corrections ε_DFT − ε_MP2 for the solute ion with 7 and 17 solvent ion pairs. The sampling is over 400 cluster configurations carved from DFT simulation trajectories. The DFT calculation is under Periodic Boundary Conditions (PBC) in a cell of size 25.4 Å.

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⁻¹. In Panel (b) of Fig. 2, the distributions for the Cl⁻ ion exhibit more closely Gaussian behavior since the fluctuation terms are closer: 25.2 kcal mol⁻¹ for uncoupled simulation with a mean value of 84.05 kcal mol⁻¹ and 27.6 kcal mol⁻¹ for coupled simulation with mean value of 21.40 kcal mol⁻¹.

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⁻¹. 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.⁴⁵ 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 configurations 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⁶⁹ reveals that there is a potential shift 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⁻ ions share the same PK contribution (24.2 kcal mol⁻¹), while the IS contribution to the Na⁺ solvation free energy is −44.7 kcal mol⁻¹ and that for the Cl⁻ ion is −43.6 kcal mol⁻¹.

The numerical results for both the Na⁺ and Cl⁻ ions are listed in Table 1. The finite-size correction term is due to the artificial effect of the Ewald potential on the free energy. It is shown to be , where Q is the net charge in unit of elementary charge e, L is the box size, ξ = −2.837 297 for a cubic lattice and = 332.063301 kcal mol⁻¹ Å⁻¹ × 10².⁴³ The summation of the excess free energy for both solute ions with 256 solvent ion pairs is −178.5 kcal mol⁻¹, which is converged to within the standard error of the mean (SEM) of ±1.1 kcal mol⁻¹ 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 configurations. To the best of our knowledge, this is the first prediction of the absolute solvation free energy for ions in a molten salt using deep learning techniques.

Table 1 The excess chemical potentials (kcal mol⁻¹) 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 γ = 0.1. The calculations were performed on the OSC Pitzer cluster (Dual Intel Xeon 6148s Skylakes 2.4 GHz and 192 GB RAM). AIMD used 400 cores (4 nodes) and NNIP-MD used 40 cores (1 node)

Solute	μ ^PKX	μ ^ISX	μ ^LRX	FS	μ ^exX	AIMD	NNIP-MD
Na^+	25.4(0.2)	−44.8(0.2)	−142.8(0.7)	−28.8	−191.0(0.8)
Cl^−	25.4(0.2)	−43.8(0.2)	63.6(0.7)	−28.8	16.4(0.7)
NaCl-64					−174.6(1.1)	1.5 × 10^−3	5.5 × 10^−7
Na^+	24.5(0.3)	−44.6(0.2)	−149.2(0.7)	−23.1	−192.4(0.8)
Cl^−	24.5(0.3)	−43.7(0.2)	58.0(0.7)	−23.1	15.7(0.8)
NaCl-128					−176.7(1.1)	3.6 × 10^−3	2.9 × 10^−7
Na^+	24.2(0.2)	−44.7(0.2)	−156.0(0.7)	−18.4	−194.9(0.8)
Cl^−	24.2(0.2)	−43.4(0.2)	54.0(0.7)	−18.4	16.4(0.8)
NaCl-256					−178.5(1.1)	4.0 × 10^−3	2.7 × 10^−7
Na^+	24.3(0.2)	−44.4(0.2)	−158.8(0.7)	−14.7	−193.6(0.8)
Cl^−	24.3(0.2)	−43.6(0.2)	48.4(0.7)	−14.7	14.4(0.8)
NaCl-512					−179.1(1.1)	6.6 × 10^−3	2.5 × 10^−7

---

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.⁴² 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⁻¹. Consequently our QCT calculation at the DFT level over-estimates the excess Gibbs binding free energy by roughly −15 kcal mol⁻¹.

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 ρ_o = 15.61(nm)⁻³ and T₁ = 1150 K (NVT ensemble). The ideal gas of ions are changed into the NPT ensemble with pressure ρ_o = 1 bar and temperature of 1150 K. Then the temperature of the ideal gases is reduced to T₀ = 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₁ = 1150 K. The change of free energy from the NVT ensemble (ρ_o,T₁) to the NPT ensemble (p_o,T₁) of the liquid molten NaCl(l) is neglected.⁴³ 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⁻¹.

In a previous study by Gray-Weale et al.,⁴² it was noted that the DFT calculation over-estimates the Na atom solvation free energy in liquid Na by −19 kcal mol⁻¹ and the total Gibbs free energy change for the redox reaction Na(l)+ Cl₂(g) ⇌ NaCl(l) by −15 kcal mol⁻¹. 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.⁴⁵

Herein to estimate the correction for the free energy based on the underlying DFT simulations, we implement higher-level density-fitted MP2 theory⁷⁰ calculations over 400 cluster configurations with the basis set aug-cc-pvdz^71,72 in the Psi4 package.⁷³ The cluster configurations 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 Δμ^ex is given by

where the first term is the interaction energy correction and the second term is the fluctuation correction. The subscript ε_DFT indicates that the sampling is over configurations 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 flipped when inserting the data into eqn (13).)

Inclusion of the D2 correction over-estimates the interaction energy by −6.4 kcal mol⁻¹ on average. As mentioned above, the over-binding between ions is also indicated by the 7% over-estimation of the molten NaCl density at 1150 K.⁴² The PBC interactions (related to the Ewald potential^45,62,63) for the simulation cell of size 16.0 Å contributes about 10.0 kcal mol⁻¹ to the interaction energy correction for the solute ion pair, while in the cell of size 25.4 Å it contributes about 2.0 kcal mol⁻¹. The distributions of the interaction energy corrections under PBC in the cell of size 25.4 Å are presented in Panel (d) of Fig. 2.

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⁻¹. Assuming the dispersion interaction energy is proportional to the solvent ion density (which is approximately a constant outside the first 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⁻¹ for Na⁺ and −8.8(0.8) kcal mol⁻¹ 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⁻¹ or +20.3 kcal mol⁻¹ in the desired DFT → MP2 direction. The Ewald correction in this direction is −1.6 kcal mol⁻¹.

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 σ 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	−1.9(0.5)[2.3]	−0.7(0.4)[3.5]	−3.1(0.4)[4.7]	−1.9(0.4)[3.9]
2) DFT + D2	−8.7(0.6)[2.8]	−6.5(0.5)[3.5]	−10.4(0.5)[5.3]	−7.9(0.3)[4.8]
3) DFT + D2 PBC 16.0	−15.9(0.6)[4.2]	10.8(0.5)[4.4]	−31.3(0.6)[9.6]	23.6(0.7)[7.0]
4) DFT + D2 PBC 25.4	−10.5(0.6)[2.8]	−2.5(0.4)[3.5]	−16.7(0.5)[3.4]	0.0(0.4)[3.4]
5) D2, (2)-(1)	−6.8(0.8)	−5.8(0.6)	−7.3(0.6)	−6.0(0.5)
6) PBC 16.0,(3)-(2)	−7.2(0.8)	17.3(0.7)	−20.9(0.8)	31.5(0.8)
7) PBC 25.4,(4)-(2)	−1.8(0.8)	4.0(0.6)	−6.3(0.7)	7.9(0.5)

---

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 fluctuation 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⁻¹ 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⁻¹ for the solute ions. Compared to the experimental reference value of −163.5 kcal mol⁻¹, the calculated total excess chemical potential of −164.8(2.2) kcal mol⁻¹ provides strong initial validation of the methodology.

4 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 significantly to roughly 1 kcal mol⁻¹. Due to the over-estimation of the magnitude of the attractive interaction energy 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⁻¹ to the free energy is required to make a quantifiable 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 fission products) in the liquid phase.^40–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.

Author contributions

Yu Shi: conceptualization, methodology, software, validation, investigation, data curation, visualization, writing – original draft. Stephen T. Lam: methodology, validation, writing – review&editing. Thomas L. Beck: conceptualization, methodology, supervision, writing – review&editing, project administration, funding acquisition.

Conflicts of interest

There are no conflicts 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).

Acknowledgements

We acknowledge NSF grants CHE-1565632 and CHE-1955161 at the University of Cincinnati for financial support of this research. The computations were performed at the Ohio Supercomputer Center and the Advanced Research Computing Center in University of Cincinnati. Y. Shi acknowledges the support of the College of Arts and Sciences at the University of Cincinnati. S. Lam is supported by the Department of Energy, Office of Nuclear Energy, Nuclear Energy University Program under award no. DE-NE0009204. Research at Oak Ridge National Laboratory is supported under contract DE-AC05-00OR22725 from the U.S. Department of Energy to UT-Battelle, LLC. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE ASCR Office of Science User Facility supported under Contract DE-AC05-000R22725.

References

1. C. Le Brun, J. Nucl. Mater., 2007, 360, 1–5.
2. D. F. Williams and K. T. Clarno, Nucl. Technol., 2008, 163, 330–343.
3. S. Jerome, A. Michel, B. Ondrej, D. Sylvie, F. Olga, G. Veronique, H. Daniel, H. David, I. Victor, K. J. Leen, L. L. andMerle Lucotte Elsa, U. Jan, Y. Ritsuo and Z. Dai, Prog. Nucl. Energy, 2014, 77, 308–319.
4. H. Zhang, J. Baeyens, J. Degreve and G. Caceres, Renewable Sustainable Energy Rev., 2013, 22, 466–481.
5. U. Pelay, L. Luo, Y. Fan, D. Stitou and M. Rood, Renewable Sustainable Energy Rev., 2017, 79, 82–100.
6. M. Mehos, C. Turchi, J. Vidal, M. Wagner, Z. Ma, C. Ho, W. Kolb, C. Andraka and A. Kruizenga, Concentrating Solar Power Gen3 Demonstration Roadmap, National renewable energy laboratory technical report, NREL/TP-5500-67464, 2017.
7. F. G. Edwards, J. E. Enderby, R. A. Howe and D. I. Page, J. Phys. C: Solid State Phys., 1975, 8, 3483–3490.
8. S. Biggin and J. E. Enderby, J. Phys. C: Solid State Phys., 1984, 17, 977–985.
9. H. Ohno and K. Furukawa, J. Chem. Soc., Faraday Trans. 1, 1981, 77, 1981–1985.
10. G. J. Janz and N. P. Bansal, J. Phys. Chem. Ref. Data, 1982, 11, 505–693.
11. A. Di Cicco, M. Minicucci and A. Filipponi, Phys. Rev. Lett., 1997, 78, 460–463.
12. J. E. Mayer, J. Chem. Phys., 1933, 1, 270–279.
13. M. L. Huggins and J. E. Mayer, J. Chem. Phys., 1933, 1, 643–646.
14. F. Fumi and M. Tosi, J. Phys. Chem. Solids, 1964, 25, 31–43.
15. M. Tosi and F. Fumi, J. Phys. Chem. Solids, 1964, 25, 45–52.
16. M. Sangster and M. Dixon, Adv. Phys., 1976, 25, 247–342.
17. M. Wilson and P. A. Madden, J. Phys.: Condens. Matter, 1993, 5, 2687–2706.
18. N. Galamba and B. J. Costa Cabral, J. Chem. Phys., 2007, 126, 124502.
19. A. Bengtson, H. O. Nam, S. Saha, R. Sakidja and D. Morgan, Comput. Mater. Sci., 2014, 83, 362–370.
20. H. O. Nam and D. Morgan, J. Nucl. Mater., 2015, 465, 224–235.
21. X. Li, J. Song, S. Shi, L. Yan, Z. Zhang, T. Jiang and S. Peng, J. Phys. Chem. A, 2017, 121, 571–578.
22. H. Nam, A. Bengtson, K. Vortler, S. Saha, R. Sakidja and D. Morgan, J. Nucl. Mater., 2014, 449, 148–157.
23. H. Guo, J. Li, H. Zhang, T. Li, J. Luo, X. Yu, S. Wu and C. Zong, Chem. Phys. Lett., 2019, 730, 587–593.
24. F. Wu, S. Roy, A. S. Ivanov, S. K. Gill, M. Topsakal, E. Dooryhee, M. Abeykoon, G. Kwon, L. C. Gallington, P. Halstenberg, B. Layne, Y. Ishii, S. M. Mahurin, S. Dai, V. S. Bryantsev and C. J. Margulis, J. Phys. Chem. Lett., 2019, 10, 7603–7610.
25. S. Roy, F. Wu, H. Wang, A. S. Ivanov, S. Sharma, P. Halstenberg, S. K. Gill, A. M. Abeykoon, G. Kwon and M. Topsakal, et al. , Phys. Chem. Chem. Phys., 2020, 22, 22900–22917.
26. S. T. Lam, Q.-J. Li, J. Mailoa, C. Forsberg, R. Ballinger and J. Li, J. Mater. Chem. A, 2021, 9, 1784–1794.
27. S. Tovey, A. Narayanan Krishnamoorthy, G. Sivaraman, J. Guo, C. Benmore, A. Heuer and C. Holm, J. Phys. Chem. C, 2020, 124, 25760–25768.
28. Q. Li, E. Kucukbenli, S. Lam, B. Khaykovich, E. Kaxiras and J. Li, Cell Rep. Phys. Sci., 2021, 2, 100359.
29. S. T. Lam, Q.-J. Li, R. Ballinger, C. Forsberg and J. Li, ACS Appl. Mater. Interfaces, 2021, 13, 24582–24592.
30. A. Rodriguez, S. Lam and M. Hu, ACS Appl. Mater. Interfaces, 2021, 13, 55367–55379.
31. J. K. Startt, C. Deo and S. S. Raiman, Modelling Salt-Material Interactions in Molten Salt Reactors: Opportunities and Path., Oak ridge national laboratory technical report, 2019.
32. J. J. McMurray, T. Besmann and S. Utlak, Status of the molten salt thermodynamic database, MSTDB, Oak ridge national laboratory technical report, 2019.
33. J. W. McMurray and T. M. Besmann, in Thermodynamic Modeling of Nuclear Fuel Materials, ed. W. Andreoni and S. Yip, Springer International Publishing, 2020, pp. 2335–2363.
34. T. M. Besmann and J. Schorne-Pinto, Thermo, 2021, 1, 168–178.
35. J. Schorne-Pinto, J. A. Yingling, M. S. Christian, A. M. Mofrad, M. A. A. Aslani and T. M. Besmann, ACS Omega, 2022, 7, 362–371.
36. C. Gueneau, N. Dupin, L. Kjellqvist, E. Geiger, M. Kurata, S. Gosse, E. Corcoran, A. Quaini, R. Hania, A. Smith, M. Piro, T. Besmann, P. Turchi, J. Dumas, M. Welland, T. Ogata, B. Lee, J. Kennedy, C. Adkins, M. Bankhead and D. Costa, Calphad, 2021, 72, 102212.
37. C. Bale, E. Belisle, P. Chartrand, S. Decterov, G. Eriksson, A. Gheribi, K. Hack, I.-H. Jung, Y.-B. Kang, J. Melancon, A. Pelton, S. Petersen, C. Robelin, J. Sangster, P. Spencer and M.-A. Van Ende, Calphad, 2016, 54, 35–53.
38. V. R. Manga, S. Bringuier, J. Paul, S. Jayaraman, P. Lucas, P. Deymier and K. Muralidharan, Calphad, 2014, 46, 176–183.
39. M. Salanne, C. Simon, P. Turq and P. A. Madden, J. Phys. Chem. B, 2008, 112, 1177–1183.
40. M. A. Pounds, M. Salanne and P. A. Madden, Mol. Phys., 2015, 113, 2451–2462.
41. C. Kwon, S. H. Noh, H. Chun, I. S. Hwang and B. Han, Int. J. Energy Res., 2018, 42, 2757–2765.
42. M. Liu, P. Masset and A. Gray-Weale, J. Electrochem. Soc., 2014, 161(8), E3042–E3048.
43. T. L. Beck, M. E. Paulaitis and L. R. Pratt, The Potential Distribution Theorem and Models of Molecular Solutions., Cambridge University Press, 2006.
44. Y. Shi and T. L. Beck, J. Chem. Phys., 2013, 139, 044504.
45. Y. Shi and T. L. Beck, Proc. Natl. Acad. Sci., 2020, 117, 30151–30158.
46. L. R. Pratt and D. Asthagiri, in Free Energy Calculations-Theory and Applications in Chemistry and Biology, Springer Berlin HeidelbergNew York, 2006, pp. 323–352.
47. V. Weber, S. Merchant and D. Asthagiri, J. Chem. Phys., 2011, 135, 181101.
48. T. L. Beck, J. Stat. Phys., 2011, 145, 335–354.
49. Y. Shi and T. Beck, J. Phys. Chem. B, 2017, 121, 2189–2201.
50. T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt and F. Schiffmann, et al., J. Chem. Phys., 2020, 152, 194103.
51. H. Wang, L. Zhang, J. Han and E. Weinan, Comput. Phys. Commun., 2018, 228, 178–184.
52. L. Zhang, J. Han, H. Wang, R. Car and E. Weinan, Phys. Rev. Lett., 2018, 120, 143001.
53. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ‘t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al., Comp. Phys. Comm., 2022, 271, 10817.
54. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
55. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
56. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
57. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B, 1996, 54, 1703–1710.
58. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
59. S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473.
60. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
61. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
62. K. Leung and M. Marsman, J. Chem. Phys., 2007, 127, 154722.
63. S. W. de Leeuw, J. W. Perram, E. R. Smith and J. S. Rowlinson, Proc. R. Soc. A: Math. Phys. Eng. Sci., 1980, 373, 27–56.
64. G. Martyna, M. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
65. J. A. Morrone and R. Car, Phys. Rev. Lett., 2008, 101, 017801.
66. D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980, 2014, 1–15.
67. M. W. Chase Jr, J. Phys. Chem. Ref. Data, 1998, Monograph 9, 789–792.
68. I. Barin, in Thermochemical Data of Pure Substances, Third Edition, VCH Verlagsgesellschaft mbH, VCH Publishers, Inc., Weinheim, New York, NY, 3rd Edition, 1995, p. 1118.
69. Y. Shi, C. C. Doyle and T. L. Beck, J. Phys. Chem. Lett., 2021, 12, 10310–10317.
70. U. Bozkaya, J. Chem. Phys., 2014, 141, 124108.
71. D. E. Woon and T. H. Dunning, J. Chem. Phys., 1993, 98, 1358–1371.
72. B. P. Prascher, D. E. Woon, K. A. PeterSon, T. H. Dunning Jr and A. K. Wilson, Theor. Chem. Acc., 2011, 128, 69–82.
73. R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett, A. E. DePrince, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remigio, R. M. Richard, J. F. Gonthier, A. M. James, H. R. McAlexander, A. Kumar, M. Saitow, X. Wang, B. P. Pritchard, P. Verma, H. F. Schaefer, K. Patkowski, R. A. King, E. F. Valeev, F. A. Evangelista, J. M. Turney, T. D. Crawford and C. D. Sherrill, J. Chem. Theory Comput., 2017, 13, 3185–3197.

---