The impact of hydrogen valence on its bonding and transport in molten fluoride salts

Abstract

Interest in molten salts has increased significantly over the last decade due to their potential application in various clean-energy technologies including hydrogen generation, solar heat storage, advanced fission nuclear power plants, and compact fusion energy systems. In nuclear fission and fusion power plants, high heat capacity molten salts allow operation at high temperature and atmospheric pressure, which could dramatically increase efficiency, reduce capital cost, and enable passive safety features. In many of these systems, the hydrogen isotope is of particular importance due to its ability to corrode structural materials as ³H⁺ in fluoride salts, and its potential to cause significant radioactive release as diffusive ³H⁰, which are cited as key barriers to technological deployment. Yet, the chemistry and transport behavior of the hydrogen species remain poorly understood due to the difficulties in handling toxic salts and radioactive materials. Here, using ab initio molecular dynamics, we present a coupled examination of hydrogen speciation in the most common prototypical salts 66.6% LiF–33.3% BeF₂ (Flibe) and 46.5% LiF–11.5% NaF–42% KF (Flinak). Using extensively validated calculations on the local structure and dynamics, we find significant difference between ³H⁰ and ³H⁺ transport behaviors that are usually overlooked. We find that ³H⁰, which always exists as H₂, diffuses 3–5 times faster than ³H⁺, which can be ascribed to hydrogen bonding and complexation in solution. This work explains contradicting experimental results and provides useful species transport data for designing hydrogen capture and corrosion control systems for molten salts.

---

1 Introduction

Interest in high-temperature molten salts has greatly increased in the last decade due to their application in a variety of emerging clean energy technologies including thermal storage, hydrogen production, batteries, fuel cells, and nuclear power.^1–5 In many of these salt systems, corrosion remains a central challenge.⁶ This is particularly true for halide salts where corrosion impurities are soluble in the salt, so a passivation layer between the structural materials and working fluid cannot be established. In halide salts, atmospheric impurities such as hydrogen, oxygen, moisture and hydroxides are thermodynamically unstable and are consequently chemically reduced by the corrosion of metallic components of structural materials.^7,8 The presence of impurities has been shown to cause high rates of initial corrosion, resulting in significant mass loss and material degradation.⁹

Fluoride molten salts are of particular interest in the development of next-generation nuclear reactors, since they offer a combination of high heat capacity, high boiling point, and low neutron absorption.¹⁰ While existing water-cooled nuclear power generates more than 50% of the carbon-free electricity in the U.S. today, total nuclear capacity is projected to steadily decrease towards 2050 and new nuclear installations will be limited by concerns of safety, and economic viability.^11,12 A step-change improvement in nuclear technology is therefore required for nuclear power to significantly contribute to decarbonization of the electrical grid. Thus, there has been a significant push towards the development of advanced molten salt nuclear systems, which could enable greatly increased heat removal, thermodynamic efficiency, and operation at near-atmospheric pressure. As a result, the use of molten salts could enable design simplification, economic competitiveness, and inherent safety via passive heat rejection.

Molten-salt cooled advanced reactors of all types have proposed the use of lithium-based salts since the addition of lithium depresses the melting point, making them easier to handle.^13–15 Among candidate salts, 66.6LiF–33.3BeF₄ (Flibe) and 46.5LiF–11.5KF–42NaF (Flinak) are the most popular choices for their superior heat transfer properties.^16,17 One of the greatest known challenges however, is that radioactive tritium (³H) is generated by the neutronic reaction ⁶LiF + n → 4He + ³HF.¹⁸ The presence of ³HF drives up the redox potential in the system,⁷ as hydrogen is less electropositive than lithium, which causes corrosion of metallic species such as chromium via the equilibrium reaction: 2³HF_(g) + Cr_(s) ↔ CrF_2(soln) + ³H_2(g), forming diatomic tritium gas. The relative concentration of the Mⁿ⁺/M (M = Cr, Fe, Ni or other structural metal elements) and HF/H₂ couples depend on the temperature and redox potential, which is established by the Nernst equation. The diffusion of hydrogen/tritium into structural materials can cause hydrogen-induced embrittlement, as well as tritium permeation and tritium inventory problems that are nuclear safety concerns. Furthermore, corrosion products like CrF₂ have the potential to precipitate in the system, which leads to heat exchanger fouling and the presence of these dissolved species can also change the physical, thermodynamic and chemical properties of a molten salt mixture.⁸ In order to minimize oxidation in structural steels, the salt must be sufficiently reducing such that the vast majority of the tritium is converted to H₂, with a HF to H₂ ratio less than 0.01. On the other hand, the potential must remain sufficiently high (oxidizing) to prevent precious metal precipitation, resulting in a narrow operating window [U_lower, U_upper] (U_upper − U_lower) ∼ 0.1 V for Flibe) at reactor operating temperatures. Charge neutral H₂ in coolant fluid and atomic H in some metals are highly diffusive and must be continuously removed from the system and captured to minimize the potential for a significant radiological release.¹⁹

The removal of tritium from salt is a coupled problem of redox chemistry and transport that requires understanding and accurate prediction of the in situ behavior of possible tritium species. Differences in diffusion and reaction rates can result in differential species removal rates. Depending on the system design, this can change the relative quantities of the H oxidation states present in the salt, in which ½ H₂(molten-salt) ↔ H⁺(molten-salt) + e⁻(U), where the left-hand side is one-half of a charge-neutral molecule that is solvated in the liquid molten salt, and the right-hand side is a proton/triton cation living in the same liquid phase with dynamic coordinations of fluorides, plus an excess quantum mechanical free electron living in the metal electrode of Fermi energy E_F = −eU. At thermodynamic equilibrium, the concentration quotient [H⁺]/[H₂]^½ in the molten-salt phase can shift the voltage U of the metal.^20,21 By metal–salt interactions (corrosion or deposition), U can also affect the concentrations of other charge-neutral or charged species (impurities or fission products) in the melt, which can produce unwanted changes in performance. Thus, the valence-dependent transport and chemical properties of tritium must be precisely quantified and understood. Yet, little reliable experimental data is available due to the costly handling of toxic salts and radioactive tritium. Thus, tritium control has been repeatedly emphasized as a bottleneck in the road to developing any reactor using lithium-based molten salts.^22,23

In this work, we examine the chemical behavior of tritium species in both prototypical salts Flibe and Flinak. First-principles simulation is used to interpret experimental results and predict salt properties, elucidating the chemical structure and transport properties required to predict tritium evolution in the system. This work establishes a link between tritium chemistry and transport behaviors over the thermodynamic conditions for a typical operating reactor. By comparing disparate experimental results with predictions from first principles, we highlight experimental methods and provide further guidance for measuring radioactive salts. Lastly, our validated calculations also generate useful engineering data that can be directly used in the design of tritium or salt handling systems.

2 Tritium structure and transport in molten salt

Since few experiments were performed and the redox chemistry in the system is usually not controlled, identifying the particular species of interest (HF or H₂) remains difficult. For tritium transport properties, the reported values spread over orders of magnitude with great uncertainty.^24–29 The experimental uncertainty depends on the technique (capillary, membrane, pulsed-gradient nuclear magnetic resonance), but can generally be expected to be in the ±20% range.^30–33 For structural analysis, X-ray and neutron diffractions are used. In practice, these methods should be coupled to simulation and other experimental methods to fully resolve features of multi-component systems.^34–36

In this work, tritium local structure, chemistry and transport are investigated using ab initio simulation with realistic dynamics and thermochemistry.^37,38 Dilute tritium as ³H⁺ and ³H⁰ are separately introduced into the prototypical molten salts. The simulation and analysis protocols are described in detail in the methods section. In addition to tritium-in-salt systems, calculated properties of various fluoride salts including LiF, KF, LiF–KF, BeF₂, NaF, Flibe and Flinak are compared and supported by experimental data, increasing confidence in simulation data. These salts were chosen due to their chemical similarity with the prototypical salts of interest, availability of data, and relevance in fuel processing, spent fuel recycling, cooling advanced fission reactors, or as potential for tritium fuel production in fusion devices.^39–42 The results of these supporting calculations are discussed in ESI.†

3 Results

3.1 Structure and chemistry of tritium in fluoride salts

3.1.1 Radial coordination of fluorine in Fibe and Flinak. To examine the local structures of both possible tritium oxidation states (of ³H⁺ and ³H⁰), Fig. 1 shows the local coordination of fluorine and the H–F radial distribution functions in Flibe and Flinak. In Flibe, tritium as ³H⁺ shows a distinct peak with a maximum at the radial distance of 0.9 Å, which is close to the radii of covalent bond at 0.88 Å.⁴³ In comparison, the H–F peak for ³H⁰ is wide and shallow with first peak maximum located at 2.7 Å, indicating limited coordination between H⁰ and F. This suggests limited ionic dissolution of the ³H⁰ in Flibe. The coordination number N(r) of fluorine to tritium is calculated by the integral of the RDF, which is also shown in Fig. 1. This further confirms the peak shown by the RDF, where a plateau can be seen for ³H⁺ but not for ³H⁰. Over the trajectory, the coordination number of fluorine to tritium at each time step was tallied based on the cutoff distance of the first minimum of the RDF. It was found that the ³H⁺ was coordinated to two fluorine atoms for 61% of the simulation and one fluorine atom for 39% of the simulation.

Fig. 1 Radial distribution functions and running coordination number of ³H–F at prototypical operating temperature of 973 K in (a) Flibe and (b) Flinak. The running coordination number is calculated by integration of the RDF.

For Flinak, H–F coordination is very similar to that of Flibe. The H–F RDF, where tritium has a + 1 oxidation state shows a distinct peak at 1.1 Å, with a calculated coordination number of 1.83, which indicates a coordination of either 1 or 2 atoms. Coordination tallying shows that the tritium was coordinated to two fluorine atoms 95% of the time, and only one fluorine for 5% of the simulation. The coordination number and peak distances are summarized in Table 1.

Table 1 First peak radius and coordination number of ³H⁰/³H⁺ and F⁻ pairs in Flibe and Flinak at 973 K (design temperature for typical reactor). The coordination number is determined by integration of the first shell of the RDF. Numbers are excluded for ³H⁰ due to lack of clear solvation shell

Salt	Ion pair	First peak radius (Å)	Coordination number
Flibe	H^+–F^−	0.9	1.76
	H^0–F^−	2.7	—
Flinak	H^+–F^−	1.1	1.85
	H^0–F^−	2.7	—

---

3.1.2. Tritium coordination complexes and reactions in Flibe. While coordination counting and RDF analysis gives a general indication of the nearest neighbor structure ³H⁺ and ³H⁰, it does not provide more detailed information about chemical transformation and speciation. In order to understand the chemical reaction and time-dependent coordination in extended systems, a graph theoretical approach combined with a hidden Markov model is used to automatically detect the generation and destruction of new complexes along the trajectory.^44,45

For ³H⁰ in Flibe, no reactions were found and tritium existed solely in H₂ with a H–H interatomic distance of 0.76 Å at 973 K. The dissolution of ³H⁺ is confirmed in Fig. 2, which shows the interatomic distance between a single ³H³H⁺ atom and different fluorine atoms in the mixture. ³H⁺ becomes coordinated to fluorine with an average interatomic distance of 1 Å, and jumps between different fluorine atoms during the simulation. In addition, tritium can coordinate to two fluorine atoms simultaneously as shown between t = 15–24 ps (the overlapping time period on 1 and 2 in Fig. 2a) with fluorine atoms ‘1’ and ‘2’ (see Fig. 2b for the snapshots). In contrast, ³H⁰ remains un-coordinated to fluorine and the interatomic distance is always greater than 1.5 Å. Thus, ³H⁺ appears much more reactive in Flibe than ³H³H⁰. This suggests that ³H⁺ has a higher solubility than ³H⁰ in Flibe, which agrees with the relative magnitudes of their Henry's coefficients in literature.^46,47

Fig. 2 (a) ³H³H⁺ in Flibe, interatomic distance between tritium and fluorine atoms in the solution. The hopping H⁺ occurs in sequence 1–5 shown below. Step 1 to 2 clearly show the existence of the ³HF₂⁻ complex where the ³H⁺ is bound to two fluorine atoms. (b) Visualization of tritium hopping between two fluorine atoms labelled ‘1’ and 2’. The simulation temperature is 973 K.

We found ³H⁺ exists in various chemical forms, including HF, HF₂⁻, BeF₄HF²⁻, and Be₂F₇HF³⁻, in the fractions shown in Fig. 3d. A snapshot of the simulation is shown in Fig. 3a, numerical representation of molecules is shown in Fig. 3b, and an example reaction is shown in Fig. 3c. Simple fluorides HF and HF₂⁻ were the most common states with mean life-time τ = 4.4 and 1.8 ps at 973 K respectively, while tritium bound to BeF₄²⁻ (τ = 0.43 ps) tetrahedral and tetrahedral polymer chains were found to exist less than 15% of the simulation time (60 ps). The “other” category in Fig. 3d consists of rare (single-event) structures making up the remainder of the simulation, which includes free tritium ion or tritium complexes with larger BeF complexes such as Be₃F₁₀HF⁴⁻ (τ = 0.33 ps) or Be₄F₁₃HF⁻⁵ (τ = 0.097 ps), which existed for less than 0.1% of the simulation. With increasing temperature, free tritium and other transition states become more common, occupying 5% of the simulation at 973 K, and up to 20% of at 1373 K. The average lifetime of H⁺ complexes reduces with increasing temperature with τ = 1.4, 1.0 and 0.8 ps for 973, 1173 and 1373 K respectively.

Fig. 3 (a) Snapshot of HF in Flibe, (b) equivalent representations of Be₂F₇ corner-sharing tetrahedral as an adjacency matrix and graph of nodes and vertices used to find chemical structures in system, (c) relationship between chemical structures found across spatial and temporal dimensions to construct reaction coordinates, (d) fraction of each molecule at 973 K, 1173 K, and 1373 K.

Reactions with H⁺ involving the most common complexes are shown in Table 2, all of which involve the combination with HF to an existing atom or molecule in the solution. The first reaction shows that the formation of the HF₂⁻ complex from F⁻ ion and HF combination with ΔE = −0.4 eV calculated for the isolated molecules. For this reaction, the rate of association and dissociation of F⁻ with HF increases with increasing temperature, likely dependent on the amount of free fluorine that is available from the surrounding mixture. The second reaction involves a HF molecule bonding to a BeF₄²⁻ tetrahedral complex with ΔE = −0.6 eV. Reaction two is more than twice as frequent at 973 K than 1173 K and 1373 K, which corresponds to the availability of BeF₄HF²⁻ at lower temperatures as shown in Fig. 3d. The last reaction shows the bonding of an HF molecule to two corner-sharing BeF₄²⁻ tetrahedral complexes with ΔE = −0.6 eV. The ‘other’ reactions all involve less common species, such as transition complexes swapping F⁻ ions with solution, or HF breaking and joining with larger Beryllium complexes: Be₃F₁₀HF⁴⁻ ↔ BeF₃⁻ + Be₂F₇HF³⁻. The frequency of these reactions increases with increasing temperature due to the increase in transient high-energy structures that are discovered at higher temperatures.

Table 2 Tritium (³H⁺) reactions in Flibe found during simulations at 973, 1173 and 1373 K

Reaction		Count
		T = 973 K	T = 1173 K	T = 1373 K
1	F + HF ↔ HF2^−	22	32	35
2	BeF4^2− + HF ↔ BeF4HF^2−	54	16	20
3	Be2F7^− + HF ↔ Be2F7HF^−	8	18	3
Other		50	59	166

---

This analysis suggests high ³H⁺ reactivity to fluorine and beryllium mixture in Flibe, which are in agreement with the RDFs analysis and further confirm that solvation of ³H³H⁺ is significantly stronger than that of ³H⁰.

3.1.3. Tritium coordination complexes and reactions in Flinak. Examination of tritium coordination in Flinak reveals that ³H⁺ is coordinated to two fluorine atoms for most of the simulation. At 973 K, the tritium remains bonded with the same fluorine atoms (|d⃑_H–F| = 1.18 Å) for the simulation interacting with ions as HF₂⁻. At the higher temperatures, the HF₂⁻ is able to exchange F⁻ with the solution and the tritium becomes more dissociative. The most common states for ³H⁺ are shown in Fig. 4, of which the most common structure is HF²⁻ (τ = 12 ps at 973 K) which existed for 80–90% of the 60 ps simulation, followed by NaF₄HF³⁻ (τ = 0.14 ps), and NaF₅HF⁴⁻ (τ = 0.25 ps), which were found to be present for 1.5–5% of the simulation. The ‘other’ category of the simulation again includes rare atomic structures such as HF (τ = 0.16 ps), NaF₃HF²⁻ (τ = 0.09 ps) appearing at 1373 K.

Fig. 4 (a) Tritium molecules in Flinak, (b) fraction of each molecule at 973 K, 1173 K, and 1373 K.

In Flinak, all the cations are monovalent and form fewer complex molecules, in contrast to Flibe where longer chains of Be–F complexes are found. Li or K do not form detectable reaction products, indicating free lithium and potassium ions in solution. Most sodium atoms in the simulation were also freely dissociated; only ∼11% of the Na found in a polyhedral configuration (3% as NaF₄, 8% as NaF₅). Further, Na complexes usually have lifetimes on the order of 10² fs and readily react or exchange with fluorine in the solution. As a result, relatively few H⁺ reactions were found in the solution. A handful of reactions found included HF₂⁻ combination reactions with NaF complexes such as: HF₂⁻ + NaF₄³⁻ ↔ NaF₅HF⁴⁻ (ΔE = −0.5 eV), and NaF₅HF⁴⁻ + F⁻ ↔ HF₂⁻ + NaF₅⁴⁻ (ΔE = −0.1 eV). At the higher temperatures of 1173 K and 1373 K, some dissociative behavior was observed for HF₂⁻ with the fluorine ejection reaction HF₂⁻ ↔ HF + F⁻ (ΔE = 0.5 eV) observed 4 times.

In Flinak, the behavior of ³H⁰ is similar to that in Flibe, where tritium exists purely as H₂ with |d⃑_H–H| = 0.77 Å at 973 K. This again agrees with the result of the radial distribution function in Fig. 1, which shows weak solvation of hydrogen with fluorine. Thus, in both clean Flibe and Flinak, ³H⁰ shows little chemical interaction with surrounding ions compared to ³H⁺.

3.2 Tritium diffusivity in FliBe and FLiNaK

In this section, the calculation of diffusivity and activation energy of hydrogen is discussed for clean Flibe and Flinak and compared to existing data. It should be noted that the form of tritium measured in experiments is not explicitly known because metallic impurities and corroding materials can form halides in the salt and change the chemical form of tritium being measured as described in Section 1. In this study, the isotopic effect of hydrogen is expected to be small, especially relative to the differences between experiments.^48,49

3.2.1 Tritium diffusivity in Flibe. The diffusivity and activation energy of hydrogen isotopes has been measured by four studies over the temperature range of 1020–1325 K.^24–27Fig. 5 compares the results of the different experiments against the simulation of tritium as ³H⁺ and ³H⁰.

Fig. 5 Experiments versus simulation. Calderoni used tritium ³H gas diffusing across a 1-D nickel membrane, Oishi used ³H created by salt irradiation and measured diffusion in a capillary, Anderl used ²H diffusion in a cylindrical cell, and Nakamura used ¹H in a cylindrical cell. The shaded region represents 95% confidence interval of regression.^24–27

The calculated diffusivity for ³H⁰, as a result of H₂ molecular diffusion, was found to be on average 3–4 times larger than that of ³H⁺. Such a large difference in diffusivity is related to the coordination of ³H⁺ with beryllium fluoride complexes and fluorine in the mixture. The calculated diffusivity and activation energies are shown in Table 3. The diffusion equation for two oxidation states are expressed:

Table 3 Calculated diffusivities of tritium as ³H⁺ and ³H⁰ in Flibe

Temperature (K)	^3H^+ diffusivity × 10^9 (m^2 s^−1)	^3H^0 diffusivity × 10^9 (m^2 s^−1)
973	3.0 ± 0.5	6.3 ± 1.4
1173	7.0 ± 1.3	16.3 ± 3.3
1373	11.1 ± 1.6	27.0 ± 4.1

---

The calculated activation energies for ³H⁰ and ³H⁺ are similar, 40.0 ± 3.0 kJ mol⁻¹ and 36.0 ± 2.0 kJ mol⁻¹ respectively. Compared to the experimental data, the simulation results are within error bounds of Calderoni's and Oishi's experiments. In Oishi's experiment,²⁵ tritium is generated as ³H⁺ from salt irradiation. However, some of the tritium is likely to chemically react with the Ni container to generate ³H₂. In Calderoni's experiment,²⁴ hydrogen purge gas is used, which likely isotopically exchanges with ³HF to form ³H¹H in the solution. Further, both experiments purified the salt by hydro-fluorination and thus may contain some proportion of ¹H molecules which could react or exchange with ³HF or ³H₂. Thus, in both experiments, the diffusivity measurement likely corresponds to a combination of tritium species and oxidation states. The activation energy calculated using Oishi's data was 36.0 ± 6.0 kJ mol⁻¹, and the activation energy from Calderoni was found to be 42.0 ± 7.0 kJ mol⁻¹, which are similar to the computational values.

In the other two studies, data collection was limited (2–3 experimental measurements), and data significantly diverged from simulation and previous experiments. Both used diffusion in a similar cylindrical diffusion cell. In Anderl's experiment,²⁶ diffusivity measured was 0.8 × 10⁻⁵ cm² s⁻¹, which is 80% lower than computed ³H⁰ diffusivity and 40% lower than ³H⁺ calculated from eqn (1) and (2). In Anderl's experiment, impurity analysis found 600 ppm of oxygen impurity in the salt before the start of the experiment. This combined with any potential oxygen ingress could cause reactions with HF via: BeO + 2HF → BeF₂ + H₂O. In this case, the measured diffusivity could actually be that of H₂O or other intermediate reaction complexes. Further, their analytical model used diverged from their experimental data over time. This suggests a systematic error in the transport model, possibly due to increasing air ingress or other interfacial effects.²⁶ Similarly, Nakamura's experiment under-predicted by the same order of magnitude at 823 K.²⁷ The experimental activation energy was also lower than the calculated values at 25.2 ± 8 kJ mol⁻¹. In Nakamura's study, the lower diffusivities and activation energies could be caused by potential impurities (not measured) or temperature gradients in the experimental cell.

While previous studies have found that impurities can significantly increase in corrosion and change the thermochemical properties of salt, the specific impact of different impurities on tritium is not well understood.^8,9 Here, additional simulations were performed for tritium and Flibe in the presence of a single unit of Cr (dominant corrosion product) and H₂O (atmospheric moisture) at 973 K. In these cases, tritium diffusivity decreased by 25% and 15% respectively. From structure analysis, tritium was found to coordinate with the impurities (Cr–³H and ³H–O bonding), which likely resulted in the observed differences in diffusivity relative to pure Flibe. The structures discovered can be found in the ESI in Fig. S7.† Similarly, it was found that the introduction of corrosion impurities (CrF, HF) can alter the diffusivities and activation energies of Li and F ions in Flibe shown in Fig. S6.† These results support the idea that impurity effects can impact experimental results summarized in Fig. 5 and provides further evidence that impurities can significantly impact chemical and transport properties. However, many more atmospheric and fission-generated impurities can exist simultaneously, which can have unpredictable effects in real systems. As such, the role of impurities in a variety of chemical forms, oxidation states, and conditions warrants a more detailed investigation in future studies.

3.2.2 Tritium diffusivity in Flinak. The diffusivity of hydrogen isotopes in Flinak has been measured by only three experimental studies, which all used ¹H as shown in Fig. 6.^27–29 Similar to the relative behavior of tritium in Flibe, ³H⁰ diffuses faster than ³H⁺. The average diffusivity of ³H⁰ is 4× larger than that of ³H⁺ across the temperature range from 973 K to 1373 K (Table 4). The diffusion coefficients are expressed:

Fig. 6 Comparison of hydrogen diffusivity between ab initio simulations and three experiments. Experiment by Zeng and Fukada diffuse H₂ across a 1-D nickel membrane, while experiment by Nakamura was done by radial diffusion across a cylindrical Monel window. The highlighted region represents the 95% confidence interval from regression of the collected data.^27–29

Table 4 AIMD diffusivities of tritium as ³H⁺ and ³H⁰ in Flinak

Temperature (K)	^3H^+ diffusivity × 10^9 (m^2 s^−1)	^3H^0 diffusivity × 10^9 (m^2 s^−1)
973	2.9 ± 0.5	10.6 ± 2.3
1173	6.7 ± 0.7	23.4 ± 3.1
1373	8.25 ± 1.25	31.0 ± 5.1

---

The activation energies of both oxidation states are similar at 29.0 ± 6.4 kJ mol⁻¹ for ³H⁰ and 29.0 ± 9 kJ mol⁻¹ for ³H⁺ in Flinak. From analysis of chemical structure, the lower diffusivities of ³H⁺ can be attributed to the fact that the molecule exists as a larger and more chemically reactive molecule HF₂⁻ in the melt, whereas ³H⁰ exits as H₂. The lower activation energies of tritium in Flinak compared to Flibe indicate a lower diffusion barrier. This can be rationalized in terms of chemical structures in both salts; the formation of BeF₄²⁻ chains in Flibe obstructs transport while relatively free ions dominate in the Flinak solution. Like Flibe, the redox condition in Flinak can be either reducing or oxidizing depending on the impurities in the salt.

Fukada suggested that their results measured H⁰, which showed the most similarity to the calculated diffusivity of ³H⁰ compared to other studies. In this experiment, tests were performed to ensure low corrosion of nickel due to contact with Flinak. Further, both pH measurement and mass spectroscopy were performed during the diffusion measurement to confirm H₂ measurement was up to a sensitivity of 1 ppm and HF up to sensitivity of 1 ppb. Here, the simulations again confirm the importance of accurately confirming the redox and impurity state. The activation energy calculated from their data from 800 K to 1000 K was 43.0 ± 4.7 kJ mol⁻¹, which is higher than the simulated values. This is likely due to the higher temperatures used in simulation. Fukada found that, the activation energy for diffusion increased at lower temperatures due to increased mass transfer resistance at interfaces caused by increased H⁺ interaction. Thus, these results show good agreement with our simulations in general.

In contrast, corrosion impurities and chemical species were not measured explicitly in the other two studies. In both studies, this is combined with a limited number of experimental data points, yielding high uncertainty in the overall result. The activation energy calculated by Nakamura was found to be near-zero to slightly negative −5.4 ± 5.7 kJ mol⁻¹, which has not been observed in any other species transport in molten salts.³⁰ Here, we note that Nakamara's Flinak study uses the same experimental procedure as their Flibe study discussed in Section 3.2.1, both of which show significant deviations from simulations in this work and other experiments in the literature. This suggests a systematic error in the study, possibly caused by the measurement protocol (cylindrical diffusion with analytical model), or the use of hydrogen interacting materials (Monel). This effect, along with impurity interaction as examined for Flibe in Section 3.2.1, should be explored more thoroughly since the diffusivity and activation energies differ dramatically from other experimental data and these simulations.

3.3 Implications of coupled hydrogen diffusion and chemistry in molten salt

In this work, tritium chemistry and transport are examined for the two possible tritium oxidation states H⁰ and H⁺ in prototypical salts Flibe and Flinak. Accurate prediction of these properties is critical for the design of salt systems that can contain hydrogen from moisture ingress or generated by irradiation. In many systems, the tritium transport in salt could limit overall transport due to relatively slow diffusion through salt, particularly as H⁺ ion. Further, quantifying large differences in chemistry and transport behavior of different chemical forms (D⁰_H:D_H⁺ ≅ 2–5) is important in understanding and controlling the redox condition, and consequently, corrosiveness of the salt. Conventionally, the differences between the diffusion of H⁰ and H⁺ in molten salt has been ignored due to a lack of information. In some cases, it has even been assumed that D_H+ > D⁰_H, which could cause significant error in predicting transport behavior.^20,50 From this work, it is clear that experiments that controlled for speciation, impurities, and corrosion produced more reliable data than those that did not, which was supported by simulation results. Further, the simulations provide insights that allowed understanding of transport behavior at the level of local atomic structure, which would otherwise be difficult to infer from experimental data alone.

4 Conclusions

Through extensive first-principles molecular dynamics simulations, we discovered that hydrogen transport in molten fluorides depends on its oxidation state and the salt's redox condition. We examined the bonding and connectivity of ³H⁺ and ³H⁰ and established the link between chemistry and transport properties in prototype salts Flibe and Flinak. In addition, we accurately predicted the local structure (±0.1 Å peak distance, and ± 0.4 coordination) and ionic diffusivities (within experimental uncertainty, ±20%) of various binary, ternary and quaternary molten salt fluorides. Chemical analysis was performed using a combination of local radial structure, and an automated search of molecular graphs to monitor temporal changes in extended chemical structures.

Chemical-structural analyses suggest that, in fluoride salt, ³H⁺ shows soluble and dissociative behaviors in both Flibe and Flinak. Specifically, ³H⁺ coordinates with BeF₄²⁻ tetrahedral corner-sharing polymer-like structures in Flibe, while, in Flinak, ³H⁺ exists as HF₂⁻ with some association to NaF polyhedra. In contrast, ³H⁰ remained bound as a gaseous dimer H₂ in both Flibe and Flinak. As a result, the diffusivity of ³H⁰ is 3–4× and up to 5× higher than that of ³H⁺ in Flibe and Flinak, respectively. Further, the diffusion activation energies in Flibe (40.0 ± 3.0 kJ mol⁻¹ for ³H⁰ and 36.0 ± 2.0 kJ mol^{−1 3}H⁺) are higher than those in Flinak (29.0 ± 6.4 kJ mol⁻¹ for ³H⁰ and 29.0 ± 9 kJ mol⁻¹ for ³H⁺) due to the formation of BeF networks that lower overall ionic transport in Flibe, in comparison to Flinak where ions are predominately free. In both salts, the ³H⁺ activation energy is similar to ³H⁰, showing fairly weak association with similar characteristic to hydrogen bonding between hydrogen and fluorine atoms.

In interpreting the limited experimental results with wide discrepancies spanning more than an order of magnitude, our ab initio simulations provide a deeper understanding of experimental sources of error like impurity concentration, surface effects, and corrosion. Furthermore, we provide atomistic level insight into the interplay among oxidation states, structures and transport properties of tritium in both Flibe and Flinak, as well as useful engineering data that is difficult to obtain experimentally. The relative differences between H⁰ and H⁺ can greatly impact the design of impurity removal and corrosion control systems in both nuclear and non-nuclear salt-based systems. In advanced nuclear reactors, diffusive radiological release and corrosion are two of the most prominent challenges that are tightly coupled through tritium transport and chemistry. This work elucidates hydrogen behavior that must be accounted for in future experiments and modeling of molten salts.

5 Methods

5.1 Materials and salt systems

A range of fluoride systems were simulated and compared against experimental structure and transport data over a range of temperatures from 775–1373 K. These systems include LiF, KF, NaF, BeF₂, LiF–KF, LiF–BeF₂ and LiF–NaF–KF, which can be found in the ESI.† Tritium as ³H⁰ and ³H⁺ were added as a single H₂ and HF as they would be found in prototypical salts Flibe (66.6% LiF–33.3% BeF2) and Flinak (46.5% Li–11.6% NaF–41.9% KF). The Flibe consisted of 10 Be, 40 F, and 20 Li and Flinak consisted of 20 Li, 5 Na, 18 K, and 43 F to make 46.5% Li–11.6% NaF–41.9% KF.

In each system, the cell size was chosen to match the experimental density of the salt. For systems where accurate density is not known, it can be estimated by running a series of simulations to fit an equation of state, or by running the isobaric-isothermal ensemble (NPT) until equilibrium is reached. The relatively small system sizes are chosen for computational efficiency. While periodic effects might be expected for smaller system sizes, it has been shown by simulating larger cells of up to 200 atoms that properties like diffusion and local structure are generally converged on the smaller system sizes.³⁷

5.2 First-principles molecular dynamics simulations

In all of the salts, the system was initialized by randomizing the positions using Packmol software.⁵¹ The system was equilibrated under the canonical ensemble (NVT) for at least 10 picoseconds at least 200 K above the melting point. The production runs lasted for at least 60 picoseconds sampled at 1 fs time step. The temperature was controlled with a Nosé–Hoover thermostat using a Nosé mass corresponding to a period of 40 fs. Born–Oppenheimer ab initio molecular dynamics (AIMD) was performed with the Vienna Ab initio Simulation Package.^52,53 Periodic boundary conditions were used in all directions. The AIMD calculations were performed with density functional theory (DFT) using a generalized gradient approximation (GGA) functional for electron exchange and correlation energy. The projector augmented wave (PAW) Perdew–Burke–Ernzerhof (PBE) functionals used were Li_sv, Be, F, Na_pv, K_sv and H. An energy cutoff of 600 eV was used for plane wave expansion, and the Brillouin zone was sampled at the Γ-point only to improve computational speed. The parameters chosen have been shown to yield good convergence within 2 meV per atom.⁵⁴

5.3 Atomic coordination and complexation analysis

Chemical analysis was performed by two methods. The first is simple and involves the examination of coordination based on the radial distribution function (RDF). The coordination of a central atom is found by taking the sum of the number of neighbor atoms within a distance cutoff. The cutoff is specified using the first minimum of the RDF of each species pair. The coordination numbers of each element with respect to the others are tallied to produce a time-averaged distribution. The averaged coordination number over the trajectory is determined by the integration of the RDF g(r):⁵⁵where, N_A–B is the nearest neighbor coordination number, ρ_B is the number density of species B, and r_m is the radius of the first minimum of g(r). By performing the search on coordination distribution, immediate local structures are easily found.

The second method allows for the detection of extended molecular structures and reactions by using graph-theory based searching algorithms followed by statistical filtering for noise removal. The method follows that of Martinez,^45,56 which has previously been used to study reactions and transport in organic systems and surfaces.^44,57–59 At each step in the trajectory, the system is represented by an atomic graph with nodes N and edges E. The nodes represent individual atoms, while the edges represent bonds between atoms, defining the graph G_t as below:

G_t = (N:{A₁, …, A_n}, E:|R_ij ≤ R^C_ij) (6)

where R_ij is the Euclidean distance between atom i and j, A₁⋯A_n contain information on the atom's elemental species and unique index in the simulation. R^C_ij is the radial cutoff distance to register atoms as being bonded, which in this case is defined based on the atomic species of i and j. In previous work, the bond length was estimated by the sum of covalent radii.⁴³ However, this can be significantly different than the actual length of associated atoms in an ionic liquid. Thus, the cutoffs are estimated from the bond lengths found in the simulation RDFs: R^C_ij = 1.4R̄_ij where R̄_ij is the average bond length determined from the radial distribution function, and a factor of 1.4 is applied to encapsulate thermal noise due to vibrations at high temperatures. At each time step, the overall graph G_t can be separated into its component connected subgraphs g_m, which represent individual molecules. The subgraphs are found using the breadth-first search algorithm.⁶⁰

For each molecular subgraph found in the simulation, a binary time series E^m(t) is constructed valued at 1 or 0 is constructed to indicate whether the subgraph exists (value of 1) at each time step in the simulation. The observed signal E, is modeled using a two-state Hidden Markov Model (HMM) in which there exists an underlying lower-frequency signal X, representing the true signal. The HMM is defined by the joint probability:⁶¹

where the output probability matrix containing P(E_i|X_i) and the transition probability matrix containing P(X_i|X_i−1) are parameterized:

The transition probability of switching states, t_p was set between 10⁻⁶ to 10⁻⁷, which is markedly low to provide a stronger noise filter that removes false transitions in the observed signal. This is necessary due to the high-temperature liquid state causing high-frequency fluctuations in the system. The Viterbi algorithm is used to find the likely sequence X that produced observation E.

With the filtered sequences X^m for each molecule, chemical reactions are then found. The generation of a new product is characterized by a 0 → 1 transition in the sequence occurring at time τ_p and the destruction of a reactant is characterized by a 1 → 0 transition at τ_r. At each product time, the nearest reactants and products in time are located and a chemical balance is performed, checking for stoichiometric conservation among the reactants and products. If a reaction is not found, the search window for reactants and products is iteratively expanded, until chemical balance can be achieved or a maximum search window Δτ = 150 fs is reached. If no atom-conserving reaction can be found for a generated product, it is disregarded as unfiltered noise. The reaction energies are determined by performing geometry optimized DFT calculations on the reactants and products extracted from the trajectory.

5.4 Calculation of diffusivity

The diffusion coefficient is taken as self-diffusivity at the dilute limit. This applies to the single hydrogen atom in the molten salt simulation cells, and to experiments where tritium concentrations exist at the ppm level.^46,47 The diffusivity is calculated from the slope of the mean squared displacement as a function of time by the Einstein relationship:where the slope is determined for and averaged over subsets n_t of the total trajectory, and over all the atoms of the same elemental type:where r_i(t) is the coordinates of atomic i at time t, and dt is the length of time spanned by each subset in time. The total trajectory spanned more than 60 000 time steps (60 ps), of which more than 15 equal-sized sub-trajectories were sampled. Sufficient spacing was used between time origins to eliminate correlation between diffusion evaluations. The standard error and confidence intervals were determined using the block-averaging method.⁶²

Author contributions

S. L, C. F., R. B, and J. L conceived and guided the overall study. S. L and Q. J performed data generation and analysis. S. L and J. M. developed the graph-theoretic code for chemical analysis of extended structures. All authors contributed to the preparation of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

SL acknowledges funding support from the National Science and Engineering Research Council of Canada, and the Shanghai Institute of Applied Physics. QJL and JL acknowledge support by the Department of Energy, Office of Nuclear Energy, Nuclear Energy University Program (NEUP) under Award Number DE-NE0008751.

6 References

1. Y. Wang, et al., Current understanding and challenges of solar-driven hydrogen generation using polymeric photocatalysts, Nat. Energy, 2019, 4(9), 746–760.
2. S. Frangini and A. Masi, Molten carbonates for advanced and sustainable energy applications: Part II. Review of recent literature, Int. J. Hydrogen Energy, 2016, 41(42), 18971–18994.
3. T. Voisin, A. Erriguible and C. Aymonier, A new solvent system: Hydrothermal molten salt, Sci. Adv., 2020, 6(17), eaaz7770.
4. S. Bell, T. Steinberg and G. Will, Corrosion mechanisms in molten salt thermal energy storage for concentrating solar power, Renewable Sustainable Energy Rev., 2019, 114, 109328.
5. S. Delpech, C. Cabet, C. Slim and G. S. Picard, Molten fluorides for nuclear applications, Mater. Today, 2010, 13(12), 34–41.
6. K. Sridharan and T. R. Allen, Corrosion in Molten Salts. Elsevier Inc., 2013.
7. J. Zhang, et al., Redox potential control in molten salt systems for corrosion mitigation, Corros. Sci., 2018, 144, 44–53.
8. J. Zhang, Impurities in Primary Coolant Salt of FHRs: Chemistry, Impact, and Removal Methods, Energy Technology, 2019, 7(10), 1–13.
9. S. S. Raiman and S. Lee, Aggregation and data analysis of corrosion studies in molten chloride and fluoride salts, J. Nucl. Mater., 2018, 511, 523–535.
10. S. Delpech, C. Cabet, C. Slim and G. S. Picard, Molten fluorides for nuclear applications, Mater. Today, 2010, 13(12), 34–41.
11. E. Information and E. Information, Electricity generation from natural gas and renewables increases as a result of lower natural gas prices and declining costs of solar and wind renewable capacity, making these fuels increasingly competitive, 2020.
12. D. Gielen, F. Boshell, D. Saygin, M. D. Bazilian, N. Wagner and R. Gorini, The role of renewable energy in the global energy transformation, Energy Strategy Reviews, 2019, 24, 38–50.
13. A. L. Qualls, Advanced Reactor Technology Program Molten Salt Reactor Campaign FY 2018 Summary, Oak Ridge National Laboratory, Oak Ridge, Tennessee, ORNL/SPR-2018/987, 2018.
14. C. Forsberg and P. F. Peterson, Basis for Fluoride Salt-Cooled High-Temperature Reactors with Nuclear Air-Brayton Combined Cycles and Firebrick Resistance-Heated Energy Storage, Nucl. Technol., 2016, 196(1), 13–33.
15. B. N. Sorbom, et al., ARC: A compact, high-field, fusion nuclear science facility and demonstration power plant with demountable magnets, Fusion Eng. Des., 2015, 100, 378–405.
16. D. F. Williams, L. M. Toth, and K. T. T. Clarno, Assessment of Candidate Molten Salt Coolants for the Advanced High-Temperature Reactor (AHTR), Oak Ridge National Laboratory, Oak Ridge, Tennessee, ORNL/TM-2006/12, 2006.
17. O. Beneš and R. J. M. Konings, Molten salt reactor fuel and coolant. Elsevier Inc., 1st edn, 2012, vol. 3.
18. C. W. C. W. Forsberg, et al., Tritium Control and Capture in Salt-Cooled Fission and Fusion Reactors: Status, Challenges, and Path Forward, Nucl. Technol., 2017, 197(2), 119–139.
19. S. Guo, J. Zhang, W. Wu and W. Zhou, Corrosion in the molten fluoride and chloride salts and materials development for nuclear applications, Prog. Mater. Sci., 2018, 97, 448–487.
20. H. Qin, C. Wang, S. Qiu, D. Zhang, W. Tian and G. H. Su, Study of tritium transport characteristics in a transportable fluoride-salt-cooled high-temperature reactor, Int. J. Energy Res., 2018, 42(4), 1536–1550.
21. S. T. Lam, et al., Tritium Management and Control Using Carbon in a Fluoride-Salt-Cooled High-Temperature Reactor, Fusion Sci. Technol., 2017, 71(4), 644–648.
22. D. F. Williams and P. F. Britt, Technology and Applied R&D Needs For Molten Salt Chemistry, in Molten Salt Chemistry Workshop at Oak RidgeNational Laboratory, 2017.
23. D. E. Holcomb, G. F. Flanagan, G. T. Mays, W. D. Pointer, K. R. Robb, and G. L. Yoder, Fluoride Salt-Cooled High-Temperature Reactor Technology Development and Demonstration Roadmap, Oak Ridge National Laboratory, Oak Ridge, Tennessee, ORNL/TM-2013/401, 2013.
24. P. Calderoni, P. Sharpe, M. Hara and Y. Oya, Measurement of tritium permeation in flibe (2LiF-BeF2), Fusion Eng. Des., 2008, 83(7–9), 1331–1334.
25. J. Oishi, H. Moriyama, S. Maeda and Y. Asaoka, Tritium recovery from molten LiF-BeF2 salt, Fusion Eng. Des., 1989, 8, 317–321.
26. R. A. Anderl, et al., Deuterium/tritium behavior in Flibe and Flibe-facing materials, J. Nucl. Mater., 2004, 329–333, 1327–1331.
27. A. Nakamura, S. Fukada and R. Nishiumi, Hydrogen isotopes permeation in a fluoride molten salt for nuclear fusion blanket, J. Plasma Fusion Res., 2015, 11, 25–29.
28. S. Fukada and A. Morisaki, Hydrogen permeability through a mixed molten salt of LiF, NaF and KF (Flinak) as a heat-transfer fluid, J. Nucl. Mater., 2006, 358(2–3), 235–242.
29. Y. Zeng, et al., Development of the Tritium Transport Analysis Code for the Thorium-Based Molten Salt Reactor, Nucl. Technol., 2018, 203, 48–57.
30. G. J. Janz and N. P. Bansal, Molten Salts Data: Diffusion Coefficients in Single and Multi Component Salt Systems, J. Phys. Chem. Ref. Data, 1982, 11(3), 505–693.
31. I. Farnan and J. F. Stebbins, High-Temperature29Si NMR Investigation of Solid and Molten Silicates, J. Am. Chem. Soc., 1990, 112(1), 32–39.
32. A. L. Rollet, V. Sarou-Kanian and C. Bessada, Measuring self-diffusion coefficients up to 1500 K: A powerful tool to investigate the dynamics and the local structure of inorganic melts, Inorg. Chem., 2009, 48(23), 10972–10975.
33. C. S. Johnson, Diffusion ordered nuclear magnetic resonance spectroscopy: Principles and applications, Prog. Nucl. Magn. Reson. Spectrosc., 1999, 34(3–4), 203–256.
34. G. W. Neilson, A. K. Adya and S. Ansell, 8 Neutron and X-ray diffraction studies on complex liquids, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2002, 98(1), 271–320.
35. L. B. Skinner, et al., Structure of Molten CaSiO3: Neutron Diffraction Isotope Substitution with Aerodynamic Levitation and Molecular Dynamics Study, J. Phys. Chem. B, Nov. 2012, 116(45), 13439–13447.
36. H. E. Fischer, A. C. Barnes and P. S. Salmon, Neutron and x-ray diffraction studies of liquids and glasses, Rep. Prog. Phys., 2006, 69(1), 233–299.
37. S. Filipe Sousa, P. Alexandrino Fernandes and M. João Ramos, General Performance of Density Functionals, J. Phys. Chem. A, Aug. 2007, 111(42), 10439–10452.
38. A. Jain, Y. Shin and K. A. Persson, Computational predictions of energy materials using density functional theory, Nat. Rev. Mater., 2016, 1(1), 1–13.
39. O. Beneš, Thermodynamics of molten salts for nuclear applications, Institute of Chemical Technology, Prague, PhD Thesis, 2008.
40. M. Salanne, C. Simon, P. Turq, N. Ohtori, and P. A. A. Madden, Modeling of Molten Salts, in Molten Salts Chemistry, Elsevier Inc., 1st edn, 2013, no. i, pp. 1–16.
41. C. Bessada and A. L. Rollet, In Situ Spectroscopy in Molten Fluoride Salts. Elsevier Inc., 2013.
42. Y. Iwadate, Raman Spectroscopy and Pulsed Neutron Diffraction of Molten Salt Mixtures Containing Rare-Earth Trichlorides: Trial Approaches from Fundamentals to Pyrochemical Reprocessing, Elsevier Inc., 1st edn, 2013.
43. B. Cordero, et al., Covalent radii revisited, J. Chem. Soc., Dalton Trans., 2008, 21, 2832–2838.
44. M. Döntgen, M. D. Przybylski-Freund, L. C. Kröger, W. A. Kopp, A. E. Ismail and K. Leonhard, Automated Discovery of Reaction Pathways, Rate Constants, and Transition States Using Reactive Molecular Dynamics Simulations, J. Chem. Theory Comput., 2015, 11(6), 2517–2524.
45. L. P. Wang, A. Titov, R. McGibbon, F. Liu, V. S. Pande and T. J. Martínez, Discovering chemistry with an ab initio nanoreactor, Nat. Chem., 2014, 6(12), 1044–1048.
46. P. E. Field and J. H. Shaffer, The solubilities of hydrogen fluoride and deuterium fluoride in molten fluorides, J. Phys. Chem., 1967, 71(10), 3218–3222.
47. A. P. Malinauskas and D. M. Richardson, The Solubilities of Hydrogen, Deuterium, and Helium in Molten Li2BeF4, Ind. Eng. Chem. Fundam., 1974, 13(3), 242–245.
48. R. Mills and K. R. Harris, The effect of isotopic substitution on diffusion in liquids, Chem. Soc. Rev., 1976, 5, 215–231.
49. R. Mills, Search for isotope effects in heavy-ion diffusion in liquids, Nature, 1975, 256(5519), 638–639.
50. J. D. Rader, M. S. Greenwood Paul and W. W. Humrickhouse, Verification of Modelica-Based Models with Analytical Solutions for Tritium Diffusion, Nucl. Technol., 2018, 203(1), 58–65.
51. L. Martinez, R. Andrade, E. Birgin and J. Martinez, Packmol: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2008, 30, 2158–2164.
52. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys.c, 1996, 54(16), 11169–11186.
53. G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6(1), 15–50.
54. H. O. Nam and D. Morgan, Redox condition in molten salts and solute behavior: A first-principles molecular dynamics study, J. Nucl. Mater., 2015, 465, 224–235.
55. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, New York, NY, USA, Clarendon Press, 1989.
56. L. P. Wang, R. T. McGibbon, V. S. Pande and T. J. Martinez, Automated Discovery and Refinement of Reactive Molecular Dynamics Pathways, J. Chem. Theory Comput., 2016, 12(2), 638–649.
57. L. Agosta, E. G. Brandt and A. P. Lyubartsev, Diffusion and reaction pathways of water near fully hydrated TiO2surfaces from ab initio molecular dynamics, J. Chem. Phys., 2017, 147(2), 024704.
58. J. P. Mailoa, et al., A fast neural network approach for direct covariant forces prediction in complex multi-element extended systems, Nature Machine Intelligence, 2019, 1, 471–479.
59. C. W. Park, M. Kornbluth, J. Vandermause, C. Wolverton, B. Kozinsky and J. P. Mailoa, Accurate and scalable multi-element graph neural network force field and molecular dynamics with direct force architecture, 2020, pp. 1–33.
60. C. Cormen, T. H. Leiserson, C. E. Rivest and R. L. Stein, Introduction to Algorithms, MIT Press, Cambridge, MA, 2001.
61. S. J. Russell and P. Norvig, Artificial Intelligence: A Modern Approach, Pearson Education, 2003.
62. D. C. Rapaport and D. C. R. Rapaport, The Art of Molecular Dynamics Simulation. Cambridge University Press, 2004.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0ta10576g

---