**3. Results and Discussion**

In this study, the pair correlation function *g*(*r*) is generated by the VASPKIT code [29] in order to characterize the structure of the liquid metal system. The function *g*(*r*) is formalized as

$$\log(r) = \frac{1}{4\pi r^2} \frac{1}{N\rho} \sum\_{i=1}^{N} \sum\_{j=1 \atop i \neq i}^{N} < \delta \left(r - |r\_i - r\_j|\right) > \tag{1}$$

where *r* is the radial distance, *N* is the total number of atoms in the modeling system and *ρ* is the average density of the system. The purpose of normalizing the *g*(*r*) function by the density is to ensure that the radial distribution approaches unity for the long radial distance. Following Equation (1), the partial pair correlation function between two elements A and B can be written as

$$\log\_{AB}(r) = \frac{1}{4\pi r^2} \frac{N}{N\_A N\_B \rho} \sum\_{i \in A}^{N} \sum\_{\substack{j \in B, \ j \neq i}}^{N} < \delta \left(r - |r\_i - r\_j|\right) > \tag{2}$$

Figure 1 shows the total pair correlation function *g*(*r*) of liquid Na at 700 K, 1100 K and 1600 K after a 60 ps simulation. The primary peak of *g*(*r*) is located between 3 Å and 4 Å, which agrees well with Li et al's AIMD results [21]. It is interesting that our present simulation shows that there is a small peak before the primary peak at a relatively low temperature, and that this small peak disappears at 1600 K. This small peak was also experimentally observed in the *g*(*r*) function of liquid cesium when the temperature was below 400 K [30]. Figure 1 also demonstrates that the primary peak decreases and broadens obviously as the temperature increases. This phenomenon was also reported by Bickham and his collaborators [31].

**Figure 1.** Total pair correlation function of liquid sodium at different temperature conditions.

The self-diffusion of Na atoms in the liquid phase is also evaluated. The diffusivity (*D*) can be calculated based on the Einstein–Stokes equation:

$$\lim\_{t \to \infty} < \mathbb{R}^2(t) > = 6Dt \tag{3}$$

here, < *R*2(*t*) > is the averaging mean-square displacement of Na. When recording the data, the system is equilibrated for 30 ps first, and then another 30 ps of simulation is run for collecting data.

According to the present study, the Na diffusivities are 13.2 × <sup>10</sup>−<sup>5</sup> cm2/s, 21.4 × <sup>10</sup>−<sup>5</sup> cm2/s and 32.6 × <sup>10</sup>−<sup>5</sup> cm2/s at 700 K, respectively. According to the result of Li et al, the Na self-diffusion coefficient at 723 K was 12.8 × <sup>10</sup>−<sup>5</sup> cm2/s [21], which is close to our present study. Based on experimental data, Meyer fitted an equation to predict the Na self-diffusion coefficient as [32]

$$D = 1.01 \times 10^{-3} \text{ } \exp\left(-\frac{10, 174}{RT}\right) \tag{4}$$

It is worth noting that, in Meyer's work, the experimental data were collected below the temperature of 500 K. According to Equation (4), the Na self-diffusion coefficient is 19.2 cm<sup>−</sup>2/s at 700 K, which is higher than the present and previous AIMD results. However, all computational results yield a magnitude order of 10−<sup>5</sup> cm2/s.

The liquid Na system with Mo or Re atoms is modeled by replacing Na atoms with solute atoms. The effect of the solute concentration on the diffusion and clustering behaviors are investigated in the present study. Here, replacing one Na atom with a solute atom in the 15 × <sup>15</sup> × 15 Å3 computational domain corresponds to a concentration of 0.5 mol/L, replacing two Na atoms with solute atoms corresponds to 1 mol/L and replacing four Na atoms with solute atoms corresponds to 2 mol/L.

Figure 2 shows the partial pair correlation function of one solute atom in the computational domain (0.5 mol/L). It is found that the primary peak of the Mo-Na pair is lower than that of the Re-Na pair. Compared with the Re-Na pair, the primary peak position of the Mo-Na pair also shifts to the right. Therefore, it can be inferred that the solute Re atom can coordinate with more Na atoms than the solute Mo atom. The Bader charge analysis [33] is also performed in order to understand the interaction between the solute and the solvent from the aspect of the electronic structure. It is found that the net charge on the Mo atom is approximately −2.2~−3.0 |e|, and the net charge on the Re atom is approximately −2.8~−3.4 |e|. There is more negative charge transferred from the solvent atoms to the Re atom, which leads to the stronger Re-Na bonds.

**Figure 2.** Partial pair correlation function of 0.5 mol/L solute atoms at the temperature of (**a**) 700 K, (**b**) 1100 K and (**c**) 1600 K.

At the relatively higher concentration, the clustering of solute atoms is observed. Snapshots in Figure 3 show the final structures of the two solute atoms in liquid Na at different temperatures after 60 ps AIMD simulations. It is interesting to find that the solute atoms can stabilize as a dimer in the liquid Na. Figure 4 demonstrates the Mo-Mo distance evolution during the AIMD simulation. The dimer is dynamically stable once it forms, and the average Mo-Mo bond length is 1.79 Å, 1.85 Å and 1.91 Å at 700 K, 1100 K and 1600 K, respectively. Here the Mo-Mo bond length is calculated based on the last 1000 AIMD steps. At 700 K, the dimer forms after 38 ps, and the increase in temperature can speed up the formation of the dimers. It should be noted that Figure 4 cannot reflect the dimerization rate because each case is only run once.

The thermodynamic stability of the dimer is also estimated by the binding energy, which is expressed as

$$E\_b = 2\mathbb{E}\_{M\_1 Na\_{N-1}} - \left(\mathbb{E}\_{Na\_N} + \mathbb{E}\_{M\_2 Na\_{N-2}}\right) \tag{5}$$

Here, *ENaN* is the energy of the pure liquid Na without solute atoms, *EM*1*NaN*−<sup>1</sup> is the liquid Na system with one solute atom, *EM*2*NaN*−<sup>2</sup> is the energy of the system with two solute atoms and the subscript "*N*" represents the number of atoms in the model. All energetic terms are the average value of the last 1000 simulation steps. Table 1 shows the binding energy of dimers at different temperature conditions. The positive binding energy indicates that forming dimers is an exothermic reaction and that dimers are thermodynamically stable at the temperature range of 700 K to 1600 K.

**Table 1.** Biding energy of Mo2 and Re2 dimers at different temperature.


**Figure 3.** Snapshots of two solute atoms (1 mol/L) in liquid Na at different temperatures after 60 ps AIMD simulation: (**a**–**c**) Mo atoms in liquid Na at 700 K, 1100 K and 1600 K, respectively; (**d**–**f**) Re atoms in liquid Na at 700 K, 1100 K and 1600 K, respectively. Violet spheres, yellow spheres and green spheres represent Na atoms, Mo atoms and Re atoms, respectively.

**Figure 4.** Mo-Mo distance evolution for two Mo atoms in the liquid Na system at different temperatures.

It is worth noting that the binding energy increases as the temperature increases. The Bader charge analysis is performed to calculate the net charge on dimers shown in Figure 5. The net charge on Mo2 is −3.6 |e| at 700 K, −2.4 |e| at 1100 K and −2.1 |e| at 1600 K. As discussed in Figure 2, the solute atom captures electrons from surrounding Na atoms, and fewer Na atoms coordinate with the solute atom when increasing the temperature. Therefore, the Mo2 dimer at a higher temperature gains less of a net charge. It is worth

noting that the binding energy of the dimer is proportional to the net charge (Figure 5). For the Mo2 molecule, the lowest unoccupied molecular orbital (LUMO) is empty doubly degenerate *δ*∗ <sup>4</sup>*<sup>d</sup>* antibonding orbitals. The migration of electrons from the Na atoms to the Mo2 dimer will occupy these antibonding orbitals. As the occupation of the antibonding orbitals increases, the Mo-Mo interaction weakens, which leads to a lowering of the binding energy. The net charge of the Re2 dimer is −4.1 |e| at 700 K, −3.9 |e| at 1100 K and −2.3 |e| at 1600 K. As with the Mo2 dimer, the binding energy of the Re2 dimer is also proportional to the net charge. It is worth noting that Re 5*d* orbitals have one more electron than Mo 4*d* orbitals. When forming a neutral Re2 dimer, its *δ*∗ <sup>5</sup>*<sup>d</sup>* antibonding orbitals are halfly occupied, which indicates that the bond of neutral Mo2 is weaker than Re2. In addition, the Re2 dimer can capture more electrons from the Na solvent than the Mo2 dimer. Hence, the occupation of the antibonding orbitals of Re2 in the liquid Na is much higher than Mo2, which results in the much lower binding energy. Taking Re2 in liquid Na at 700 K, for instance, the net charge on this diatomic molecule is −4.1 |e|. In this case, the antibonding *δ*∗ <sup>5</sup>*<sup>d</sup>* orbitals are fully occupied, and the antibonding *π*<sup>∗</sup> <sup>5</sup>*<sup>d</sup>* orbitals are also halfly occupied, which lead to the lowest binding energy of 3.83 eV.

**Figure 5.** The relationship between the binding energy of the dimer and the net charge.

The clustering behaviors of four solute atoms (corresponding to the concentration of 2 mol/L) in the liquid Na are investigated. The solute atoms can be four Mo atoms, four Re atoms or a mixture of two Mo atoms and two Re atoms. Figure 6 shows the total pair correlation function of the liquid system after a 60 ps simulation. For models including four Mo atoms, they still exhibit the typical liquid characterization, as shown in Figure 1. However, obvious peaks located around 2 Å can be found in Figure 6e–i. These peaks can be attributed to the segregation of solute atoms [21]. The atomic configurations of liquid systems after the 60 ps AIMD simulation are shown in Figure 7. It is interesting that the Mo solute and Re solute exhibit relatively different behaviors. Four Mo atoms form two Mo2 dimers, but cannot form a stable Mo4 tetramer. It should be noticed that these two dimers interact with each other and diffuse together (see videos in the Supplementary Material).

For the Re solute, Re4 is observed at 1100 K and 1600 K after a 60 ps simulation. Li et al. also found a cerium tetramer (Ce4) in the liquid Na using the AIMD simulation [21]. In their study, forming Ce4 was attributed to the rejection of the liquid phase to the solute atoms. However, our simulation has demonstrated that there are attractive interactions between solvent Na atoms and Re atoms. It is known that the outmost shell of a Na atom is 3*s* orbital, whereas the outmost shell of a Re atom is 5*d* orbital. The energy level of the former is much lower than the latter. According to the linear combination of the atomic orbitals–molecular orbitals (LCAO-MO) theory, molecular orbitals are always formed by atomic orbitals with a small energy difference. Therefore, the Re-Re attractive interaction is stronger than that of the Re-Mo pair, which leads to the forming of the large Re cluster in the liquid Na. For the mixed solute, the formation of Mo2Re2 is also observed, as shown

in Figures 6g–i and 7g–i. The spin states of dimers and tetramers are also investigated. It is found that the Mo species does not show spin polarization. For the Re species, the magnetic momentum of Re2 is around 0.8 *μ<sup>b</sup>* per Re atom, whereas Re4 does not show an obvious spin polarization.

**Figure 6.** Total pair correlation function of four solute atoms in liquid Na at different systems: (**a**–**c**) four Mo atoms in liquid Na at 700 K, 1100 K and 1600 K, respectively; (**d**–**f**) four Re atoms in liquid at different temperatures; (**g**–**i**) two Mo atoms and two Re atoms in liquid Na at different temperatures.

**Figure 7.** Structures of four solute atoms in liquid Na after 60 ps AIMD simulation: (**a**–**c**) four Mo atoms in liquid Na at 700 K, 1100 K and 1600 K; (**d**–**e**) four Re atoms in liquid Na at different temperatures; (**g**–**i**) two Mo atoms and two Re atoms in liquid Na at different temperatures.

In the present work, the accumulation behavior of six Mo or Re atoms is investigated at 1600 K. It is found that Mo atoms will form three interacting Mo2 dimers, whereas Re atoms form an octahedral-like Re6 cluster.

Understanding mass transport is essential for evaluating the performance and service life of alkali metal high-temperature heat pipes. In the present study, the diffusivities of single solute atoms and small clusters are calculated and shown in Figure 8. As expected, the diffusivity of any species is dependent on the temperature, and a higher temperature leads to a faster diffusion. In addition, the increase in cluster size can reduce the diffusivity. The diffusivities of the Mox (x = 1, 2, 4) clusters are approximately two times as high as those of the Rex clusters, because the atomic weight of Mo (95.94 a. u.) is only half of the atomic weight of Re (186.2 a. u.). The diffusivity of a single Re atom is even lower than the Mo2 cluster. As discussed above, there are not any stable Mo4 clusters formed during the AIMD simulation, but two Mo2 clusters can interact with each other and move together. The interaction between the two Mo2 dimers can enhance the diffusion effects

at low-temperature conditions, as shown in Figure 8. At 700 K, the diffusivity of one Mo2 dimer is 2.53 × <sup>10</sup>−<sup>5</sup> cm2/s, whereas the diffusivity of two interacting Mo2 dimers is 2.81 × <sup>10</sup>−<sup>5</sup> cm2/s. At 1100 K, the diffusivity of two interacting Mo2 dimers is also approximately 0.3 × <sup>10</sup>−<sup>5</sup> cm2/s higher than that of a Mo2 dimer.

**Figure 8.** Diffusivity of Mo and Re clusters in liquid Na at different temperatures.

#### **4. Conclusions**

In this study, AIMD simulations are performed in order to understand the clustering and diffusion behaviors of Mo and Re atoms in liquid Na. For a single solute atom, the Re-Na interaction is stronger than the Mo-Na interaction. For models with two solute atoms, both Mo2 dimers and Re2 dimers are found in the Na solvent. The binding energy of the dimer is proportional to the occupation of the antibonding orbitals of the dimer. The higher occupation always leads to a lower binding energy. For models with four solute atoms, four Mo atoms form two interacting Mo2 dimers, but a stable Mo4 tetramer is not observed. However, for the pure Re solute or mixed solute, a stable tetramer is observed. The diffusivity of the Mo species is much faster than that of the Re species, and the diffusivity decreases as the cluster size increases.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/met11091430/s1, Video S1: Trajectory of four Mo atoms in Na at 1600 K; Video S2: Trajectory of four Re atoms in Na at 1600 K.

**Author Contributions:** Z.L. and M.M. concept this work; Z.L. performs simulations, collects data and writes the original draft; all authors analyze results and polish the manuscript. W.L. and H.D. have contributed writing–review and editing manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The President Foundation of China Academy of Engineering Physics (No. YZJJLX2018002); Fundamental Research Funds for the Central Universities.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work is financially supported by the President Foundation of China Academy of Engineering Physics (No. YZJJLX2018002). Z. Liu also thanks Fundamental Research Funds for the Central Universities.

**Conflicts of Interest:** The authors declare no conflict of interest.
