X-ray scattering reveals ion clustering of dilute chromium species in molten chloride medium

Abstract

Enhancing the solar energy storage and power delivery afforded by emerging molten salt-based technologies requires a fundamental understanding of the complex interplay between structure and dynamics of the ions in the high-temperature media. Here we report results from a comprehensive study integrating synchrotron X-ray scattering experiments, ab initio molecular dynamics simulations and rate theory concepts to investigate the behavior of dilute Cr³⁺ metal ions in a molten KCl–MgCl₂ salt. Our analysis of experimental results assisted by a hybrid transition state-Marcus theory model reveals unexpected clustering of chromium species leading to the formation of persistent octahedral Cr–Cr dimers in the high-temperature low Cr³⁺ concentration melt. Furthermore, our integrated approach shows that dynamical processes in the molten salt system are primarily governed by the charge density of the constituent ions, with Cr³⁺ exhibiting the slowest short-time dynamics. These findings challenge several assumptions regarding specific ionic interactions and transport in molten salts, where aggregation of dilute species is not statistically expected, particularly at high temperature.

---

Introduction

Simple inorganic molten salts have emerged as important media for the development of modern concentrating solar power (CSP) plants and next-generation molten salt nuclear reactors (MSRs) having significant potential to provide sustainable, carbon-free energy for the future.^1–4 Among different salt formulations, mixtures of potassium chloride (KCl) and magnesium chloride (MgCl₂) hold special promise for these technologies, owing to their stability under high-temperature regimes and their exceptional thermal characteristics.^5–9 In addition, chloride salts are abundant and therefore relatively inexpensive. However, they can be very corrosive to metal alloy materials under the extreme, real-world operating conditions (873–1073 K) of CSP and MSR systems.^10–12 A major issue is related to chromium being primarily depleted from the alloys upon contact with the molten KCl–MgCl₂, leading to the formation of dilute CrCl₃ in the melt.¹⁰

Despite extensive investigations of molten salts in the past, the only insights regarding how low-concentration Cr³⁺ metal ions behave in the high-temperature ionic medium have been achieved through optical spectroscopy.^13,14 The absence of long-range order, such as that present in a single crystal or powder sample, limits the usefulness of conventional in-lab X-ray diffraction equipment using low-photon-energy sources, while the low metal content limits the information available from Raman spectroscopy (ESI Fig. 1†). X-ray absorption spectroscopy (XAS) is a powerful tool to investigate metal ion speciation,^15–18 but provides reliable information only about the first coordination sphere, and recent XAS experiments were unable to determine the local structure of Cr³⁺ in a KCl–MgCl₂ molten salt mixture at the high temperatures (>973 K)¹⁹ that are relevant to technological processes in CSP and MSR applications. The accurate identification and understanding of coordination environments and dynamic structural behaviour of chromium species could provide useful insights into the processes of corrosion and transport in molten salts. At a more fundamental level, this knowledge will broaden our general understanding of metal ion solvation in unconventional ionic media at high temperature.

In this study, we demonstrate a highly integrated approach combining theory, simulations and high-energy synchrotron X-ray scattering experiments to elucidate the atomic-scale structure and dynamics of dilute Cr³⁺ species in the KCl–MgCl₂ molten salt at 1073 K. In order to enhance the information about the local structure of Cr³⁺ in a matrix of strong X-ray scattering ions, we applied a differential pair distribution function (dPDF) method. This technique used in conjunction with ab initio molecular dynamics (AIMD) simulations, reverse Monte Carlo (RMC) modeling and theory enabled us to interpret and reveal information about local Cr³⁺ environments. While the dPDF technique has been used to investigate uranyl complexes²⁰ and CO₂ (ref. 21) in aqueous solutions as well as solvation shells around various nanoparticles,^22–24 to the best of our knowledge, this is the first instance of its use for metal ions dissolved in other highly ionic media under extreme temperature conditions. In addition, to understand the local dynamics of Cr³⁺ as compared to that of other ions in such complex medium, we have formulated and used a hybrid transition state (TS)-Marcus theory approach. Using this rate theory model, we found that dynamical processes in the melt are highly sensitive to the charge density of the constituent ions, with Cr³⁺ showing the slowest short-time dynamics. Perhaps surprisingly for such a high-temperature low-Cr³⁺ concentration system, our studies found that octahedrally coordinated Cr³⁺ ions tend to dimerize, forming persistent structures.

Results and discussion

Synchrotron X-ray diffraction measurements were initially performed on the salt mixture of 50 : 50 mol% KCl–MgCl₂. To allow for the observation of structural changes with the increase of temperature, the data were collected at 298 K, 673 K (below the melting point), and 1073 K (molten state) using a customized furnace (ESI Fig. 2†) specifically designed to hold samples in quartz capillaries at high temperature.²⁵ Our measured coherent scattering intensities provide the X-ray structure factor, S(Q), where the momentum transfer, Q, is 4π sin(θ)/λ (2θ is the scattering angle and λ is the incident X-ray wavelength). The Fourier transform of S(Q) gives pair distribution functions [G(r) or D(r) = 4πρ₀rG(r), where ρ₀ is the ionic number density] that are related to the probability of finding ion pairs with a given separation distance r^26,27 (see ESI† for S(Q), G(r) and D(r) definitions). Fig. 1a shows the evolution of D(r) as a function of temperature. Consistent with the increasing disorder and the loss of crystalline symmetries in the solid salt above the melting point (∼760 K), features in D(r) broaden upon heating.

Fig. 1 Temperature-dependent PDF results. (a) X-ray reduced pair distribution functions, D(r), for the KCl–MgCl₂ salt mixture (50 : 50 mol%). (b) KMgCl₃ (Pnma) and CrCl₃ (P3₁12) crystal structures with highlighted Mg–Cl and Cr–Cl average bond lengths in Å. Color scheme: K, red; Mg, blue; Cr, purple; Cl, green. Vertical dashed lines in (a) are guides to the eye.

Based on the relative X-ray weighting factors (ESI Fig. 3†) and the available orthorhombic KMgCl₃ crystal structure (Fig. 1b) corresponding to the room temperature phase of KCl–MgCl₂, the first peak in the D(r)s is primarily dominated by short-range Mg–Cl ionic correlations. The experimentally determined Mg–Cl average bond length contracts on heating from 2.50 Å at 298 K to 2.47 Å at 673 K, achieving the shortest distance of 2.37 Å in the molten state, whereas longer-range order correlations show the opposite trend consistent with thermal expansion (Fig. 1a), and at 1073 K they significantly broaden or disappear as expected in the molten state. What may at first seem counterintuitive, the observed shrinking of the cation–anion close-contact distances is quite typical for inorganic salts and is associated with a decrease in the coordination numbers on melting, often leading to shorter nearest-neighbour bond lengths than those observed in the solid state.^25,28–30 The characteristically broad peaks in D(r) at the highest temperature associated with the liquid state make the interpretation of scattering results difficult, requiring more sophisticated procedures in the data analysis.

To emulate the presence of dilute metal ion species in the molten ionic medium, we added a small amount of CrCl₃ (5 mol%) to the KCl–MgCl₂ salt mixture. Because the concentration of Cr³⁺ in the sample is low, information in S(Q) and the total PDF about its interactions with other species in the matrix will be masked by the KCl–MgCl₂ solvent. This problem is compounded by the fact that interionic distances between Cr³⁺ and other species can be quite similar to those between solvent species. Therefore, in order to capture the small variation in S(Q) and G(r) caused by the presence of dilute Cr³⁺ species, we performed high-energy X-ray diffraction measurements for the molten salt mixture before and after addition of CrCl₃ (Fig. 2).

Fig. 2 X-ray scattering measurements of the KCl–MgCl₂ (black) and CrCl₃–KCl–MgCl₂ (red) molten salt mixtures at 1073 K. (a) X-ray structure functions S(Q). (b) X-ray pair distribution functions G(r) generated from the S(Q) patterns in (a), and the corresponding differential PDF dG(r) (blue), which is offset for clarity. Circles represent the experimental measurements; solid lines are from the RMC modeling (due to a nearly perfect match the discrepancy between the measured datapoints and the RMC model is smaller than the line thickness); dashed lines are from the AIMD simulations. The inset shows a snapshot from the RMC modeling highlighting ionic clustering of Cr³⁺ species (enlarged purple spheres) and representative chromium chloride dimers present in the melt.

As expected, the obtained S(Q)s and G(r)s show very subtle differences for the CrCl₃-loaded and pristine salt systems (Fig. 2). Close examination reveals an enhancement of the first peak intensity at ∼2.37 Å, suggesting that Cr–Cl first neighbor distances nearly coincide with those of Mg–Cl in the melt. These results also imply that Cr³⁺ is surrounded by an octahedral configuration of Cl⁻ anions, since the Cr–Cl bond lengths do not seem to change in the molten state as compared to the crystalline form of octahedral CrCl₃ in Fig. 1b. This evidence is further supported by our high-temperature optical absorption spectroscopy results (ESI Fig. 4†). The absorption spectrum consists of two bands in the visible region that correspond to the ⁴A₂ → ⁴T₂ (⁴F) and ⁴A₂ → ⁴T₁ (⁴F) transitions for the peaks at 12 000 cm⁻¹ and 18 000 cm⁻¹, respectively, as well as a charge transfer band corresponding to the ⁴A₂ → ⁴T₁ (⁴P) transition, consistent with the spectrum of dilute Cr³⁺ in the LiCl–KCl eutectic melt where octahedral coordination was previously identified.¹⁴

The differential PDF, dG(r), obtained by direct subtraction of the reference PDF of the KCl–MgCl₂ molten salt from that of the CrCl₃-loaded mixture is shown in Fig. 2b. We performed AIMD simulations to interpret the differential PDF results. Fig. 2a shows very good agreement between the experimental and simulated S(Q) beyond 2 Å⁻¹. At lower Q values, a good match is not expected due to the limited size of our AIMD boxes. RMC simulations, starting from the AIMD snapshot configuration where the box size has been quadrupled, completely resolved this discrepancy, providing a perfect match between the experiment and the RMC model (Fig. 2). The comparison of radial distribution functions, g(r), in ESI Fig. 5† indicates a good agreement between the results from RMC and AIMD simulations, despite the limitation of the latter technique to reproduce the S(Q) features at low Q. This provides some confidence in the analysis of the dPDF at short distances based on the first principles results. In addition, close agreement of the AIMD simulations with the experimental X-ray G(r) allows full isolation of all ionic correlations (Cr–Cl, Cr–Cr, Cr–Mg, Cr–K, etc.) via the Fourier transformation of the partial subcomponents of S(Q) (ESI Fig. 6†). This can also provide valuable information as to which pair interactions contribute to the differential PDF at specific distances (ESI Fig. 7†). It is worth mentioning, however, that while computational G(r) functions reproduce the experiment well, the match in the differential PDF is apparently worse, highlighting the complexity of extracting data from the differential technique, which makes use of small changes in the overall signal upon introduction of dilute Cr³⁺ species.

Fig. 3a shows the partial subcomponents of dG(r) as defined from the subtraction of partial subcomponents from the two samples in the first principles simulations. One may see that the first and most prominent peak in dG(r) mainly results from two competing contributions, a positive peak from Cr–Cl and a negative peak at similar distances from Mg–Cl correlations. This positive peak is very simple to interpret: in a differential PDF calculation, terms that only exist in the CrCl₃-loaded sample are equal to the corresponding partial subcomponent of G(r), i.e. dG_Cr–Cl(r) = G_Cr–Cl(r). In other words, the first peak in dG(r) shows the location of the closest contact between Cr³⁺ and Cl⁻.

Fig. 3 Interpretation of the differential PDF. (a) Partial subcomponents of dG(r) from AIMD simulations. (b) Differential coordination number d(CN) for the solvent–solvent correlations, showing the difference between the cumulative coordination number in the CrCl₃ loaded system and the pristine KCl–MgCl₂ molten salt.

Instead, solvent–solvent type terms that appear in both systems, and are of the same order of magnitude, require more nuance in our analysis and interpretation. Fig. 3b displays the differential coordination number, d(CN), of selected ions. A positive d(CN) means that the cumulative sum of ions of certain class up to a distance r surrounding a central ion is larger in the Cr³⁺-containing sample, whereas a negative differential coordination number means that this quantity is smaller. We see for example, that whereas the depletion in the case of dG_Mg–Cl(r) below 2.5 Å is consistent with a depletion in the coordination number up to that distance (and also with a depletion in the actual computed radial distribution function g_Mg–Cl(r) at that distance (ESI Fig. 8†)), this is not the case for the K–Cl correlation at ∼3 Å. Here, consistent with the behavior of dG_K–Cl(r), the first peak in g(r) for this ionic pair is slightly lower for the system with CrCl₃, but the coordination number of Cl⁻ around K⁺ is larger (Fig. 3b). In other words, a negative/positive peak in the partial subcomponents of dG(r) for solvent–solvent components, cannot be simply interpreted as a depletion/increase in coordination. This is because the addition of a solute (CrCl₃) results in a concomitant change in the molar volume and hence in the number density of the different species. This phenomenon is particularly complicated in molten salts where a solute, in our case Cr³⁺, necessarily comes combined with Cl⁻ counter ions that are both part of the solute and the solvent. Instead, the experimental observable dG(r) depends on mole fractions of the ions and their effective number of electrons.

The second peak in dG(r) is dominated by Cl–Cl interactions between 3 and 3.5 Å. In this case, the enhancement is both in dG_Cl–Cl(r) and the coordination number (Fig. 3). It mostly represents a shift in the g_Cl–Cl(r) towards shorter distances (ESI Fig. 8†) due to the strong Cr–Cl attractive interactions when dilute Cr³⁺ ions are present in the melt. This behavior is consistent with the Cl–Cl contraction observed recently in neutron diffraction experiments for the binary CrCl₃–NaCl melt with much higher CrCl₃ content (22 mol%).³¹

A perhaps unexpected result from our PDF fitting using RMC is that Cr³⁺ ions tend to aggregate, mainly forming [Cr₂Cl₁₀]⁴⁻ dimeric species in the KCl–MgCl₂ molten salt even at low concentration and high temperature (see inset in Fig. 2b). Our AIMD simulations, starting from the initial configuration where Cr³⁺ ions were deliberately placed far away from each other, confirm Cr³⁺ clustering after about 10 ps. To gain additional insight into the observed ion clustering and coordination environments, we explored free energy landscapes that describe the distribution and relative stability of various coordination structures in the melt. As can be seen from the free energy profiles in Fig. 4, Cr³⁺ can be surrounded by 4–6 chloride anions in the first solvation shell, with the octahedral coordination clearly being the most favourable. The 5- and 4-coordinated states are rarely accessible due to rather high free energy barriers. This is consistent with our optical spectroscopy measurements, showing the persistence of spectral features and hence the octahedral geometry of the Cr³⁺ species over a wide range of temperatures (ESI Fig. 4†). In contrast, the 5- and 4-coordinated states of the solvent Mg²⁺ ions are almost equally probable. The free energy barrier between them is less than 1 kcal mol⁻¹, which is well below the thermal energy at 1073 K (2.13 kcal mol⁻¹) affording frequent interconversions of these states. Unlike highly charged Cr³⁺ and Mg²⁺, K⁺ ions exhibit a single 9-coordinated state, i.e. 9 Cl⁻ around a K⁺ ion (Fig. 4). Addition of Cr³⁺ ions slightly modifies the free energy profiles for the solvent cations. The 4-coordinated state becomes less stable for Mg²⁺, indicating an increase in the average coordination number. This is in agreement with our results in Fig. 3b, showing positive d(CN) values for Mg–Cl at the cutoff distance (3.5 Å) which we use to define the average coordination number for Mg²⁺. The free energy profile for K⁺ shifts slightly to a higher coordination number as well (Fig. 4), consistent with the corresponding d(CN) for K–Cl. This can be likely attributed to the high charge density of Cr³⁺ capable of polarizing Cl⁻ significantly, which results in weaker Cl–Cl repulsive interactions and thus allows more Cl⁻ ions in the first coordination shell of the cations.

Fig. 4 Free energy as a function of coordination number of Cr³⁺, Mg²⁺, and K⁺. Solid and dashed lines indicate the KCl–MgCl₂ salt mixtures with and without Cr³⁺, respectively.

The polarization effect caused by Cr³⁺ ions is better understood by analyzing the free energy surfaces for the chloride shared cation–cation configurations. The 2D-free energy surfaces in Fig. 5a show that the coordination environment of Cr³⁺ remains octahedral regardless of whether it makes close contact with different Cr³⁺ or Mg²⁺ ions. However, Mg²⁺ adopts a slightly higher coordination number (between 5 and 6) when contacting with Cr³⁺, while switching back and forth between the 4- and 5-coordinated states as it moves away from the highly charged Cr³⁺. In other words, Mg²⁺ draws more Cl⁻ ions in its first solvation shell when forming a dimeric species with Cr³⁺. The free energy surfaces in Fig. 5b demonstrate that Cr³⁺ and Mg²⁺ ions in their most probable [Cr₂Cl₁₀]⁴⁻ and [CrMgCl₉]⁴⁻ dimeric configurations (Fig. 5d) are linked together through two bridging Cl⁻ anions. In addition, the existence of chloride-shared dimeric Cr–K configurations is rare according to our calculations (Fig. 5b), hence stable dimers can only be formed between multivalent cations.

Fig. 5 2D-free energy surfaces highlighting (a) most probable coordination number of Cr³⁺ (CN_Cr) and Mg²⁺ (CN_Mg) across a range of Cr–Cr and Cr–Mg distances; (b) the chloride-shared and chloride-separated states for Cr–Cr, Cr–Mg, and Cr–K dimeric configurations; (c) number of Cr³⁺ ions coordinating with different Cr³⁺ and Mg²⁺ ions, and number of Mg²⁺ ions coordinating with Cr³⁺ and different Mg²⁺ ions via chlorides, with the respective dimers, trimers, and tetramers pointing to the region where they form; (d) representative chloride-shared [Cr₂Cl₁₀]⁴⁻ and [CrMgCl₉]⁴⁻ dimeric configurations present in the CrCl₃–KCl–MgCl₂ melt. Color bar represents the contour levels in the 2D-free energy surfaces.

One may expect that the high charge density of Cr³⁺ can also induce the formation of more complex aggregates, such as trimers or even longer chains. However, only the Cr–Cr dimeric species were found in the melt, likely due to the low concentration of Cr³⁺ and high temperature of the system. The 2D-free energy surface in Fig. 5c (left panel) justifies this, showing that even the formation of Cr–Cr–Cr trimers is highly improbable. The same conclusion can be drawn for Cr–Mg–Cr trimers, which are unlikely to exist (Fig. 5c, middle panel). However, free energetically metastable Cr–Mg–Mg trimers and Cr–Mg–Mg–Mg tetramers can still form (Fig. 5c, right panel) due to the relatively high concentration of Mg²⁺ ions in the KCl–MgCl₂ molten salt. Whereas this analysis can only reveal correlations that are shorter-range, based on our experience studying KCl–MgCl₂ mixtures,²⁵ it is likely that observations in Fig. 5c (right panel) imply that dilute Cr³⁺ is becoming a component in Mg²⁺ networks.

In addition to the structural diversity, the dynamics of ions with different charge density can be highly heterogeneous in the multicomponent molten salt system. To examine this, we analyzed the mean square displacement (MSD) of each ion over time. Fig. 6 shows that the large ions with low charge density, K⁺ and Cl⁻, move at faster rates than the smaller ions with high charge density, Mg²⁺ and Cr³⁺. Typically, the MSD analysis can provide diffusion coefficients, but it requires information about long-time displacements that is lacking in our 40 ps AIMD simulation time. Nevertheless, the short-time dynamics of these ions are still critically useful and can be correlated with ion exchange and clustering kinetics.^32,33 We therefore study these processes around the solute Cr³⁺, the solvent Mg²⁺ and K⁺ cations using a hybrid TS-Marcus theory framework, which allows us to elucidate the origin of heterogeneous local dynamics in our multicomponent ionic system.

Fig. 6 Mean square displacement (MSD) of ions as a function of time, indicating faster dynamics in the order: K⁺ > Cl⁻ > Mg²⁺ > Cr³⁺. Dots represent the MSD data and solid lines are linear fits.

Conventional rate theory-based studies of ion pairing and solvation typically utilized the interionic/intermolecular distance as the reaction coordinate to determine the transition rates between the reactant and product states.^34–37 However, recent studies indicate that inclusion of a second reaction coordinate associated with solvent fluctuations is necessary to provide a better understanding of the underlying mechanisms and rates.^38–41 These studies show that the fluctuating coordination number or the electric field are both reasonable collective variables that can be used to represent solvent fluctuations. We find the electric field on a solute to be a better reaction coordinate because it is determined by considering the contribution from the entire solvent, whereas the coordination number only accounts for the local solvent environment. Herein, considering that the reactant and product states can be resolved as well-separated states in the electric field space, we formulate a hybrid TS-Marcus framework (see ESI† for more details) to determine the reactant-to-product transition rates describing ion exchange and clustering processes.

To investigate local dynamics around Cr³⁺, we focused on three specific processes: chloride exchange around Cr³⁺ and the dissociation of Cr–Cr and Cr–Mg dimers. For the case of chloride exchange, we obtained the free energy surface as a function of the Cr–Cl distance and the electric field that Cl⁻ experiences along the direction (Fig. 7a). Analogously for the cases of Cr–Cr (Fig. 7b) or Cr–Mg (Fig. 7c) dissociation, the free energy is for the Cr–Cr or Cr–Mg distance and the electric field that the Cr³⁺ or Mg²⁺ ion experiences along the or direction. When Cl⁻ is located within the first coordination shell of Cr³⁺ representing the reactant state of the chloride exchange process (Fig. 8a), the Cl⁻ ion experiences a strong field. However, as soon as another Cl⁻ moves from the bulk solution into the Cr³⁺ solvation shell driven by solvent rearrangement, the local electrostatic environment changes in a way to reduce the positive field experienced by the chromium bound Cl⁻ ion (the intermediate state in Fig. 8a). Further weakening of the Cr–Cl interaction leads to the transition toward the product equilibrium, in which the initially bound Cl⁻ completely escaped the first coordination sphere of Cr³⁺. In the product state, the Cl⁻ ion experiences a broad range of negative and positive fields distributed around the mean value of ∼0. This is due to the possibility of Cl⁻ interacting with all other ions in an isotropic environment.

Fig. 7 2D free energy surfaces, involving Cr–Cl, Cr–Cr, Cr–Mg distance and electric field, E, experienced by (a) Cl, (b) Cr, and (c) Mg with corresponding Marcus parabolas (solid lines are parabolic fits and dots are actual data), describing chloride exchange around Cr³⁺ and dissociation of the Cr–Cr and Cr–Mg dimers. The lines with arrows are the Marcus pathways showing how reorganization of ionic media causes electric field rearrangement and drives the exchange and dissociation events. (d) The corresponding 1D free energy profiles.

Fig. 8 Reactant, intermediate, and product states for (a) chloride exchange and (b) dissociation of dimeric complexes, indicating that the solvent rearrangement-driven motion of Cl⁻ into the solvation shell changes the local electrostatic environment, triggering the exchange/dissociation process. This change in the electrostatic environment is gleaned from (c) the lower and higher adiabatic energy surfaces (solid lines) obtained from coupled Marcus diabats (dots) in the electric field space. The reactant-to-product transition path is shown by a green arrow.

The discussed pathway of chloride exchange is depicted in Fig. 7a using an arrow on the free energy surface. For the dissociation of Cr–Cr and Cr–Mg dimers, the underlying mechanism can be explained using a similar approach as illustrated in Fig. 7b, c and 8b. In these cases, the field in the reactant state is already week due to sharing Cl⁻ ions between the two cations in the dimeric configurations, leading to the partial screening of the Cr–Cr and Cr–Mg interactions. After a full separation in the product state, their interaction strength vanishes completely.

It is worth noting that the minima representing the equilibrium close-contact (r_R) and dissociated (r_P) distances on the 2D-surfaces are identical to those on the 1D-free energy profiles along the interionic distance (Fig. 7d). To perform a Marcus theory-based analysis, we curved out a couple of slices from the 2D-free energy surfaces along r = r_R and r = r_P resulting in the reactant and product diabats in the electric field space, which we modeled as parabolas (Fig. 7 and ESI Table 1†). The coupling between these diabats generates lower and higher adiabatic free energy surfaces shown in Fig. 8c. From the perspective of Marcus theory, we find that while chloride exchange around Cr³⁺ can be described as the reactant (strong field)-to-product (near zero field) transition through the crossing region (identified as “normal region”) on the lower adiabatic surface, a nonadiabatic jump from the lower to the higher adiabatic surface through the crossing region (identified as “inverted region”) is required to drive the dissociation of the Cr–Cr or Cr–Mg dimeric structures. The crossing region corresponds to the respective intermediate configuration depicted in Fig. 8a and b, where the field gets weakened, subsequently triggering the transition to the product state.

To discern between the local solvent dynamics around the Cr³⁺ ion and the overall solvent dynamics in the molten salt system, we explored 2D-free energy surfaces and their respective Marcus diabatic and adiabatic states, and examined chloride exchange around solvent Mg²⁺ and K⁺ ions as well as dissociation of the Mg–Mg dimers (ESI Fig. 9 and 10†). While chloride exchange around Mg²⁺ is an adiabatic process, chloride exchange around K⁺ and Mg–Mg dissociation are nonadiabatic processes. Considering all the processes involving both solute and solvent ions, we find a trend that when an ion experiences a strong field in a close-contact configuration, the reactant and the product states are well-separated in the electric field space and cross at the “normal region”, exhibiting adiabatic interconversion. On the other hand, a weak field on an ion at a close-contact configuration does not separate the reactant and product states well and is responsible for these states crossing at the “inverted” region. Thus, the nature of these dynamical processes is highly sensitive to the charge density of ions in the molten salt.

There are several different factors that govern the kinetics of chloride exchange and dissociation of the dimers at a given temperature. As discussed in detail in the ESI,† the reactant-to-product transition rate (given by eqn (1)) depends on the free energy barrier (W_tot), the transmission coefficients (κ_LZ) representing barrier recrossing events, and the mass-weighted “reactant volume” at the given thermal energy, 1/β.

Of all these parameters, we find that the most important ones are the free energy barriers and transmission coefficients that change much more drastically than the mass-weighted “reactant volume” as we go across different cases (ESI Table 1†). The trend in the total barriers is consistent with the strength of interionic electrostatics in the close-contact reactant state. For instance, chloride exchange around Cr³⁺ has the highest barrier of ∼10 kcal mol⁻¹, while Mg–Mg dissociation has to overcome a very low barrier of ∼2 kcal mol⁻¹. Because of the difference in barrier heights, if we consider no barrier recrossing chloride exchange around Cr³⁺ can be slower than for Cr–Cr and Cr–Mg dissociation (for no or minimal barrier-recrossing, κ_LZ ∼ 1, while κ_LZ ≪ 1 for a significantly large number of recrossing events). However, when barrier recrossing is taken into account, Cr–Cr separation and chloride exchange around Cr³⁺ appear to occur on the same timescale of about 32–35 ps. This is because of the low transmission coefficient (0.116) for the separation of Cr–Cr dimers. This low κ_LZ implies that dissociation of these species is normally promptly followed by recombination. This inseparability of the time scales associated with chloride exchange around Cr³⁺ and the Cr–Cr dissociation implies that dilute chromium species may exist in the melt as kinetically stable [Cr₂Cl₁₀]⁴⁻ clusters for several tens of picoseconds. In contrast, the lifetime of the Cr–Mg dimers is much shorter (∼4 ps) than the timescales for chloride exchange around Cr³⁺ and Mg²⁺ ions. This is because of the much smaller barrier and larger transmission coefficient. The dynamics of the solvent ions follow a similar trend, with the barrier height making chloride exchange around Mg²⁺ slower than chloride exchange around K⁺ and Mg–Mg dissociation. In fact, without considering barrier recrossing, the latter two processes exhibit the fastest dynamics in the molten salt system at the same timescale of 0.5 ps. Interestingly, the transmission coefficient for Mg–Mg dissociation is extremely small (κ_LZ ∼ 0.009), suggesting that rapid dissociation of Mg–Mg pairs is immediately followed by recombination providing dynamical stability of the –Mg–Mg– chains.

Conclusions

The dynamic structural behaviour of dilute Cr³⁺ species in a high-temperature ionic chloride medium was investigated by means of high-energy X-ray scattering and AIMD simulations. In the molten KCl–MgCl₂ mixture, the first solvation shell of Cr³⁺ comprises 6 Cl⁻ anions, with computational results indicating persistence of the octahedral environment at least up to 1073 K. Therefore, the six-fold coordination of Cr³⁺, typically observed in the solid state, also occurs in the liquid chloride mixture, which can be explained by the high crystal field stabilization energy for the octahedral Cr³⁺ complex.⁴² Beyond simply identifying the Cr³⁺ coordination environment, our findings challenge the paradigm of molten salts being a “sea of cations in a sea of anions,” where aggregation of dilute species is not expected due to high temperature conditions. Instead, the results of our studies indicate the formation of dinuclear chloride-shared Cr–Cr clusters with a relatively long lifetime in the melt. These fundamental insights into the putative mechanisms of ion exchange and clustering were only accessible when our X-ray scattering experiments and AIMD simulations were put into context using the framework of a hybrid TS-Marcus theory, highlighting the importance of understanding the structure, dynamics, and energetics of the molten salt composition at the atomic level. Furthermore, our analysis shows that Cr³⁺ species possessing high charge density exhibit the slowest short-time dynamics among all constituent ions in the melt. Hence, during Cr dealloying/corrosion it is possible that concentration gradients will occur at the alloy-molten salt interface resulting in Cr³⁺-containing clusters, like those observed in this study, or even larger networks. These structures may be significantly different from those in the bulk and so will be the local chemical potentials conceivably affecting the dealloying rate. We expect that results presented here stemming from multiple experimental, theoretical and computational techniques will foster discussion on the issue of kinetics and aggregation in purely ionic media, with pertinence to novel molten salt reactors and emerging concentrated solar power technologies.

Experimental section

Sample preparation

The MgCl₂ (anhydrous, >98%) was purchased from Sigma-Aldrich and further purified by fractional distillation described previously.^25,29 The KCl (anhydrous, 99.999%) and CrCl₃ (anhydrous, 99.9%) were both purchased from Sigma-Aldrich and used for experiments as received. Great care was taken so that the salts were maintained under a vacuum or inert atmosphere for the duration of sample preparation. Distilled MgCl₂ and KCl were both added to a fused silica test tube inside an inert atmosphere glovebox (<0.5 ppm O₂) at a 50 : 50 molar ratio. A vacuum valve compression fitting was then attached to the open end of the tube and the valve closed. This mixture was taken outside of the glovebox, evacuated to a pressure of 1.0 × 10⁻³ torr, and placed in a furnace at a temperature of 850 °C. Dynamic vacuum and temperature were maintained for 15 minutes to fuse the MgCl₂–KCl mixture. The fused mixture was then cooled to room temperature under dynamic vacuum, sealed, and returned to the glovebox. Proper masses of both the MgCl₂–KCl and CrCl₃ salts were mixed to obtain the corresponding 5.0 mol% CrCl₃ in the 50 : 50 MgCl₂–KCl. This mixture was then fused together under dynamic vacuum following the same procedure as described for the MgCl₂–KCl. The obtained salt compositions were then crushed to fine powder, added to separate thin-walled quartz capillaries (Charles Supper Co., 1.5 mm O.D., 0.010 mm wall thickness), evacuated to a pressure of 1.0 × 10⁻³ torr, and flame sealed. Additional details of the density measurements are provided in the ESI.†

Synchrotron X-ray scattering experiments

Collection of total X-ray scattering data for PDF analysis took place at the 28-ID-1 beamline of the National Synchrotron Light Source II (NSLS-II), Brookhaven National Laboratory. A horizontally focusing side bounce monochromator was used to deliver X-rays with the energy of 74.4 keV (0.1667 Å) and a beam of cross-sectional area 0.25 × 0.25 mm². An amorphous 2D silicon-based area detector (PerkinElmer XRD 1621, 200 × 200 micron pixels) was positioned about 20.8 cm from the sample to collect total scattering intensity data, which yielded a useable Q range up to 24 Å⁻¹. A Ni powder standard was used to calibrate the sample to detector distance, beam center, and detector tilt and rotation. Measurements were performed in transmission mode using a customized furnace (ESI Fig. 2†) described in our previous study.²⁵ The temperature was set using an Omega Engineering CN7800 series temperature controller. An Omega stainless-steel-sheathed K-type thermocouple was used to control a 24 volt DC TDK Lambda power supply. The furnace was temperature-calibrated using a K-type thermocouple inserted into a quartz tube (Charles Supper Co., 1.5 mm O.D., 0.010 mm wall thickness).

The X-ray scattering data were collected for the duration of 30 min for each sample and empty capillary. To minimize the overexcitation of pixels in the 2D detector, the duration of beam exposure was shorter for the crystalline samples than for the liquid salts at 1073 K. In addition, dark current scans were collected and subtracted from the raw X-ray patterns, allowing to remove residual intensity and prevent pixel overexcitation.⁴³ After data collection, the GSAS-II program⁴⁴ was used to radially integrate the 2D detector images to 1D diffraction patterns. The appropriate masks were applied to remove dead pixels and the beamstop from data integration. PDFgetX2 (ref. 45) was then used to subtract the background signal (empty quartz capillary and air scattering) and to perform additional corrections (sample self-absorption, multiple scattering, and inelastic Compton scattering) following standard procedures.^26,27 The obtained data were subsequently normalized to the average electron density given by the weighted sum of the atomic X-ray form factors, yielding the total structure function, S(Q), which was then converted to PDF via Fourier transformation using a Q range of 0.5–18 Å⁻¹ (see the ESI† for the corresponding equations).

Optical spectroscopy measurements

All salt handling operations and optical spectroscopy measurements were performed inside an Ar atmosphere glovebox maintained at less than 0.1 ppm moisture and O₂ concentrations. The absorption spectra were measured with the Cary 5000 spectrophotometer, equipped with deuterium and tungsten lamp covering the ultraviolet-visible-near infrared region over the wavelength range of 2300 to 200 nm. High temperature rated fiber optic cables procured from Fiberguide Inc. were used to couple the spectrophotometer (located external to the glovebox) to the custom designed furnace located outside the glovebox. To record the absorption spectra in transmission mode, molten salt samples were contained in a quartz cuvette (optical path length 1 cm, FireflySci) placed on a multi-cell rack designed with graphite; the furnace setup is described in our previous study.⁴⁶ The wavelength scale of the spectrophotometer over the range of 880 to 280 nm was verified and calibrated using a NIST certified reference material consisting of an aqueous solution of didymium perchlorate that is permanently sealed by heat fusion in a high quality Far UV quartz cell. The salt samples were loaded in cleaned and dried quartz cuvettes. The total mass for each sample was approximately 4 g. A blank sample of base salt composition was used to calibrate the baseline of the instrument prior to measurement of the spectra at each temperature setup point.

Raman spectroscopy experiments

High temperature Raman studies were carried out using a Horiba T64000 spectrometer. A diode laser of 532 nm was focused through a 50× objective of the microscope. The power of the laser was 100 mW. A Linkam heating stage TMS1100 was used for measuring Raman signal of molten salts at 1073 K. The stage was connected to a nitrogen flow (15 ml min⁻¹) to keep the samples (50–80 mg) in alumina crucibles under an inert atmosphere during the measurements. The acquisition time was 10 s for each sample. The Raman shift was measured from 100 to 1200 cm⁻¹, the range that should cover the expected Raman signal of MgCl₂–KCl molten salt mixtures.^47,48

Ab initio molecular dynamics simulations

The initial structures for AIMD simulations were first pre-equilibrated by running force field-based simulations using the polarizable ion (PIM) model. The initial configuration for 50% MgCl₂ + 50% KCl without CrCl₃ was taken from our previous AIMD simulation that employed a different density functional.²⁵ We used the average equilibrium density from the PIM model (ρ = 1.57 g cm⁻³) that is in close agreement with the experimental value (ρ = 1.56 g cm⁻³).⁴⁹ Since there are no published PIM parameters for CrCl₃, we constructed an initial system containing AlCl₃, more specifically 6 Al³⁺, 60 Mg²⁺, 60 K⁺, and 198 Cl⁻ ions. We performed force field-based simulations of this system using the PIM model^25,50 in the constant temperature–pressure ensemble, wherein we first heated the system at 1600 K at 1 bar for 100 ps and then equilibrated it at 1073 K keeping the same pressure for 500 ps. We extracted a configuration from this trajectory that has the experimental density of 1.65 g cm⁻³ (a total of 324 ions in a box with length of 22.3812 Å) and replaced Al with Cr before starting the AIMD simulations.

AIMD simulations were carried out using the Quickstep module of CP2K 6.1 (ref. 51) utilizing the PBE exchange–correlation functional^52–54 in conjunction with Grimme's D3 dispersion correction.⁵⁵ The MOLOPT-DZVP basis set⁵⁶ in conjunction with Goedecker–Teter–Hutter (GTH) pseudopotentials⁵⁷ was applied for all ions, together with a density cutoff of 600 Ry. All DFT calculations were performed with the Hubbard + U term (U_eff = U − J = 3.36 eV). Simulations were carried out using a time step of 0.5 fs and the temperature was set to 1073 K using a Nosé–Hoover thermostat^58,59 with a velocity rescaling time constant of 1.0 ps. We obtained a 55 ps long AIMD trajectory and used the last 40 ps for the structural and dynamical analysis. Spin unrestricted calculations with several options for initial magnetizations were performed before the MD step to ensure that each Cr³⁺ has the expected spin state and the ground electronic state is identified. The orbital transformation method with a FULL_ALL preconditioner and a conjugate gradient minimizer was employed for achieving and accelerating the SCF convergence. Details of S(Q), G(r), and dG(r) calculations based on the AIMD trajectory are given in the ESI.†

Reverse Monte Carlo modelling

The X-ray total scattering pattern was fitted through the RMC approach with the AIMD cell expanded to a 4 × 4 × 4 supercell, targeting an optimal match with experimental data in a metropolis manner. In our study, the RMCProfile package^60,61 with such an algorithm implemented was employed to fit the X-ray total scattering data in both real and reciprocal space. A total of 20 736 ions in a cubic box with length of 89.5248 Å was used for the RMC fitting, with ∼493 moves per ion generated and ∼56 moves per ions accepted, overall. Apart from the experimental data, only the minimum distance for each ionic pair was used as constraint to allow the system to fully relax. To guarantee the generality of the obtained RMC results, the fit was performed using 50 different configurations and an average was taken over all RMC runs to obtain partial radial distribution functions, which were compared to our AIMD results.

Author contributions

Project administration, A. S. I., C. J. M., V. S. B. Conceptualization and design of experiments, A. S. I. X-ray scattering experiments, A. S. I., M. A., S. K. G. and analysis, A. S. I. Raman measurements, A. S. I., D. S. M. Optical spectroscopy measurements, R. G. Sample preparation, P. H., D. S. M., S. M. M., S. D. AIMD simulations, V. S. B. Structure function, PDF, dPDF computation and analysis, C. J. M., S. S., W. V. K., F. W. RMC modelling, A. S. I., Y. Z. TS-Marcus theory model conceptualization, rate theory calculations and clustering analysis, S. R. Writing – original draft, A. S. I., C. J. M., S. R. All authors proofread, commented on, and approved the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported as part of the Molten Salts in Extreme Environments (MSEE) Energy Frontier Research Center, funded by the U.S. Department of Energy Office of Science. MSEE work at Iowa was supported via subcontract from Brookhaven National Laboratory (BNL). BNL and ORNL operate under DOE contracts DE-SC0012704 and DE-AC05-00OR22725, respectively. The project used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC05-00OR22725. The 28-ID-1 beamline of the National Synchrotron Light Source II was used, which is a U.S. Department of Energy (DOE) Office of Science User Facility operated for the DOE Office of Science by Brookhaven National Laboratory under Contract No. DE-SC0012704. S. S., W. V. K., F. W. and C. J. M. acknowledge the computational resources at the University of Iowa high performance computing facility. The authors thank Dr James F. Wishart (BNL) for critically reading the manuscript and offering helpful suggestions. This manuscript has been authored in part by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Notes and references

1. L. A. Weinstein, J. Loomis, B. Bhatia, D. M. Bierman, E. N. Wang and G. Chen, Chem. Rev., 2015, 115, 12797–12838.
2. C. D. Ferguson, L. E. Marburger, J. D. Farmer and A. Makhijani, Nature, 2010, 467, 391–393.
3. D. LeBlanc, Nucl. Eng. Des., 2010, 240, 1644–1656.
4. R. G. Reddy, J. Phase Equilib. Diffus., 2011, 32, 269–270.
5. J. Bulmer, Fused salt fast breeder: reactor design and feasibility study, United States Atomic Energy Commission, Technical Information Service Extension, 1957.
6. D. J. Diamond, N. R. Brown, R. Denning and S. Bajorek, Phenomena important in molten salt reactor simulations, Brookhaven National Lab.(BNL), Upton, NY (United States), 2018.
7. C. K. Ho, Sol. Energy, 2017, 152, 38–56.
8. J. C. Vidal and N. Klammer, AIP Conf. Proc., 2019, 2126, 080006.
9. X. Xu, X. Wang, P. Li, Y. Li, Q. Hao, B. Xiao, H. Elsentriecy and D. Gervasio, J. Sol. Energy Eng., 2018, 140, 051011.
10. J. C. Gomez-Vidal, npj Mater. Degrad., 2017, 1, 1–9.
11. J. M. Kurley, P. W. Halstenberg, A. McAlister, S. Raiman, S. Dai and R. T. Mayes, RSC Adv., 2019, 9, 25602–25608.
12. A. Ronne, L. He, D. Dolzhnikov, Y. Xie, M. Ge, P. Halstenberg, Y. Wang, B. T. Manard, X. Xiao and W.-K. Lee, ACS Appl. Mater. Interfaces, 2020, 12, 17321–17333.
13. H. Brookes and S. Flengas, Can. J. Chem., 1970, 48, 55–58.
14. D. Gruen and R. L. McBeth, Pure Appl. Chem., 1963, 6, 23–48.
15. E. T. Dias, S. K. Gill, Y. Liu, P. Halstenberg, S. Dai, J. Huang, J. Mausz, R. Gakhar, W. C. Phillips, S. Mahurin, S. M. Pimblott, J. F. Wishart and A. I. Frenkel, J. Phys. Chem. Lett., 2021, 12, 157–164.
16. S. K. Gill, J. Huang, J. Mausz, R. Gakhar, S. Roy, F. Vila, M. Topsakal, W. C. Phillips, B. Layne, S. Mahurin, P. Halstenberg, S. Dai, J. F. Wishart, V. S. Bryantsev and A. I. Frenkel, J. Phys. Chem. B, 2020, 124, 1253–1258.
17. C. Hardacre, Annu. Rev. Mater. Res., 2005, 35, 29–49.
18. A. L. Smith, M. N. Verleg, J. Vlieland, D. de Haas, J. A. Ocadiz-Flores, P. Martin, J. Rothe, K. Dardenne, M. Salanne, A. E. Gheribi, E. Capelli, L. van Eijck and R. J. M. Konings, J. Synchrotron Radiat., 2019, 26, 124–136.
19. J. Moon, C. Abney, D. Dolzhnikov, J. M. Kurley, K. A. Beyer, J. Wright, S. Dai, R. T. Mayes, S. S. Raiman and J. Wright, Chemrxiv, 2019 DOI:10.26434/chemrxiv.9959807.
20. L. Soderholm, S. Skanthakumar and R. E. Wilson, J. Phys. Chem. A, 2011, 115, 4959–4967.
21. H. Deguchi, Y. Kubota, H. Furukawa, Y. Yagi, Y. Imai, M. Tatsumi, N. Yamazaki, N. Watari, T. Hirata, N. Matubayasi and Y. Kameda, Int. J. Greenhouse Gas Control, 2011, 5, 1533–1539.
22. V. Kamysbayev, V. Srivastava, N. B. Ludwig, O. J. Borkiewicz, H. Zhang, J. Ilavsky, B. Lee, K. W. Chapman, S. Vaikuntanathan and D. V. Talapin, ACS Nano, 2019, 13, 5760–5770.
23. S. L. Thomä, S. W. Krauss, M. Eckardt, P. Chater and M. Zobel, Nat. Commun., 2019, 10, 1–7.
24. M. Zobel, R. B. Neder and S. A. Kimber, Science, 2015, 347, 292–294.
25. 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.
26. T. Egami and S. J. Billinge, Underneath the Bragg peaks: structural analysis of complex materials, Newnes, 2012.
27. H. E. Fischer, A. C. Barnes and P. S. Salmon, Rep. Prog. Phys., 2005, 69, 233.
28. C. Angell and D. Gruen, J. Phys. Chem., 1966, 70, 1601–1609.
29. F. Wu, S. Sharma, S. Roy, P. Halstenberg, L. C. Gallington, S. M. Mahurin, S. Dai, V. S. Bryantsev, A. S. Ivanov and C. J. Margulis, J. Phys. Chem. B, 2020, 124, 2892–2899.
30. J. C. Wasse and P. S. Salmon, J. Phys.: Condens. Matter, 1999, 11, 1381–1396.
31. Q. Li, D. Sprouster, G. Zheng, J. C. Neuefeind, A. D. Braatz, J. McFarlane, D. Olds, S. Lam, J. Li and B. Khaykovich, ACS Appl. Energy Mater., 2021, 4, 3044–3056.
32. S. Roy, F. Wu, H. Wang, A. S. Ivanov, S. Sharma, P. Halstenberg, S. K. Gill, A. M. Abeykoon, G. Kwon and M. Topsakal, Phys. Chem. Chem. Phys., 2020, 22, 22900–22917.
33. Y. Zhang and E. J. Maginn, J. Phys. Chem. Lett., 2015, 6, 700–705.
34. H. V. R. Annapureddy and L. X. Dang, J. Phys. Chem. B, 2014, 118, 8917–8927.
35. P. L. Geissler, C. Dellago and D. Chandler, J. Phys. Chem. B, 1999, 103, 3706–3710.
36. E. Pluharova, M. D. Baer, G. K. Schenter, P. Jungwirth and C. J. Mundy, J. Phys. Chem. B, 2016, 120, 1749–1758.
37. S. Roy and L. X. Dang, Chem. Phys. Lett., 2015, 628, 30–34.
38. S. Roy, M. D. Baer, C. J. Mundy and G. K. Schenter, J. Chem. Theory Comput., 2017, 13, 3470–3477.
39. S. Roy and V. S. Bryantsev, J. Phys. Chem. B, 2018, 122, 12067–12076.
40. S. Roy, M. Galib, G. K. Schenter and C. J. Mundy, Chem. Phys. Lett., 2018, 692, 407–415.
41. S. Roy, G. K. Schenter, J. A. Napoli, M. D. Baer, T. E. Markland and C. J. Mundy, J. Phys. Chem. B, 2020, 124, 5665–5675.
42. C. Cook, Jr, J. Inorg. Nucl. Chem., 1963, 25, 123–127.
43. L. B. Skinner, C. J. Benmore and J. B. Parise, Nucl. Instrum. Methods Phys. Res., Sect. A, 2012, 662, 61–70.
44. B. H. Toby and R. B. Von Dreele, J. Appl. Crystallogr., 2013, 46, 544–549.
45. X. Qiu, J. W. Thompson and S. J. Billinge, J. Appl. Crystallogr., 2004, 37, 678.
46. W. C. Phillips, R. Gakhar, G. P. Horne, B. Layne, K. Iwamatsu, A. Ramos-Ballesteros, M. R. Shaltry, J. A. LaVerne, S. M. Pimblott and J. F. Wishart, Rev. Sci. Instrum., 2020, 91, 083105.
47. S. Dai, G. Begun, J. Young and G. Mamantov, J. Raman Spectrosc., 1995, 26, 929–932.
48. C.-H. Huang and M. Brooker, Chem. Phys. Lett., 1976, 43, 180–182.
49. G. J. Janz, R. Tomkins, C. Allen, J. Downey, Jr, G. Garner, U. Krebs and S. K. Singer, J. Phys. Chem. Ref. Data, 1975, 4, 871–1178.
50. M. Salanne, L. J. Siqueira, A. P. Seitsonen, P. A. Madden and B. Kirchner, Faraday Discuss., 2012, 154, 171–188.
51. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
52. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
53. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
54. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
55. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
56. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
57. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703.
58. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
59. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
60. M. G. Tucker, D. A. Keen, M. T. Dove, A. L. Goodwin and Q. Hui, J. Phys.: Condens. Matter, 2007, 19, 335218.
61. Y. P. Zhang, M. Eremenko, V. Krayzman, M. G. Tucker and I. Levin, J. Appl. Crystallogr., 2020, 53, 1509–1518.

---

Footnote

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

---