Temperature-Dependent Properties of Molten Li2BeF4 Salt Using Ab Initio Molecular Dynamics.

Molten lithium tetrafluoroberyllate (Li2BeF4) salt, also known as FLiBe, with a 2:1 mixture of LiF and BeF2 is being proposed as a coolant and solvent in advanced nuclear reactor designs, such as the molten salt reactor or the fluoride salt cooled high-temperature reactor. We present the results on the structure and properties of FLiBe over a wide range of temperatures, 0-2000 K, from high-throughput ab initio molecular dynamics simulation using a supercell model of 504 atoms. The variations in the local structures of solid and liquid FLiBe with temperature are discussed in terms of a pair distribution function, coordination number, and bond angle distribution. The temperature-dependent electronic structure and optical and mechanical properties of FLiBe are calculated. The optical and mechanical property results are reported for the first time. The results above and below the melting temperature (∼732 K) are compared with the experimental data and with data for crystalline FLiBe. The electronic structure and interatomic bonding results are discussed in correlation with the mechanical strength. A novel concept of total bond order density (TBOD), an important quantum mechanical parameter, is used to characterize the internal cohesion and strength in the simulated models. The results show a variation in the rate of change in properties in solid and liquid phases with anomalous behavior across the melting region. The observed trend is the decrease in mechanical strength, band gap, and TBOD in a nonlinear fashion as a function of temperature. The refractive index shows a surprising minimum at 850 K, among the tested temperatures, which lies above the melting point. These findings provide a new platform to understand the interplay between the temperature-dependent structures and properties of FLiBe salt.

Introduction

Molten salts (MSs) are a special class of ionic liquids with specific properties and a diverse range of applications.[1,2] In recent years, they have received tremendous interest along with the relatively rapid development of advanced nuclear energy technologies. The most common MSs are nitrates, for use at relatively low temperatures below 560 °C, and alkali fluorides, chlorides, and hydroxides and their mixed combinations, for use at higher temperatures. They have a variety of important applications such as in thermal energy storage,[3,4] advanced nuclear reactor technologies,[5−9] concentrated solar power and battery industry,[10−14]and so forth. Relatively recently, there has been a resurgence of interest in advanced nuclear technologies with the development of several high-temperature reactor designs using molten salt as a coolant or a fuel, which have the potential to change the landscape of carbon-free energy generation with much reduced cost.[15] A fundamental understanding of the structure and properties of MSs and their interaction with the surroundings and contents is necessary to accelerate their applications in the current and emerging advanced nuclear technologies.[16,17] However, there are still many aspects of scientific issues on MS that are not well-understood.

Next-generation nuclear reactors are drawing considerable attention worldwide. Among the leading candidate is the molten salt reactor (MSR) where molten fluoride salts are being considered both as a primary coolant and as a solvent for nuclear fuel. Lithium tetrafluoroberyllate (FLiBe) is one of the fluoride salts that have been extensively used in the MSR experiment(s) as a coolant and solvent for fissile material due to its advantageous thermophysical properties.[1,2,18−27] Properties such as small absorption cross section for thermal neutrons, high thermal efficiency, low spent fuel per unit energy, high degree of passive safety, high heat capacity, atmospheric pressure operation, chemical stability at high temperature, low melting point, and high boiling point make FLiBe a primary candidate for application in MSRs.[9,22,28−31] FLiBe has a melting and boiling point of ∼732 and ∼1703 K, respectively.[31] The liquid salt is a good solvent for nuclear fuel as it can dissolve high concentrations of Th + U fuel.[24] FLiBe can also be used as a tritium breeding material in fusion reactors because of its high chemical stability.[8,19,21,32] Thus, complete information regarding the temperature-dependent atomic structure and macroscopic properties is crucial for the further advancement of FLiBe in nuclear applications.

Although the experimental efforts have been accumulating[27,33−42] in parallel with simulation,[26,43−50] the lack of fundamental understanding at the atomistic level has been an issue in the required application of FLiBe for the MSR. Recently available specialized facilities and instrumentation have enhanced the state-of-the-art experimental research on molten salts.[31,51−54] However, an understanding of the connectivity between atomistic properties with experimental characterization results is crucially important especially given that FLiBe is toxic and corrosive at high temperatures which precludes cost-effective experimental research over a wide range of temperatures.[55] The profound effect of trace impurities in molten salt on its corrosivity and contamination from corrosion products further makes experimental work extremely challenging particularly at higher temperatures.[16] Hence, the implications on the corrosion effect of FLiBe such as the formation of oxides by the reaction with water and with reactor container metallic alloys must be addressed.[31,51,56] To limit the extent of the corrosion, molten salt redox potential control has been proposed with either the addition of a metallic buffer, a soluble redox couple buffer, or overhead gas control. In addition to the redox potential, fluoroacidity, which describes the coordination of the solvent, is also an important parameter of the salt. FLiBe is a fluoroacid salt since BeF2 can accept electrons from fluoride ions and form complexes with fluorides. As an example, the atomic interaction and change in properties from the amorphous solid to liquid state need to be fully elucidated to reveal the implications on the reaction 2LiF + BeF2 → 2Li+ + BeF42– at high temperature and estimation of redox potential[52,57] since control of redox potential is essential to mitigate the corrosion by FLiBe in MSR. The formation of these monomeric fluoro complexes helps in minimizing the fluorine potential of the melt and reduces corrosion. To address these challenges, a comprehensive understanding of the structure and properties of FLiBe over a wide range of temperatures is necessary.

Experimental studies on FLiBe salt are still limited due to the high costs for such work associated with its toxicity and corrosivity. Accurate and effective computer simulation is a viable alternative to complement the experiments and provide new insights to guide the experiments. There are relatively few such studies on FLiBe and they are limited to a structure at a particular or a narrow range of temperatures. Even where classical and first-principles molecular dynamics (MD) have been used to some extent to study crystalline,[58] solid, and liquid FLiBe[26,43,44,47−50,59] using different methods and potentials, they are mostly focused on the structure, dynamical features, thermodynamic and vibrational properties, thermal neutron scattering cross section data, and so forth, whereas information on electronic, optical, and mechanical properties of FLiBe over a wide range of temperatures is still not covered.

This paper expounds the temperature-dependent structural, physical, and electronic properties of FLiBe salt at an atomistic level based on ab initio molecular dynamics (AIMD). Our motivation is to use accurate large-scale modeling to unravel some of the missing links discussed above and to provide future directions in the study of this salt. We start with a detailed investigation of the crystalline FLiBe. Then, the structures and properties of FLiBe as a function of temperature are systematically studied covering both the amorphous and the liquid phases at 10 different temperatures (0, 300, 450, 700, 730, 750, 850, 1000, 1500, and 2000 K). With the melting point being ∼732 K, these temperatures cover the solid and liquid state. The structure at 2000 K is above the boiling point of FLiBe (∼1703 K) and is treated as a superheated liquid. In this paper, we first provide a description of the computational methods and approaches and then results with detailed descriptions starting with the crystalline model followed by amorphous and liquid phases. These results are further discussed in connection with experimental findings.

Methods and Models

The accuracy of results in any atomistic simulation depends mainly on the accuracy of interactions between atoms in the system. The use of empirical force field potential in many classical MD simulations lacks sufficient accuracy compared to AIMD. It cannot, for example, predict the relevant electronic structure properties when the structure changes with temperature. Each new condition requires the fitting of a set of potential parameters which is limited in some cases due to inevitable indeterminacy. Moreover, such a potential cannot reproduce the temperature-dependent physical properties. In this context, the use of density functional theory (DFT)-based AIMD is the more accurate route to simulate the electronic structure and properties. Hence, AIMD is necessary to capture the critical details for the crystalline, amorphous, and liquid phases of FLiBe at different temperatures. Although the AIMD provides highly reliable results, it has a limitation in the case of a large system with several thousands of atoms and for long-time calculation. For such a simulation, an MD study using machine learning-based interatomic potential is shown to be more efficient with comparable accuracy to DFT.[60−63]

In this study, we use two DFT-based methods: Vienna Ab initio Simulation Package (VASP)[64,65] and orthogonalized linear combination of atomic orbitals (OLCAO).[66] VASP is used for AIMD simulation, geometry optimization, and calculation of mechanical properties. We use the projector augmented wave (PAW) method[67,68] and describe the exchange–correlation energy by Perdew–Burke–Ernzerhoff (PBE) generalized gradient approximation (GGA).[69] The Brillouin zone is sampled at Γ-point only for supercell models and k-point mesh (5 × 5 × 7) for the crystalline system. A high cutoff energy of 500 eV is set for the plane-wave basis expansion, and the electronic convergence criterion is chosen to be 10–5 eV. For the ionic relaxation, forces are converged to 10–3 eV/Å at zero pressure. The VASP relaxed structures are used as input to the OLCAO package for electronic structure and optical property evaluation. Details about the OLCAO method are presented in Supporting Information.

The initial crystal structure of FLiBe is taken from X-ray diffraction data.[70] It has a rhombohedral symmetry hexagonal lattice defined by cell lengths “a” and “c”. The amorphous and liquid FLiBe models are simulated using the melt-quenching technique using AIMD similar to our previous studies.[71−73] All simulations are carried out under NVT ensemble where the number of atoms (N), volume (V), and temperature (T) are held constant. The heat bath temperature is controlled using a Nose thermostat.[74,75] A time step of 1 fs is chosen for the ionic motion integration. Based on the unit cell crystal structure with 126 atoms (Li = 36, Be = 18, and F = 72), we construct a cubic supercell containing 504 atoms (Li = 144, Be = 72, and F = 288), with periodic boundary conditions. The initial cubic supercell is heated to 2000 K within just 2 ps. The melt is further heated for 10 ps at this elevated temperature to fully eliminate the memory effect of the initial configuration. After this, the temperature of the system is lowered to 1500, 1000, 850, 750, 730, 700, 450, and finally to 300 K. At each temperature, the system is well-equilibrated for 10 ps to ensure it overcomes the diffusive stage and losses memory of atomic position and velocity history configuration from the previous structure. The velocity autocorrelation function (VACF) (see Supporting Information, Figure S1) ensures the loss of initial velocity from the previous configuration in the simulated model. It provides information on the dynamic motion of atoms with time. It shows that in the liquid model, the VACF dies out fast implying that the ions can leave the cage made up of surrounding ions more quickly. In a solid model, the VACF has more features and they die out much more slowly compared to melt which is due to the presence of more BeF4– tetrahedrons and the ordered structure of solid FLiBe. Out of 10 ps equilibration time, the pressure and temperature fluctuations are minimized through 8 ps, and an averaged structure is taken from the last 2 ps run. A nominal cooling rate of 20 K/ps is used to anneal the structures, which is reasonable in an accurate AIMD simulation.[76] The equilibrium volume at each temperature is determined from the AIMD run of the system in the NPT ensemble at atmospheric pressure. The ensemble uses a Parrinello–Rahman method[77,78] with a Langevin thermostat[79] as implemented in VASP. In the Langevin thermostat, the temperature of the system is maintained through a modification of Newton’s equation of motion by introducing a friction coefficient and random force. The friction coefficients for all atomic and lattice degrees of freedom are taken to be 20 ps–1 which are system-dependent and can only be optimized with enough testing. In the NPT ensemble, pressure (P) is kept constant along with N and T. The final structure at 0 K is obtained from the structure at 300 K by optimizing the atomic position and cell volume to converge the force and energy to a local minimum state at 0 K. Table lists the cell parameters of the simulated model at different temperatures. The same naming of models used in Table will be applied to all figures, results, and discussion sections.

Table 1

Structural Parameters for the FLiBe Model at Different Temperaturesa

models	a (Å)	b (Å)	c (Å)	αo	βo	γo	density	volume
crystal-exp*	13.281	13.281	08.888	90.00	90.00	120.00	2.180	1357.700
crystal-sim	13.447	13.447	08.995	90.00	90.00	120.00	2.144	1408.729
T0	18.405	18.405	18.405	90.00	90.00	90.00	1.908	6234.582
T300	18.538	18.538	18.538	90.00	90.00	90.00	1.866	6370.928
T450	18.625	18.625	18.625	90.00	90.00	90.00	1.840	6460.838
T700	18.760	18.760	18.760	90.00	90.00	90.00	1.800	6602.824
T730	18.801	18.801	18.801	90.00	90.00	90.00	1.788, 1.97–2.08†	6645.732
T750	18.915	18.915	18.915	90.00	90.00	90.00	1.757, 1.95–2.05†	6767.367
T850	19.009	19.009	19.009	90.00	90.00	90.00	1.731, 1.90–2.02†	6868.752
T1000	19.214	19.214	19.214	90.00	90.00	90.00	1.675, 1.80–1.95†	7093.522
T1500	19.519	19.519	19.519	90.00	90.00	90.00	1.597	7437.105
T2000	19.597	19.597	19.597	90.00	90.00	90.00	1.571	7525.745

Volume in Å3, density in gm/cm3, Crystal-exp = experimental structure, and Crystal-sim = simulated structure. *Ref (70), †The experimental range of density is taken based on the values reported in refs (80)–[85].

Results

Crystalline FLiBe

Crystalline FLiBe has been well-studied in the past.[58,70,86−92] It has a rhombohedral symmetry of space group R3̅ (148) with a hexagonal setting. The initial cell parameters and atom position data, taken from X-ray diffraction work by Seiler,[70] are relaxed using VASP (see Section ). The lattice parameters, density, and volume of the original and final relaxed crystal are presented in Table . Compared to experimental values, our calculated cell volume is increased by 3.75% and hence the density is underestimated by 1.65%. The ball and stick representations of this crystal in two different orientations are shown in Figure . Each Li+ and Be2+ ions form tetrahedrons with F– ions. Each F– ion shares one Be and two Li as the nearest neighbors. These tetrahedra are interconnected along the c-axis in the hexagonal cell. The Be–F bond length (BL) in BeF42– varies from 1.568 to 1.578 Å and the averaged value is 1.573 Å. There exist two different types of LiF43– tetrahedra, so the Li–F bonds are distributed over a wider range from 1.874 to 1.913 Å, than Be–F bonds. These observed BLs are in close agreement with the experimental result, although slightly overestimated, which are in the range of 1.550–1.561 and 1.846–1.889 Å for Be–F and Li–F pairs, respectively.[70] The variation of BL between different types of bonds with the corresponding bond order (BO) values is presented in Figure to be discussed later. It is observed that the Be–F bonds are congregated with overlapped BO values, while the Li–F bonds are more scattered with slightly different BO values. Hence, the key difference between Be–F and Li–F bonds is that Be–F bonds are nearly identical with high BO values, whereas Li–F bonds spread in a certain range and have low BO values. The next nearest neighbor Li–F pairs are at 3.315 Å with almost zero BO.

Figure 1

Ball and stick structures of relaxed crystalline FLiBe in two different orientations (Li = green, Be = blue, and F = orange).

Figure 9

Bond order vs bond length plot for FLiBe at different temperatures. The inset shows the distribution of Be–F bonds up to the cutoff distance of 2.5 Å.

Even though numerous studies explore the geometry of the FLiBe crystal, the complete electronic structure results at the atomic level are rarely reported.[27,70,87] The calculated density of states (DOS) for this crystal is shown in Figure a. The zero of the energy is set at the top of the valence band (VB). The VB DOS has two parts: the upper and lower. The lower part lies below −19 eV and spreads in the 1 eV energy range. The F’s 2s orbital is the main contributor for the DOS in this region, while a small contribution comes from Be and Li’s 2s and 2p orbitals. The upper part which spreads over the 3.6 eV range can be divided into three segments A, B, and C separated by a small energy gap. B and C have one sharp peak, while A has one major peak followed by two minor peaks. The atom-resolved partial DOS (PDOS) shows the orbitals of Be contribute to the segments B and C only, while those of Li and F contribute to all three segments. Therefore, a distinctive feature of DOS is that Be does not have any contribution for the upper VB. The DOS spectrum in the conduction band (CB) region is rather rough with many pronounced peaks of varying intensities which is the consequence of different types of Be–F and Li–F bond pairs. The PDOS shows that F mainly contributes to the VB, while Li and Be contribute to the CB. An important feature is that the DOS for crystalline FLiBe is completely different from the amorphous model at 0 K (T0) (see Figure to be discussed later). The calculated band structure for the FLiBe crystal is shown in Figure b. This material is an insulator with a direct large band gap of 7.172 eV at the Γ zone. This value could be slightly smaller than the real band gap value as the local density approximation adopted in this calculation has been well-known to underestimate the band gap in an insulator.[93−95] A similar band structure result is reported in a DFT calculation applying GGA[87] with a band gap value of 7.599 eV, a little larger than our value. The top of bulk VB is quite flat implying a large hole effective mass. The band gap of crystalline FLiBe is 1.02 eV larger than that of the amorphous model T0.

Figure 2

Calculated (a) total and atom-resolved partial density of states and (b) band structure for crystalline FLiBe.

Figure 7

Calculated total and atom-resolved partial density of states for the simulated FLiBe model at 0, 730, and 1500 K.

The calculated total bond order density (TBOD), pair-resolved partial BOD (PBOD), and partial charge (PC) for each atom in the crystalline model are listed in Table . The TBOD for crystalline FLiBe is higher than that of T0 and other models shown in the same table. The higher TBOD value for the crystal indicates that it is internally stronger than the corresponding amorphous and liquid models. The PBOD is higher for the Be–F pair which implies that this bond is stronger than Li–F. The PC is positive for Li and Be and negative for F, so both Li and Be lose charge to F atoms as expected.

Table 2

Calculated TBOD, PBOD, Averaged Partial Charge of Ions, Band Gap (E), Refractive Index (n), and Plasma Frequency (ωp) in the Simulated FLiBe Salt

		PBOD (e/Å3)			partial charge (in e)
models	TBOD (e/Å3)	Li–F	Be–F	Li–Li	Li	Be	F	E (eV)	n	ωp (eV)
crystal	0.01959	0.00978	0.00981	0.00000	0.568	0.866	–0.501	7.172	1.264	19.71
T0	0.01718	0.00841	0.00877	0.00000	0.569	0.880	–0.504	6.152	1.211	17.34
T300	0.01662	0.00809	0.00853	0.00000	0.570	0.883	–0.506	6.047	1.208	17.25
T450	0.01617	0.00782	0.00834	0.00001	0.575	0.885	–0.509	5.875	1.207	17.38
T700	0.01559	0.00748	0.00810	0.00001	0.578	0.890	–0.510	5.788	1.206	17.23
T730	0.01552	0.00746	0.00805	0.00001	0.579	0.891	–0.511	5.636	1.203	17.02
T750	0.01516	0.00727	0.00789	0.00001	0.580	0.892	–0.513	5.612	1.201	16.42
T850	0.01489	0.00708	0.00780	0.00001	0.585	0.894	–0.514	5.598	1.199	16.62
T1000	0.01408	0.00669	0.00738	0.00001	0.587	0.898	–0.518	5.033	1.203	16.37
T1500	0.01315	0.00625	0.00688	0.00001	0.591	0.905	–0.522	4.797	1.205	16.21
T2000	0.01284	0.00619	0.00663	0.00001	0.585	0.922	–0.523	3.982	1.210	15.42

The calculated optical properties for the FLiBe crystal in the form of frequency-dependent complex dielectric functions are shown in Figure a. The imaginary dielectric function ε2 has many absorption peaks with different intensities with the onset energy at 7.17 eV. The major absorption peaks denoted by P1–P6 occur at 10.48, 12.9, 14.52, 16.34, 18.50, and 21.77 eV, respectively. These peaks in ε2 are consistent with the sharp peaks in CB and VB DOS spectra. The real part ε1 is obtained from imaginary part ε2 using Kramer-Kroning transformation[96] and it shows features similar to that of ε2. The refractive index is calculated from the square root of ε1(0) to be 1.264 which is higher than that for T0 and other models (see Table ). The energy loss function (ELF) for this crystal calculated from the inverse of the complex dielectric function is shown in Figure b. The ELF spectrum has one sharp peak followed by minor peaks on both sides. This spectrum is asymmetric with a long plateau on the right side above 30 eV. The plasmon peak identified by the peak position of the spectrum with frequency ωp lies at 19.71 eV which defines the frequency of collective excitation of valence electrons in the system. This value in the crystal is higher than that for the amorphous and liquid model listed in Table .

Figure 3

Calculated (a) optical dielectric functions and (b) energy loss function for the FLiBe crystal.

The mechanical property of a crystal depends on its symmetry. In a hexagonal setting, the elastic properties are defined by five independent elastic constants C11, C12, C13, C33, and C44, whose values in the FLiBe crystal are 75.581, 38.512, 32.115, 98.142, and 18.529 GPa, respectively. From the elastic coefficients C and related compliance ratio S, the mechanical parameters can be calculated (see Supporting Information). The Young’s modulus (E), bulk modulus (K), and shear modulus (G) derived from these elastic constants are 54.091, 50.359, and 20.474 GPa, respectively. The elastic moduli for the crystal are larger than those of the T0 model (see Table ) which implies that the internal forces and strength are higher in a crystal than in the amorphous model. The calculated Poisson’s ratio (η) and Pugh’s modulus ratio (G/K) for this crystal are 0.321 and 0.406, respectively. The G/K is an important criterion defined by Pugh to separate the brittleness (G/K > 0.57) and ductility (G/K < 0.57) of a material.[97] It follows a trend just opposite to η. The higher Poisson’s ratio (or G/K < 0.57) in the crystalline model indicates that it is more ductile or less brittle as compared to amorphous and liquid models.

Table 4

Calculated Elastic Constants and Elastic Moduli (in Gpa) for Amorphous and Liquid FLiBe at Different Temperaturesa

models	C11	C12	C44	C44′	E	K	G	η	G/K
T0	52.063	21.041	17.803	15.511	42.955	31.382	16.886	0.272	0.538
T300	49.757	19.257	17.303	15.250	41.666	29.423	16.482	0.264	0.560
T450	38.561	12.956	12.477	12.803	31.628	21.511	12.601	0.255	0.586
T700	35.773	10.254	12.413	12.760	30.789	18.760	12.552	0.226	0.669
T730	40.152	11.274	14.326	14.439	35.075	20.900	14.371	0.220	0.688
T750	30.459	10.610	13.409	09.924	29.247	17.226	12.015	0.217	0.698
T850	29.806	09.079	12.044	10.364	27.577	15.988	11.372	0.213	0.711
T1000	33.827	11.195	11.483	11.316	28.468	18.739	11.416	0.246	0.609
T1500	37.800	12.027	12.303	12.887	31.272	20.618	12.537	0.247	0.608
T2000	28.483	08.697	12.173	09.893	27.126	15.292	11.261	0.204	0.736

C44′= (C11 – C12)/2 (see Supporting Information).

Amorphous and Liquid FLiBe

We discuss our results dividing them mainly into two regions according to different phases of the system: liquid structures above the melting point (>732 K) and amorphous structures below the melting point (<732 K) including the structure T730 close to the melting point. The final equilibrated structures, in the ball and stick form, at each temperature are shown in Figure S2.

Structural Information

The structure of noncrystalline FLiBe becomes more disordered with the increase in temperature. This is illustrated in Figure for models T0, T730, and T1500, respectively. The different cell sizes used here are proportional to their volume at different temperatures. Be2+ and F– ions form a tetrahedral linkage of BeF42–, while the Li+ ions interact with these tetrahedra and balance the charge. With the increase in temperature, a number of these tetrahedrons decreases due to the elongation of Be–F bonds and they become more scattered than that in T = 0 K. In the section for bonding analysis, we will provide more detailed information on these structures rather than the discussion made here based on merely a geometric sketch. The volume increases and the density decreases when we move from T0 to the T2000 model as anticipated. The reported experimental density for T730, T750, T850, and T1000 models is in the range of 1.97–2.08, 1.95–2.05, 1.90–2.02, and 1.8–1.95 gm/cm3, respectively.[80−85] Our calculated density for T730, T750, T850, and T1000 (see Table ) is underestimated by 9.2, 9.9, 8.9, and 6.9%, respectively, of the reported values as the application of GGA overestimates the cell volume[98] in addition to the rapid cooling rate normally associated with the AIMD simulations.

Figure 4

Ball and stick and polyhedral structures of final relaxed FLiBe at three different temperatures 0, 730, and 1500 K.

Pair Distribution Function

The pair distribution function (PDF) defined as , where N is the total number of atoms and ρ is the number density, gives the probability of finding a particle at a distance r from another particle. The calculated PDFs for the simulated model at different temperatures are shown in Figure . It depicts the local structural information and distribution of ions in the presence of neighboring ions. The inset on the top panel of the figure shows a magnified portion of the first peak. The Be–F pair has a very sharp symmetric peak followed by a deep minimum implying stable and compact BeF42– clusters and a relatively well-defined first coordinated cell. The intensity of the peak decreases from T0 to T700 with a sharp decrease on moving from T0 to T300. After 700 K, the intensity slightly increases and remains almost constant throughout 850–2000 K. The peak position shifts slightly right with the increase in temperature. For the Li–F pair, the peak is broader as compared to Be–F and the distribution is somewhat asymmetric with its minimum being above zero implying a possible exchange between F– and Li+ ions. The intensity of the Li–F peak decreases with an increase in temperature. The peak position shifts slightly to the right from T0 to T450, but after that, the BL decreases slightly (by 0.01 Å), and hence, the peak position also shifts left and then remains almost constant in the 700–2000 K range. For the F–F pair, the intensity decreases, peak width increases, and the peak position shifts to the right with the increase in temperature. The average BL for Be–F, Li–F, and F–F pairs in T0 is 1.567, 1.887, and 2.572 Å, respectively. A similar result (1.564 Å) is reported for the Be–F pair in a computational study of FLiBe at 32 K.[44] From X-ray diffraction data, Vaslow and Narten[38] found that the Be–F, Li–F, and F–F pairs have relatively similar BLs 1.58, 1.85, and 2.563–3.020 Å, respectively, at 828 and 1023 K. Our calculated BLs for Be–F, Li–F, and F–F pairs at 850 and 1000 K are 1.575, 1.893, and 2.588 Å, respectively, which are consistent with the experimental finding of Vaslow and Narten.[38] In a computational study of FLiBe at 973 K, Nam et al.(47) found the Be–F, Li–F, and F–F BLs to be at 1.55, 1.88, and 2.58 Å, respectively. The Be–F, Li–F, and F–F BLs reported by Salanne et al. in a theoretical study of FLiBe at 873 K are 1.58, 1.81, and 2.61 Å, respectively.[48] The finding of a decrease in intensity and increase in width of distribution for both Be–F and Li–F pairs indicate that Be and Li solvation shells become loose with the increase in temperature. The averaged BL calculated based on the first peak of the PDF follows the order of Be–F < Li–F < F–F. The PDFs for cations pairs Li–Li, Li–Be, and Be–Be are rather spiky showing the different possible distribution of these ions. The change in peak positions in these pairs is irregular.

Figure 5

Pair distribution function for the simulated model at different temperatures. The insets in figures (a–c) show the magnified image of the first peak.

Coordination Number and Bond Angle Analysis

The first shell coordination number (CN) of Li and Be ions with the surrounding F ions is calculated by the numerical integration of PDF defined as where rc is the cutoff distance estimated from the first minimum of the corresponding PDF. The calculated CNs for Li and Be with two different cutoff distances are presented in Table . The CN value is large with a longer cutoff distance as expected. In the T0 model, Be maintains tetrahedral coordination, while the CN for Li is slightly higher than 4. Even if both Li and Be ions retain close coordination with the F atom, their local structure reflected by PDF are totally different as explained above. The average CN of Be remains constant up to 450 K and then decreases when the temperature is increased from 0 to 2000 K as the interatomic separation with F increases with the increase in temperature. In the case of Li, the CN decreases for shorter cutoff distance and increases for longer cutoff distance with the increase in temperature. However, in both cases, the decrease or increase is not regular. In an X-ray diffraction study of FLiBe, Vaslow and Narten determined the tetrahedral coordination for both Be and Li ions at 828 and 1023 K.[38] In a theoretical study of FLiBe, Nam et al.(47) calculated the CN for Be and Li to be 4 and 4.7, respectively (at 973 K), while Salanne et al.(48) calculated a CN of 4 (at 873 K) for both Be and Li ions.

Table 3

Coordination Number Distribution for Li and Be Atoms in the Simulated FLiBe Model with Two Different Cutoff Distances (CDs)

		coordination number (CN)
atoms	CD (Å)	crystal	T0	T300	T450	T700	T730	T750	T850	T1000	T1500	T2000
Li	2.80	4.00	4.24	4.31	4.24	4.26	4.33	4.19	4.30	4.10	4.08	4.14
	3.40	4.00	4.46	4.67	4.66	4.85	4.91	4.86	4.87	4.96	4.90	5.12
Be	1.80	4.00	4.00	4.00	4.00	3.98	3.97	3.96	3.94	3.79	3.60	3.32
	2.40	4.00	4.00	4.00	4.00	3.99	3.99	3.99	3.99	3.90	3.75	3.58

The distributions of bond angles F–Be–F and F–Li–F in the simulated model are shown in Figure . The F–Be–F angle for T0 is symmetric with a sharp peak centered at ∼109.5°, a perfect tetrahedral angle. The peak becomes broader, and the intensity decreases on moving from T0 to T2000. The distribution becomes wider and asymmetric and the peak position shifts to the right, on the higher angle side, with an increase in temperature. This shows systematic changes in bond angles in BeF42– tetrahedra toward lower CN consistent with increasing PC (see Section ) with the increase in temperature. The F–Li–F angle spreads over a wider range with a broader peak than the F–Be–F angle. The intensity decreases, the peak becomes broader, and the distribution becomes wider with the increase in temperature. Compared to F–Be–F, the shift of the peak position is small for F–Li–F, but the intensity drops significantly. The difference in bond angle distribution in F–Be–F and F–Li–F is due to different charges, CNs, and arrangements of Li and Be ions.

Figure 6

Distribution of (a) F–Be–F and (b) F–Li–F bond angles in FLiBe at different temperatures.

Electronic Structure

The electronic structure is the fundamental property of a material that assists to explore many physical properties. In particular, the evaluation of interatomic bonding facilitates the understanding of mechanical properties in a substance and its implication on the corrosion of materials. A thorough understanding of the electronic structure is only possible through exhaustive atomic-level quantum mechanical calculation rather than simply inferring from structural parameters. It is surprising to note that, like crystalline FLiBe, the published literature still lacks a complete result regarding the electronic structure in amorphous and liquid FLiBe. Our study fills this gap through rigorous DFT calculation using OLCAO.

The calculated total DOS data for FLiBe at different temperatures are shown in Figure S3. Overall, the DOS spectra look similar in the studied model at different temperatures. Like the crystalline model, the VB DOS is divided into two parts: lower and upper. The lower part DOS lies below −18 eV and it has one major peak. The upper part VB DOS spreads within the 0–7.5 eV energy range and it has one major and one minor peak. The minor peak becomes less prominent and starts to merge with the body of spectra from T1500. Compared to the crystalline model, the VB DOS has a fewer sharp peak, and it spreads over a wider energy range. The CB DOS is flatter and almost similar in all models. One of the striking features observed in the DOS spectra is their intensity decreases with the increase in temperature and the peak of spectra in the VB region becomes less sharp. We provide details on DOS spectra and their origination from constituent atoms by plotting the atom-resolved PDOS, shown in Figure , at 0, 730, and 1500 K to represent the three different temperature stages of the FLiBe. The PDOS shows that the deep lower VB DOS spectrum mainly comes from the F atom. It also shows that the upper VB DOS spectrum is mostly dominated by orbitals of F with an admixture of Li and Be orbitals. The Li and Be atoms mainly contribute to the CB DOS. When we move from T0 to the T2000 model, the band gap decreases continuously, which is reflected in the CB region of Figure S3. The calculated band gap values are presented in Table and their variation with temperature is shown in Figure . It is seen that the decrease in band gap is not linear, and it has a different slope in the solid and liquid regions. The vertical dashed line represents the melting point of the FLiBe. Most importantly above 850 K, where the structure is in the liquid phase, the band gap decreases sharply implying that in the melt, more ions are free and conductive.

Figure 8

Variation of TBOD and band gap as a function of temperature in FLiBe. The hollow symbols are data for the T2000 model.

Interatomic Bonding and Charge Transfer

The interatomic bonding is the most important part of the electronic structure calculation which helps to understand the internal cohesion and strength of the material. The calculated values of TBOD and PBOD are listed in Table . Like crystalline FLiBe, the Be–F pairs have a higher BOD value than Li–F which implies that Be–F bonds are stronger than Li–F. The variation of TBOD with temperature in the FLiBe system is displayed in Figure . The TBOD decreases with an increase in temperature but in a nonlinear fashion with a change in slope around the melting point. The most interesting result found is that the TBOD decreases sharply after 850 K, implying that liquid FLiBe is internally weaker than the solid phase.

We dig into the electronic structure results further by calculating the BO values for all possible bonded pairs and provide details of bonding analysis. The use of the BO value derived accurately from electronic structure calculation based on a real bond provides more information to characterize the interatomic strength rather than just using a geometrical distribution of an atom pair in a specified distance as in the PDF. The characteristics of bonding topology are vividly displayed in Figure in the form of the BO versus BL plot. It is observed that the crystalline, amorphous, and liquid FLiBe have a broad range of bonding patterns varying from very strong to weak bonds. The details of bonding in the crystalline model are discussed earlier in the paper, so in this section, we focus our discussion on the remaining models. In all models, the Be–F bonds are much stronger than Li–F bonds with a high BO value due to the presence of some partial covalent character. Collins et al. in a deformation density study of FLiBe has also shown that Be–F has high bonding density than Li–F, implying the partial covalent character of Be–F bonds.[92] In T0 and T300 models, the Be–F bonds are clustered in a narrow BL region, but as the temperature increases, they start to spread up from the T450 model. The figure in the inset displays the distribution of Be–F bonds in the range of 1.30–2.50 Å BL. It is seen that with an increase in temperature, more elongated Be–F bonds exist. The scattered BO values with varying BL at higher temperatures imply the different possibility of Be–F coordination than at lower temperatures where CN is 4. Also, for T0 and T300 models, the BO values for Be–F and Li–F pairs are separated by a narrow gap, but from the T450 model, they start to mix up and overlap each other, which increase as a function of temperature. Both Be–F and Li–F pairs spread in a wide region and become more random with the increase in temperature; however, the expansion for Be–F bonds is more systematic than Li–F. Near or beyond 3 Å, some Li–Li pairs, with a nearly zero BO value, exist which indicates that Li–Li covalent bonding is strictly prohibited in both amorphous and liquid FLiBe systems.

The calculated PC for each atom, presented in Figure S4, provides further insights into the mechanism of charge transfer between the ions in the studied system. The distribution becomes more scattered with the increase in temperature. The averaged PC values calculated for each atom in all models are plotted in Figure S5 and listed in Table . Like the crystalline model, the PC is positive for Li and Be ions and negative for F ions. Hence, the net charge transfer in the system is from Li and Be to F ions. With the increase in temperature from 0 to 2000 K, the PC increases for Be and decreases for F. The PC for Li increases from T0 to T1500 model but slightly decreases for T2000. In most of the MD calculations, the predefined same PC values are used for the simulation at different temperatures. The variation of PC with temperature in our result shows that only quantum calculation can account for such an effect precisely in the atomic interaction.

Optical Properties

The calculated real and imaginary parts of dielectric functions for the simulated models are shown in Figure S6. Compared to the crystalline model, the absorption curve ε2 in the amorphous and liquid models has entirely different features. The ε2 curve for the crystal is spiky with many peaks, while amorphous and melted models have one prominent peak, and the curve is flatter. However, for both cases, the peaks in ε2 are consistent with the corresponding peaks in DOS spectra in their CB region. An important observation is that when we move from T0 to the T2000 model, even if the feature of ε2 is almost similar, the absorption edge shifts at the lower energy end. The real part ε1 follows the same trend as ε2. The calculated refractive indices (n) for all models are listed in Table . The variation of n with temperature is shown in Figure . The value of n decreases first, reaches a minimum at 850 K, and then increases with temperature. Although the n value decreases from T0 to the T850 model, interestingly, the slope in the solid region below 700 K and across the melting point, in the region between 700 and 850 K, is different. The refractive index increases sharply from the T1500 to T2000 model. The calculated ELFs at different temperatures are presented in Figure S7. The ELF spectrum has one major peak followed by a minor peak with a long plateau on the right side. The calculated plasma frequencies (ωp) are listed in Table . It is observed that the ωp and intensity of the spectrum decrease with the increase in temperature from 0 to 2000 K.

Figure 10

Variation of the refractive index in the simulated model at different temperatures. The hollow square is the refractive index for the T2000 model.

Mechanical Properties

The understanding of temperature-dependent mechanical strength in FLiBe salt is important in the viewpoint of its application at high temperatures. In this study, we present the results on elastic and mechanical properties of FLiBe obtained using ab initio calculations. A set of linear elastic coefficients C and the mechanical properties derived from them at various temperatures are listed in Table . In the cubic cell, the independent elastic constants can be defined by three components C11, C12, and C44 which are listed in the table. In general, the amorphous and liquid structures are isotropic; however, due to the small cell size effect, some degree of anisotropy is unavoidable as indicated by slightly different values of C44 and C44′. The values of E, K, G, η, and G/K as a function of temperatures in the simulated model are shown in Figure . Overall, the values of E, K, and G decrease on moving from the T0 to T2000 model but irregularly. Noticeably, a different trend in decrease is observed in the solid and liquid region. From T300 to T700, the elastic moduli decrease sharply but from T730 to T850, the rate in decrease is slow. The elastic moduli increase on moving from T850 to T1500. An irregular change in elastic moduli values in the series on moving through T700 → T730 → T750 shows that this salt exhibits some unusual mechanical behavior across the melting point. The values of E, K, and G decrease sharply for the superheated liquid at 2000 K when the temperature is increased beyond 1500 K. The value of η decreases from T0 to T850 and increases from T850 to T1500 when the temperature is increased from 0 to 1500 K. η decreases sharply from T1500 to T2000. The G/K follows exactly the opposite nature of η.

Figure 11

Calculated (a) elastic moduli and (b) Poisson’s and Pugh’s modulus ratio at different temperatures for FLiBe salt. The hollow symbols disconnected from the curves are data for the T2000 model.

Discussion

Finally, we discuss our results presented above and highlight some exciting new findings in this study of FLiBe salt. Our results obtained from AIMD simulation provide a wealth of information on the electronic structure of crystalline, amorphous, and liquid FLiBe in a wide range of temperatures and reveal their interconnectivity with physical properties. This will facilitate the understanding of the temperature-dependent structure and properties of FLiBe which is necessary in the application of this salt to energy systems. Such data providing the information at the atomic level cannot be obtained directly by experimental means.

The crystalline FLiBe has substantially distinct structural and physical characteristics compared to amorphous and liquid models. This crystal is an insulator with a large direct band gap than the amorphous and liquid model. Even if Be and Li ions in the crystal and amorphous model at 0 K have the same (for Be) or close (for Li) CN value, their bonding pattern is entirely different. In the crystal, Be–F and Li–F bonds are concentrated in a limited range with almost overlapped BO values, but in the amorphous model, they are spread over a wide range with varying BO values. Another different result in the crystalline FLiBe is the existence of the nearest neighbor of Li–F pairs which do not exist in the amorphous and liquid model, rather these latter phases have Li–Li pairs. Additionally, the higher TBOD and elastic moduli in the crystal imply that it is internally stronger than amorphous and liquid FLiBe.

The amorphous and liquid phases of FLiBe also have distinct local structures. The PDF and bond angle distribution show that the melted structures are more random than the amorphous structures. Be mainly forms the BeF42– tetrahedra at the lower temperature; however, with the increase in temperature, the abundance of perfect tetrahedra decreases as indicated by a continuous decrease in CN of Be. The decrease in CN of Be is consistent with the increase in PC of Be with the increase in temperature. The CN values presented with two different cutoff distances show that they are higher with longer cutoff distances. It is consistent with the inclusion of more bonded pairs with the increase in BL, as shown in Figure . The CN of Be decreases for both cutoff distances, while the CN of Li decreases for shorter cutoff distance and increases for longer cutoff distance with the increase in temperature.

In MD calculations, the bonding is interpreted based on the PDF which is a purely geometric parameter. Our results provide a more comprehensive picture of the chemical bonding using the BO value derived from the quantum calculation which more accurately represents the internal strength of the studied material. The BO–BL plot in Figure provides a rich display of bonding distributions within these complex structures and allows us to gain a more complete picture of the internal distribution of the bond strengths. The Be–F bonds spread in a narrow range compared to Li–F and they have a higher BOD value, as listed in Table , implying that Be–F bonds have an ionic character as well as a partial covalent character, while Li–F bonds are mostly ionic. Calculation of PBOD reveals that Be–F bonds have a higher contribution to the overall bonding of the system than Li–F. For T0 and T300 models, Be–F and Li–F bonds have distinguishable BO values, but from T450, they start to mix up and overlap. This implies the elongation of BL in some Be–F pairs, as shown in the inset of Figure , and more randomness in the distribution of Be–F bonds with the increase in temperature. The shortening of BL in some Li–F pairs can also be seen in the same figure. However, the average BL follows the order of Be–F < Li–F, and the order of BO is just the opposite.

As emphasized above, the internal strength of the FLiBe salt is characterized by introducing the TBOD and its interconnectivity with mechanical properties is revealed. The highest TBOD for T0 indicates that this model is internally strong which is reflected in the calculated mechanical properties that T0 has high E, K, and G values. The TBOD decreases from T0 to T2000 and similarly, overall, a decrease in elastic moduli is observed on moving from T0 to T2000. However, since the elastic moduli do not depend merely on the TBOD, they do not follow the same trend as TBOD. The PC discussed in our study is another informative quantity obtained from ab initio calculation which is missing in many previous MD studies. We show that the PC of atoms in FLiBe strongly depends on the temperature. Our result shows that the PC increases for Li and Be, while it decreases for the F atom with the increase in temperature. Most of the MD simulations use the same PC value to study the temperature-dependent properties. In this context, our calculation provides more accurate results with detailed information of charge on each atom as compared to other studies.

The calculated refractive index for amorphous and liquid FLiBe is smaller than that of the crystalline model. In a rather old study,[89] the refractive index of the FLiBe crystal is reported to be ≪1.331, and there have been no recent and detailed studies in this area. In our study, the refractive index decreases first from 0 to 850 K and then increases with the increase in temperature. An interesting result found is that the refractive index is minimum in the liquid phase at 850 K. It is noteworthy that there are no existing theoretical and experimental studies on the optical properties of amorphous and liquid FLiBe. Our results can be a pathway for the exploration of the refractive index by experimental or other means. Our results show an abrupt change in the properties of FLiBe such as a sharp decrease in TBOD and band gap and an increase in elastic moduli, Poisson’s ratio, and refractive index at 850 K which indicates that there will be some structural changes in FLiBe melt at this temperature. This observation is also supported by a sharp decrease in the first shell CN for both Li and Be ions above 850 K. Overall, our study reveals a clear demarcation that separates the bonding and optical and mechanical properties of crystalline, amorphous, and liquid FLiBe. These findings may pave a new way to diagnose and distinguish these structures by experimental means.

Conclusions

This AIMD simulation performed over a wide range of temperatures enhances the knowledge of the structure and properties of FLiBe salt. One of our findings is that the crystalline and amorphous FLiBe at the same temperature (0 K) have different local structures and bonding. Crystalline FLiBe has a comparatively higher density, band gap, TBOD, refractive index, and elastic moduli than that of the corresponding amorphous and liquid models. The FLiBe crystal appears to be more ductile than its amorphous and liquid counterparts. We show the effective use of the concept of BO to analyze the internal bonding in the FLiBe system. The association between Be and F ions is stronger than that for Li and F as characterized by a higher BO value for the Be–F pair in all models. Overall, the density, band gap, TBOD, and mechanical strength decrease with the increase in temperature from 0 to 2000 K. Another significant finding is that solid and liquid phases of FLiBe have a different rate of change in properties with some anomalous behavior across the melting point. The most important discovery is that the refractive index has a minimum value near 850 K where the FLiBe is in the liquid state. To the best of our knowledge, the optical and mechanical properties of FLiBe salt over the wide temperature range are reported for the first time and we anticipate that our results will motivate future experimental works in this direction.