## 2 Materials and Methods  
A co-melt of FLiNaK and CeF\({}_{3}\) was obtained by melting four fcc crystals: LiF (9300 ions), NaF (2300 ions), KF (8400 ions), and CeF\({}_{3}\) (750 ions) at 3000 K for 1 million time steps (\(\Delta t=0.1\) fs). A system of this composition corresponded to the FLiNaK eutectic: LiF--46.5 mol.%, NaF--11.5 mol.%, and KF--42.0 mol.%, in which CeF\({}_{3}\) (15 mol.%) is dissolved. Complete mixing of the melt was achieved after holding at the same temperature for another 2 million \(\Delta t\). All these calculations were performed in the NVT ensemble. At the next stage of the calculation in the NPT ensemble, molten salt from FLiNaK and CeF\({}_{3}\) was brought to the experimental density corresponding to a temperature of 800 K and kept at this temperature. The duration of this stage was 3 million \(\Delta t\). The Nose-Hoover thermostat and barostat were used for keeping the simulated temperature and pressure constant [11; 12]. The thermostat relaxation time was \(\tau_{T}=400\)\(\Delta t\), and the barostat relaxation time was \(\tau_{P}=1000\)\(\Delta t\). Exactly the same values of the parameters \(\tau_{T}\) and \(\tau_{P}\) were used in [13]. Subsequent calculations were carried out in the NVE ensemble by sequentially increasing the temperature by a certain value \(\Delta T\), keeping the system at a new temperature, and then calculating the physical properties. Each calculation at the new temperature took 2 million \(\Delta t\). In all cases, periodic boundary conditions acting in each of the three directions were used. Long-range interactions were taken into account on the basis of Particle-Particle Particle-Mesh (PPPM) summing [14]. The Poisson equations were solved using the fast Fourier transform. The interaction of singly charged ions that make up FLiNaK was described by the rigid-ion potential of Born-Mayer-Huggins-Tosi-Fumi (BMHFT) [15]  
\[\phi=Z_{i}Z_{j}\frac{e^{2}}{4\pi\epsilon_{0}\epsilon r}+B_{ij}\exp\left[- \frac{r-\sigma_{ij}}{\rho_{ij}}\right]-\frac{C_{ij}}{r^{6}}-\frac{D_{ij}}{r^{ 8}}, \tag{1}\]  
where \(i\) and \(j\) can represent positive or negative ion, \(Z_{i}\) and \(Z_{j}\) characterize ion charges; the last two terms in expression (1) reflect the dipole-dipole and dipole-quadrupole interaction; and \(\epsilon_{0}\) and \(\epsilon\) are the permittivity of the free space and relative permeability of the studied environment, respectively. The parameters for the potential for FLiNaK represented by expression (1) are given in [16]. The QUANTUM ESPRESSO program [17; 18] was used to calculate the parameters of the Ce-Ce and Ce-X interaction potential, where X is F, Li, Na, K. Projector-augmented wave GGA-PBE scalar relativistic pseudopotentials [19; 20] were used. Moreover, the kinetic energy cutoff for wave functions and charge density corresponded to 100 and 400 Ry, respectively. The interaction of the Ce ion with other ions was created by the Coulomb and non-Coulomb contributions. When determining the Coulomb contribution, integer values of electric charges were used: +3 for Ce ion; +1 for Li, Na, and K ions; and \(-1\) for F ions. The schematic diagrams of the complexes under study are shown in Figure 1. The complexes were placed in a simple cubic lattice with a large parameter a = 30 A, so that the interaction between periodic cells was practically excluded. In other words, non-periodic isolated systems were studied. The geometric template for determining the Ce-F interaction was represented by a planar CeF\({}_{3}\) molecule, in which the cerium atom is located at the center of a regular triangle. The Ce-F distance was varied while maintaining the symmetry of the system. To determine the Ce-Ce interaction, a 3D complex consisting of two Ce atoms and six F atoms was considered. All Ce-F distances in this system were assumed to be equal to the typical distance for melting salts (2.3 A) while the Ce-Ce distance was varied.