First-principles molten salt phase diagrams through thermodynamic integration

Molten salts are a class of high-temperature ionic fluids that have recently attracted renewed interest due to their potential applications in modular nuclear reactors¹ and thermal energy-storage systems.² A molten alkali halide salt, such as LiF, or a mixture, such as FLiNaK, may be used as a coolant instead of highly pressurized water in nuclear reactors; these salts can also act as the medium in which fuel and fission products are dissolved.³ 

Accurate knowledge of the salt’s phase diagram is critical for the design of such reactors. Experimental results serve as the ultimate benchmark of these properties, and methods such as CALPHAD can be used in the design process using parameters that are fit to experimental inputs.⁴ However, such methods often use empirical functional forms for the relevant thermodynamic quantities needed to predict phase coexistence. A more accurate method would obtain the relevant quantities—specifically free energies—from a direct description of the interactions between the constituent atoms.

Predictions of thermodynamic properties of condensed phases from atomistic simulations can employ either Monte Carlo (MC)⁵ or molecular dynamics (MD)^6,7 approaches. Each method’s accuracy depends upon the treatment of the interatomic potential energy function. This potential energy function can be calculated from first principles using Kohn–Sham electronic density-functional theory (DFT)^8,9 in ab initio molecular dynamics (AIMD) simulations, or approximated by classical force fields, such as additive pairwise potentials. MD simulations for predicting bulk phase equilibria accurately typically require system sizes containing at least 500−1000 atoms, while AIMD simulations are typically limited by computational costs to 100–200 atoms and time scales of 10–100 ps. Consequently, MD predictions of phase equilibria have typically employed classical force fields.^10–12 However, this requires explicit parameterization of empirical force fields for new materials, and even for single component systems, it is limited in accuracy for smaller temperature and pressure ranges than a first-principles method.

Machine-learned interatomic potentials promise to bridge this gap between AIMD and classical MD by using highly flexible functional forms, such as neural-network potentials (NNPs), which can better reproduce the potential energy surface from ab initio results than simpler classical force-fields.¹³ Several families of NNPs are finding increasing use for MD simulations^14–16 and essentially serve to extrapolate the DFT level predictions from smaller AIMD simulations to larger-scale MD simulations. For molten salts, NNPs have been used to predict structure, diffusivity,^17–19 shear viscosity,²⁰ equations of state, heat capacity, thermal conductivity, and phase coexistence at individual state points.²¹ However, to our knowledge, systematic mapping of the solid–liquid phase boundary of an alkali halide, such as NaCl, using either AIMD or machine-learned potentials in the pressure–temperature space has not been performed yet.

The most common way to estimate phase equilibria in MD is to carry out direct simulations of coexistence of the two phases in a large interface calculation. Such a method has been used to estimate the melting point of NaCl with machine-learned potentials previously.²¹ At a given state point, the interface will typically move to expand the thermodynamically favorable phase at the expense of the less favorable phase.²² This requires simulating large interfaces with at least 10⁴ atoms over long time scales (typically nanoseconds) and must be repeated over several state points to pinpoint the coexistence point. However, at state points close to the true coexistence point, the velocity of the interface will typically be too low to reliably capture in simulations of tractable length. Such a method also depends upon the orientation of the interface setup between the solid and liquid phases.

A more accurate approach with better resolution involves calculating the free energy difference between the phases at various state points. One can use thermodynamic integration²³ (TI) in order to obtain these relevant free energies in molecular simulations. A reversible “pseudosupercritical” pathway that directly transforms the liquid to the solid phase can be employed, where the interatomic potential U(λ) is varied continuously as a function of an introduced path parameter λ so as to establish a reversible transformation between the solid and liquid phases path at a particular state point (P, T). The resulting free energy difference between phases is calculated by integrating ∫dλ(∂U/∂λ) along the pathway. In particular, such a method avoids the problem of interfaces between two separate phases and requires much smaller system sizes (on the order of 500–1000 atoms) than the aforementioned interface coexistence technique. TI simulations have recently been used with NNPs to assess solid–liquid coexistence for uranium and solvation free energy predictions.^24–26 However, to our knowledge, no similar study has yet been performed using NNPs to compare the effects of different DFT approximations on an alkali halide phase boundary, such as NaCl.

Here, we introduce an approach using TI to combine low-cost classical force fields and more complex NNPs trained to electronic structure data in order to, respectively, combine the computational cost and accuracy advantages of each kind of interatomic potential. Briefly, our method involves performing most of the complex transformations along the pseudosupercritical pathway using a cheap additive pairwise potential and an additional bulk transformation from the classical potential to the NNP in each phase to obtain the NNP melting point. From this initially determined melting point, we can extend the phase boundary in the (P, T) space using the Clausius–Clapeyron equation.

The remainder of this paper is organized as follows: In Sec. II, we specify the classical interatomic potential and NNP parameterization details. Afterward, in Sec. III, we detail the phase equilibrium approach, starting from the prediction of a single coexistence point using TI and then using the Clausius–Clapeyron equation to extend the phase boundary in the (P, T) space. Finally, we show results of our method in Sec. IV for the NaCl solid–liquid phase boundary and for NNPs trained to AIMD data with different choices of the exchange–correlation (XC) functional. We find that the predicted phase boundaries are highly sensitive to the choice of XC functional and that those functionals that explicitly use built-in treatment of dispersion interactions agree with experiment over a much wider range of temperatures and pressures than other functionals.

We perform classical MD simulations in LAMMPS²⁷ using the standard Fumi–Tosi (FT) parameters²⁸ for NaCl. This model is also referred to as the rigid ion model (RIM) within the molten salt literature,

In the present work, the long-range Coulomb part of the FT potential is treated using the damped shifted force model,²⁹ which allows for faster computation than Ewald and particle–particle–particle–mesh (PPPM) methods.

While this simple functional form allows for rapid computation of interatomic forces, it limits the accuracy of predicted properties to specific chemical environments and thermodynamic conditions used to parameterize the model.³⁰ 

Neural networks are universal function approximators³¹ that have been proven increasingly useful for describing the complex multidimensional potential energy surfaces of atomistic systems; in an NNP, the input to the neural network is a representation of the atomic coordinates and the output is an energy.

All approaches to NNPs require some means to map the local neighbor configuration of each atom onto the input features for the neural network. It is particularly important for these features to account for rotational, translational, and atomic permutation symmetries in order for the neural network to represent the potential energy landscape of the system with practical AIMD training datasets. Prevalent current approaches for generating these features (or descriptors) from atomic configurations include atom-centered symmetry functions,³² smooth overlap of atomic positions (SOAP),³³ neighbor density bispectrum,³⁴ Coulomb matrices,³⁵ and atomic cluster expansions (ACE)³⁶.

Here, we use the SimpleNN code¹⁵ for training and evaluating neural network potentials, which implements the atom-centered symmetry function approach. In this approach, the total energy is written as a sum of atomic energies, each expressed as a neural network of several symmetry functions evaluated on the local atomic configuration. These include radial functions that effectively measure the radial density of each atom type in a finite basis and angular functions that similarly measure the angular distribution of pairs of atom types in a finite basis surrounding each atom.³² We use the default set of radial G₂ and angular G₄ symmetry functions (70 in total) implemented in the SimpleNN package with a cutoff of 6 Å.

To train neural networks using the SimpleNN code, we use the built-in principal component preprocessing to mitigate linear dependence of symmetry functions and thereby accelerate the training.¹⁵ We also use adaptive sampling of the local atomic configurations based on Gaussian density functions, which increases the weight for infrequently encountered configurations in the loss function for training and improves the transferability of the resulting potential.³⁷ Finally, we find that a standard feed-forward neural network architecture with two hidden layers of 30 nodes each (30–30) proves sufficient, with a negligible reduction in training errors with deeper or wider networks.

An appropriate dataset of AIMD calculations is critical to reliably fit the parameters in the neural networks. For the present use case, we need to ensure that the NNP is able to accurately model both solid and molten NaCl across a wide range of temperatures and pressures. Table II lists the training AIMD simulations we use. We include solid configurations in the stable rock salt and metastable cesium chloride and zinc blende structures near the melting temperature and liquid configurations at ambient and high pressures spanning a range of temperatures above the melting point. High pressure simulations are necessary in order to ensure stability of the NNPs, essentially by sampling more of the repulsive regime of the potential energy surface.²¹ 

In order to compare the effects of different functionals, we repeat the AIMD simulations, NNP training, and all subsequent calculations for four different XC functionals, starting with the most frequently used Perdew–Burke–Ernzerhof (PBE) generalized-gradient approximation (GGA).³⁸ Since PBE often underbinds solids leading to larger lattice constants, we also use its version reparameterized for solids, PBEsol.³⁹ To analyze the impact of long-range dispersion interactions, we consider two variants of dispersion corrections to PBE, namely PBE D2⁴⁰ and PBE D3,⁴¹ which have been shown to be important for accurate structure prediction in molten salts.⁴² For the remainder of this paper, we refer to these four NNPs trained to different XC functionals as NNP-PBE, NNP-PBEsol, NNP-PBE D2, and NNP-PBE D3, respectively.

Each AIMD simulation is started from an equilibrated configuration of 64 atoms using the Fumi–Tosi potential, with this size sufficient to get atomic environments extending to the symmetry function cutoff and to require no Brillouin zone sampling in the DFT. Note that this box length is larger than the ionic screening lengths, the crystal lattice size, and the cutoff length used for the NNP and should, therefore, be sufficient for training the NNP; the classical MD simulations used for property predictions with the NNPs below use much larger boxes. The AIMD simulations are performed using the open-source JDFTx software,⁴³ using a Nosé–Hoover thermostat and barostat, a time step of 1 fs, with configurations extracted for the dataset every 10 fs. Configurations are extracted at this frequency in order for a balance between obtaining sufficient training data and tractable computation time for each individual simulation. Each AIMD simulation is run for 2 ps in the NPT ensemble, except for the high pressure sweep, which consists of four snapshots chosen along a classical MD compression simulation and then simulated in AIMD as an NVE ensemble for 0.2 ps each. With enough snapshots, the NNP learns the potential energy landscape sufficiently to handle atoms approaching one another in MD to remain stable. We use a plane wave basis with kinetic energy cutoffs of 20 and 100 hartrees for wavefunctions and charge densities, respectively, as recommended for use with the Garrity-Bennett-Rabe-Vanderbilt (GBRV) ultrasoft pseudopotential set,⁴⁴ and converge the wavefunctions to an energy threshold of 10⁻⁷ hartrees at each time step. All subsequent simulations using the NNP are performed in LAMMPS.²⁷ 

We validate the trained NNP in two ways—by comparing its forces and energies to those generated via AIMD and by comparing its radial distribution functions (RDFs) to those generated by the FT potential at three different state points. This allows us to check if the forces and energies learned by the neural network are sufficient to capture the relevant structures of each phase needed for the later phase boundary calculations.

Figure 1(a) shows the correlation and errors between the NNP and AIMD (DFT) energies and forces. The energy errors are all within 26 meV/atom, which is thermal energy k_BT at ambient temperature.

The resulting structures predicted by the NNPs are checked by running larger calculations with each potential on 4096-atom simulation cells of the solid and liquid phases at different state points displayed in Fig. 1(b). In the rock salt phase, NNP-PBE predicts a less dense phase than the other potentials (RDF peaks shifted slightly to the right), which is expected due to that functional’s well-known underbinding in solids; however, the RDFs for the other NNPs largely overlap with FT. In the liquid phase, the RDFs for all four NNPs appear to overlap with FT at both ambient and high pressures, indicating that each potential has learned the appropriate structure.

We also implement a three-fold cross-validation strategy to assess the impact of training errors of each NNP on the final phase boundary results presented later in this paper. For each exchange–correlation functional, we fit three more NNPs to 2/3 of the AIMD training data, excluding a different 1/3 each time, and restart training with different randomly initialized weights in the network. The different subsets are selected to evenly span the range of configurations in the training data. A random sampling is not used to select the subsets of training data, in order that none of the subsets lose out any relevant information in the training data needed to maintain a stable potential.

All subsequent calculations for the phase boundaries use four NNPs for each functional—three NNPs using 2/3 of the training data and one NNP using all the training data. Subsequent estimates along the phase boundary are reported as the mean of the predictions made by each of these four NNPs, and error bars are reported as the standard deviation of these predictions.

The most commonly used technique for estimating the solid–liquid coexistence in MD simulations is to directly set up an interface between a crystalline configuration and a liquid configuration and allow the system to equilibrate at a specified state point.⁴⁵ If the state point is far from the phase boundary, the thermodynamically favorable phase grows at the expense of the less favorable one, whereas close to the phase boundary, the interface does not move appreciably.

We test the direct interface coexistence method as an initial benchmark for the ambient-pressure melting temperature of the FT potential and the NNP trained to the PBE functional (hereafter referred to as NNP-PBE). For each direct interface simulation, we start with a large orthorhombic supercell of NaCl on the order of 10⁴ atoms and convert the middle region to a liquid by running a high-temperature (3000 K) NVT simulation for 30 ps, while the remainder of the system is excluded from time-integration during this initial melt. We subsequently reset velocities of the entire system and run NPT simulations to equilibrate the system at each candidate temperature at ambient pressure. We monitor the fraction of atoms in the crystalline phase as a function of time with Steinhardt’s q₆ order parameter⁴⁶ and repeat this at several candidate melting temperatures. Figure 2(a) displays a representative snapshot of the section of the system with the interface and the fraction of atoms in the solid phase as a function of time at different temperatures for the FT potential.

An inherent drawback of this method is that any assessment of a coexisting state point depends upon the orientation of the interface setup between the solid phase and the liquid phase. We repeat these simulations with solid surfaces oriented along the (100), (110), and (111) facets and track the rate of change of solid fraction as a function of temperature for each, displayed in Figs. 2(b) and 2(c) for the Fumi–Tosi and NNP-PBE potentials, respectively.

While a rough value of T_m can be estimated at 1070 ± 30 K and 875 ± 40 K for the FT and NNP-PBE potentials, respectively, we see in Fig. 2(b) that there is a qualitative difference in behavior for each facet and, hence, that such a method may have limitations for estimating coexistence of the bulk of a material with high accuracy. Moreover, the velocity of the interface is essentially a product of a mobility term and a driving force, and this driving force is proportional to |T − T_m|. Close to the melting temperature, the velocity of the interface may be too slow to observe in a tractable simulation on the order of nanoseconds.²² 

Predicting coexistence from the free energy difference between phases as a function of thermodynamic variables is more a robust method and can be achieved with significantly smaller systems than the ones necessary in direct interface coexistence simulations above.⁴⁷ 

The resulting ensemble averages $⟨∂U∂λ⟩$ depend upon the way λ is introduced to the interatomic potential energy function. Since the free energy is a state variable, it does not matter whether the pathway or the intermediate states are physically realistic and only the reversibility of the path is important; such simulations are often referred to as alchemical simulations in the community.

We can employ a “pseudosupercritical” pathway¹⁰ (PSCP) to obtain the free energy difference between the solid and liquid phase at a single state point (P, T) for a given interatomic potential. This pathway involves a direct transformation of the liquid to the solid phase at an arbitrary temperature, allowing one to sidestep the long simulation trajectories that would typically be needed to sample a nucleation event. We can also use TI to transform from relatively expensive NNPs to cheaper additive pairwise potentials, such as FT potential at the start and end of said pathway. This allows us to obtain final results that only depend on the NNP, even though most of the computation over the pathway is performed with pair potentials. Once an initial coexistence point is found, the phase boundary can be extended in the (P, T) space by the integration of the Clausius–Clapeyron equation⁴⁸ using data from simulations, for which numerous techniques are available in the literature.

So, our overall scheme for predicting the first-principles solid–liquid phase boundary of NaCl using NNPs (trained to DFT) can be broken down into the following three major steps:

1. Use a pseudosupercritical pathway at ambient pressure with the FT potential to calculate the solid–liquid free energy difference as a function of temperature ΔG_sl(T) to estimate that model’s melting temperature $TmFT$⁠.

2. Use a thermodynamic cycle involving TI from NNP to FT to obtain $TmNNP$⁠.

3. Extend the phase boundary T_m(P) by integrating the Clausius–Clapeyron equation.

We detail each of the steps in Secs. III B 1–III B 3 and note that the second and third steps are repeated for NNPs trained to different DFT functionals and with different training sets for the cross-validation and error estimation strategy discussed in Sec. II B 4.

The melting temperature of NaCl at ambient pressure has been predicted previously for the Fumi–Tosi potential using different TI approaches. The initial usage of the PSCP method, proposed by Eike et al., yielded $TmFT=(1089±8)$ K.¹⁰ More recent studies using this same pathway by Aragones et al.^49,50 but run in both directions (solid to liquid and liquid to solid) yielded an estimate of $TmFT=1082±13$ K. An alternative approach proposed by Anwar et al., involving separate pathways transforming the solid to a harmonic crystal and the liquid to an ideal gas, yielded $TmFT=(1064±14)$ K.¹¹ 

Both usages of the PSCP computed the solid–liquid free energy difference ΔG_sl at a single arbitrary temperature and used analytical corrections based on the solid and liquid equations of state to obtain a correction factor to predict the melting point.¹⁰ We adapt the PSCP method and make it more robust by running it independently at multiple temperatures to directly obtain ΔG_sl(T) and extract T_m from its zero-crossing, which reduces the possibility of systematic errors in analytic corrections.

We note that we only run this pathway in a single direction, from the liquid to the solid phase. At a single state point, this pathway consists of the following four steps, displayed in Fig. 3(a):

1. Deform liquid from its equilibrium volume to the equilibrium volume of solid at same (P, T), with free energy $ΔAdeform=−∫VLVSPdV$⁠.

2. Scale down Fumi–Tosi interaction potential U_FT to ηU_FT, with free energy given by Eq. (2) applied to U(λ) = (1 − λ)U_FT + λ(ηU_FT). Using η = 0.1, this transforms the ionic liquid to a weakly interacting liquid, amenable for the next step of transformation to the relevant crystal structure.

3. Switch on a tethering potential U_tether, consisting of attractive Gaussian potentials of the form $U(r)=−Ae−Br2$⁠, with A = 2.0 eV and B = 1.1 Å⁻² at each crystal site, which interacts with the corresponding species of atoms. This path has net potential, U(λ) = (1 − λ) (ηU_FT) + λ(ηU_FT + U_tether), with free energy given by Eq. (2) and transforms the weakly interacting liquid to a “tethered crystal.”

4. Restore the original Fumi–Tosi potential and simultaneously switch off the tethering potential. This path has net potential, U(λ) = (1 − λ) (ηU_FT + U_tether) + λU_FT, with free energy given by Eq. (2) and transforms the “tethered crystal” to the ionic crystal.

Adding the free energy from these four steps yields the net Helmholtz free energy difference between phases, ΔA_sl. We can then calculate the Gibbs free energy difference ΔG_sl = ΔA_sl + PΔV_sl, where ΔV_sl is the corresponding change in volume at this state point (P, T).

We perform each of these steps in a cell with 256 Na–Cl ion pairs, in the NVT ensemble using the Nosé–Hoover thermostat in LAMMPS. The simulations use a time step of 1 fs and converged within 25 ps for each of 50 λ values in steps 2 and 4 and within 50 ps for each λ in step 3 above. The final configuration at each λ point is used as the initial configuration for the simulation at the next λ point in order to ensure a smooth transformation along the pathway.

We repeat this entire process to compute ΔG_sl(T) for several temperatures ranging from 1030 to 1140 K and fit an ensemble of kernel ridge models fit to different resamplings of the data to extract the zero-crossing with an error estimate, displayed in Fig. 3(b). We thereby estimate the Fumi–Tosi melting point, $TmFT=1060.8±5.5$ K, which is consistent with the previous estimate from Ref. 11 but is lower than the one from Refs. 10, 49, and 50. As mentioned earlier, both these papers only ran the PSCP at a single temperature and subsequently used analytical corrections from the equations of state to estimate $TmFT$⁠. Our estimate may be closer to the method employing transformations of each phase to its reference state because we run the PSCP at multiple temperatures to directly obtain the Gibbs free energy difference between phases as a function of temperature. We use our estimate of the ambient-pressure $TmFT$ as the starting point to determine the NNP melting temperatures.

Once we have an estimate for the Fumi–Tosi melting point, we can use TI to obtain an estimate for the NNP melting point. We could adapt the approach mentioned in Sec. III B 1 by simply adding on a bulk transformation between the NNP and the Fumi–Tosi potential at the start and end of the pseudosupercritical pathway [Fig. 3(a)] to obtain $ΔGslNNP(T)$⁠; however, this would involve running multiple equilibrium simulations with the NNP at intermediate λ points at every temperature—typical NNP simulations are on the order of 100 times slower than a FT simulation, so we propose an alternative approach that allows more rapid convergence for $TmNNP$⁠.

Starting from an initial guessed value for $TmNNP$⁠, we run the following thermodynamic cycle separately for each phase:

• Convert NNP to Fumi–Tosi potential through TI, with ΔA computed using Eq. (2) applied to U(λ) = (1 − λ)U_NNP + λU_FT. At a fixed volume, this changes the pressure from P for NNP equilibrium to P′ for FT equilibrium. Consequently, this step yields ΔG = ΔA + V(P′ − P).

• Change pressure (at fixed T) of Fumi–Tosi solid/liquid from P′ back to P.

• Change temperature (at fixed P) of Fumi–Tosi solid/liquid to $TmFT$⁠, where $ΔGslFT=0$ by definition.

Figure 4(c) shows a representative dU/dλ curve obtained for the NNP to FT connection with ten λ points. Note the near perfect linearity of dU/dλ with λ, indicating that this TI can be performed with very few λ points, possibly even with just three λ points at 0, 0.5, and 1. Once again, this indicates the feasibility of the present approach to keep most of the computation at the cheaper classical potential level, requiring very few calculations using NNPs. The main requirement is that the classical potential is just accurate enough to predict a liquid and a solid phase that remain stable through the TI paths shown above.

We then converge from this straight-line approximation to find the point where ΔG_sl = 0, using a method very similar to the calculation of the NNP melting point in Sec. III B 2. We run a compression step from P → P′ at constant T and then iteratively run heating steps from T → T′ at constant P′ until convergence, using Eq. (3) at each iteration to obtain the correction factors. Note that these are essentially the last two steps of the thermodynamic cycle in Fig. 4(b).

We find that this approach converges within three iterations for each step in pressure, with a convergence criterion of |ΔT| < 4 K, for all the interaction potentials used in this study. This approach is similar to the coexistence-line free-energy difference integration method¹² but distinct from Gibbs–Duhem integration.⁵¹ 

The techniques developed in Sec. III allow for mapping of the solid–liquid phase boundary in the pressure–temperature space for any interatomic potential, including classical potentials, such as FT potential for NaCl, and machine-learned potentials, such as NNPs. Figure 5 shows the comparison of the NaCl phase boundaries predicted by NNPs trained to four different DFT XC functionals against experimental measurements and the FT potential.

First, note that while the classical potential is accurate for the melting point at ambient pressure compared to experiment, it deviates from experiment at higher pressures, consistent with previous results with this potential.¹⁰ Note that even the slope dP/dT of the coexistence line is incorrect near ambient pressure, indicating that the error stems from either the predicted enthalpy difference or molar volume difference between the phases, as indicated by the Clausius–Clapeyron equation. Table III shows that the FT potential is reasonably accurate for the enthalpy difference and solid volume but overestimates the liquid volume and thereby results in a smaller dP/dT than experiment.

The NNP-PBE potential leads to a consistently lower melting point than experiment for all pressures, as shown in Fig. 5(a). This is expected given the tendency of PBE to underestimate binding in solids generally. In particular, for NaCl, PBE predicts an almost 2% larger lattice constant for the crystal than experiment and underestimates the atomization energy by 6% (Table III), leading to the $∼15$% underestimation of the melting point.

The PBEsol functional is a reparameterization of PBE that restores the correct gradient expansion of the correlation energy, generally improving performance for solids.³⁹ This fixes the lattice constant of the crystal (⁠$<0.2$% error), but the atomization energy is still underestimated by 3%. Correspondingly, the NNP-PBEsol melting point predictions shown in Fig. 5(b) are slightly improved compared to the PBE case but are still substantially lower than experiment at all pressures.

Only NNPs trained to DFT that includes dispersion corrections predict melting points that agree reasonably with experiment, as shown in Figs. 5(c) and 5(d). This was also recently pointed out for the ambient-pressure melting point in a previous study.²¹ The dispersion-corrected PBE D2 variant⁴⁰ has the lowest error for this ambient-pressure melting point, but the PBE D3 variant⁴¹ exhibits better accuracy overall for the entire range of pressures considered in this study.

Table III indicates that both PBE D2 and PBE D3 are less accurate than PBEsol for the lattice constant at lower temperatures but are more accurate for the atomization energy. Both these dispersion-corrected functionals are more accurate for the solid and liquid volumes near the melting point, so the relative underbinding of these functionals for the perfect crystal becomes less important for solids at higher temperatures and for the liquid phase. PBE D3 has both the closest molar volumes and enthalpy difference compared to experiment among all the potentials considered here (including the FT potential), resulting in the most accurate ambient pressure dP/dT and in phase boundary across the pressure–temperature space. As shown in Table IV, this leads to the best accuracy overall across the considered pressures for PBE D3, followed by the (empirical) Fumi–Tosi model and then the other DFT functionals investigated here.

We introduced a computational approach to efficiently predict ab initio level solid–liquid phase boundaries in molten salts using a combination of machine-learned potentials and thermodynamic integration. We used NNPs trained to DFT with different exchange–correlation functionals in order to compare the accuracy of different DFT methods for the thermodynamics of molten salts, with error bars on all predictions using ensembles of NNPs trained to different AIMD data.

Most importantly, we tailored the thermodynamic integration approach to carry out most of the simulations using low-cost classical potentials, with NNPs used only in the final connection. Critically, once this approach is converged, the final result depends only on the NNP interaction potential, even though we used the lower-level classical potential in all intermediate steps of the path connecting the solid and liquid. Overall, this approach makes it much more tractable to explore molten salt equilibria with accuracy ultimately limited by the first-principles methods underlying the NNPs.

In particular, for the melting of NaCl, we show that treatment of long-range dispersion interactions in the DFT exchange–correlation functional is critical, with PBE D3 yielding the overall highest accuracy for solid–liquid coexistence across a wide range of pressures. We note that this finding is consistent with Ref. 42 that found PBE D3 to yield the highest accuracy for the local structure of molten chloride salts, among the XC functionals investigated in that study.

The overall approach described here was prototyped using NaCl as a model system but is applicable for any single-component system with an interatomic potential that is accurate enough to exhibit a stable solid and a liquid phase in the relevant temperature range. Importantly, the thermodynamic integration approach removes the dependence of the final results on this potential and allows prediction via NNPs using any underlying first-principles method. Future work can extend such free-energy methods combining NNPs and thermodynamic integration to predict phase boundaries in multicomponent systems and solubility limits.

The supplementary material includes the pseudocode of all relevant algorithms and a complete set of neural network potential files, LAMMPS, and Python scripts for reproducing the results and analysis presented here.

This work was supported by funding from the DOE Office of Nuclear Energy’s Nuclear Energy University Program (NEUP) under Award No. DE-NE0008946.

The authors have no conflicts to disclose.

T.S. and K.F. contributed equally to this work.

Tanooj Shah: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Kamron Fazel: Data curation (lead); Formal analysis (equal); Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Jie Lian: Conceptualization (equal). Liping Huang: Conceptualization (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Yunfeng Shi: Conceptualization (equal); Methodology (equal); Resources (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Ravishankar Sundararaman: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Project administration (equal); Resources (equal); Software (lead); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.