Computer Simulation of Liquid Pu

1. Introduction

Direct measurements of liquid Pu properties are complicated by the metal's radioactivity and high chemical activity. Most of the properties of liquid Pu were determined by computer modeling, including the ab initio calculation method. The interparticle potential was determined in [1] in the MEAM form. The number of properties was investigated using computer calculations and by applying the theory of liquids. Some properties of liquid Pu were nevertheless measured experimentally. These are density, energy, heat capacity, speed of sound, compressibility, viscosity, and surface tension.

The density of liquid Pu was measured in the range up to 1200 K in [2,3]:

d = 17.567 – 1.451^.10^-3(T-273) g/cm³

(1)

At 950 K, this formula yields a density of 16.58 g/cm³. According to the reference book [4], at 938 K the figure obtained is close to 16.50 g/cm³.

Bearing in mind the physical structure of EAM and MEAM, we must calculate the energy U of the models with respect to a system of non-interacting atoms at rest. The energy of such a system can be chosen equal to zero. The energy of a substance is easy to calculate if its heat of evaporation ΔH is known. To calculate the energy U, it should be accounted that during the transition to the ideal gas state ΔU = ΔH – RT. Subtracting the kinetic energy of an ideal gas (3/2)RT, we find the change in energy during the transition of a substance to the state of non-interacting atoms at rest ΔU = ΔH – (5/2)RT. Consequently, the energy of a substance with respect to a gas at rest (i.e., to the origin) is equal to U = - ΔH + (5/2)RT. These estimates are not entirely accurate, since the properties of Pu gas deviate from the properties of an ideal gas [5]. According to [5,6], the heat of vaporization of Pu at the melting point is 336.6 ± 1.42 kJ/mol. Hence, the energy of liquid Pu near the melting point T_m is U = -317.7 kJ/mol. The report [5] provides data on the enthalpy of Pu at temperatures up to 2500 K.

An important characteristic of liquid Pu is its heat capacity C_p = 41.2 J/mol/K at the melting point T_m, which is significantly higher than the classical value 3R = 24.94 J/mol/K [5]. Already at 298 K, the heat capacity of solid Pu (35 J/mol/K) exceeds the classical value. This excess is due to the heat capacity of conduction electrons and the internal degrees of freedom of atoms. Monatomic Pu gas in the range of 400 - 500 K has a heat capacity of 2.76R, which is close to the classical value of (5/2)R. However, already in the range of 900-1000 K, the heat capacity of the gas is 3.56R, and in the range of 2400-2500 K C_p = 7.4 R [5]. Here, the internal degrees of freedom of atoms must be considered in calculations.

In the work [7] measurements of electrical resistance, thermal expansion, enthalpy, compressibility, Grüneisen coefficient, sound speed, and heat capacity ratio C_p/C_v were carried out at temperatures up to 4000 K. An increase in sound speed was noted at temperatures up to 2000 K, and a decrease in sound speed with further heating. The reason considered is the delocalization of 5f electrons. Accordingly, the bulk modulus and Grüneisen coefficient change with a maximum at 2000 - 2200 K. With a decrease in density, the C_p/C_v ratio increases (up to 1.40 at 4000 K).

All calculations of the self-diffusion coefficients of liquid Pu were performed using the molecular dynamics (MD) method. In the case of shear viscosity, both a conventional experiment (damping of torsional oscillations of a cylinder with liquid) and the MD method were used. Data on self-diffusion [8,9] and shear viscosity [9,10,11,12,13] have been published. In [9], the self-diffusion coefficients and viscosity of Pu were calculated at temperatures up to 1400 K. The surface tension σ of liquid Pu was measured in [14]. In the range of 913-1200 K, the value of σ = 550 – b (T – 913) mN/m, where b = 0.08-0.15.

In this work, a model of liquid metal is constructed based on the known pair correlation function (PCF) and the density of liquid Pu. The EAM interparticle potential is reconstructed and computer modeling of Pu is performed at temperatures of 500 - 5000 K.

2. Materials &Methods

Interparticle potential of Pu. In the early work of Wittenberg [15] the hard sphere model was applied. The compressibility of the metal could be correctly described, but viscosity and surface tension values were underestimated. Then, the potential of the modified embedded atom model (MEAM) was proposed for liquid Pu [1]. The parameters of this potential were found (but not given) considering the known properties of the liquid phase and six crystalline phases of Pu. This potential was applied for molecular dynamics (MD) modeling of Pu. It allows obtaining agreement with experiment for some important properties of the models: energy (3.80 eV/atom), melting temperature (1000 K), bulk moduli for α, β and γ phases of Pu. The decrease in volume upon melting was also confirmed.

Further analysis of the properties of liquid Pu was carried out using the MEAM potential [8]. The models contained 2048 atoms in the main cube. The pair correlation function (PCF) of liquid Pu at 950 K, as well as the self-diffusion coefficient D and shear viscosity η at 900 – 1300 K were calculated. At the melting temperature (913 K), the values of D = 2.3^.10^-6 cm²/s and η = 2.52 mPa^.s were found to be significantly underestimated in relation to the experimental values (D ~10^-5 cm²/s and η = 6.30 mPa^.s [12,13]). Apparently, the MEAM potential [1], calculated under conditions of insufficient initial data, required clarification.

A case of successful calculation of the interparticle potential at a significant lack of experimental data can be considered in the search of interparticle potential of liquid uranium [16]. Pair contribution to the EAM potential was chosen as the Morse potential, and the parameters of the embedding potential were selected using the dependence of the density of liquid on temperature. In this paper, an attempt to construct a model of a liquid metal based on the known pair correlation function (PCF) and density was realized. The interparticle EAM potential was calculated and computer modeling of Pu was carried out at temperatures of 500 - 5000 K.

Interparticle potential EAM. Several methods for constructing models of liquids and amorphous substances based on known diffraction data (structural factors or pair correlation functions) were previously proposed [17]. The main results were previously obtained for simple liquids with short-range pair potentials. A similar problem is presented by the case of Pu. To calculate the interparticle potential of liquid Pu, it was decided to use its pair correlation function (PCF) near the melting point, calculated in [8]. We used the Schommers algorithm [18], which allows constructing a liquid model with simultaneous restoration of the interparticle potential.

There are few data on the structure of liquid Pu. No direct diffraction experiment could be found due to the difficulties of working with fissile materials. However, the pair correlation function g(r) (PCF) of Pu was calculated by the MEAM method at 950 K [8] (r is the interparticle distance). We scanned the PCF graph of liquid Pu [8] at 950 K (see Figure 1). Hereinafter, we will call this function PCF-1. It has the usual appearance for simple liquids, with the height of the first peak slightly less than 3.0, and with a coordinate of 3.175 Å. We denote the particle number density by n = N/V (N is the number of particles in a volume V, n = 0.04091 at/ų). If one plots the dependence of g(r) on the value of rn^1/3 at 950 K, the graph for Pu almost coincides with the case of liquid Al at 943 K [8].

Additional information on the structural characteristics is provided by checking Pu for the feasibility of the main topological criteria of the PCF shape. The structural features of a liquid metal can be analyzed using relationships describing the shape of the PCF of simple liquids. As shown in [20], in the case of simple one-component liquids with pair interaction, there is a relationship between the coordinate of the 1st peak of the PCF r₁, the height of this peak g(r₁), the minimum interparticle distance r_min in this liquid and the volume per 1 atom v₁ = V/N. We define the quantities d = (V/N)^1/3 (N is the number of atoms in volume V), ρ₁ = r₁/d, ρ_min = r_min/d. The maximum possible value of ρ_min in one-component disordered systems with central interaction is 1.0673 [20]. We introduce the notation y = ρ_min^-1 - 1.0673^-1. Analysis of 32 computer models with different pair potentials, and real structures of liquid metals allows us to obtain the relationship [20]:

g_c(ρ₁) = 1 + 0.600/y^0.77 + 0.0001/y³

(2)

Deviations from this relationship are small and on average do not exceed ±5%. Thus, the ratio of the peak heights S = g(ρ₁)/g_c(ρ₁) for a real PCF to the "theoretical" value should be S = 1 ± 0.05 in the case of simple liquids and amorphous substances. Deviations to the smaller value should be topologically forbidden, and to the larger value - can be allowed for systems with non-central interaction. Note that the left (descending) branch of the 1st PCF peak of the system, i.e. the minimum interparticle distance, plays a large role in calculations. This characteristic of the PCF plays an important role in the construction of models of liquid metals, but it is the least known from experiment.

From this position, we will consider the PCF-1 of liquid Pu obtained by the MEAM method in [8]. The parameters of this PCF are given in Table 1. Judging by the values of S from Table 1, the structure of liquid Pu with such parameters as those of the PCF [8] at a temperature of 950 K, is topologically realizable, but a centrally symmetrical structure with such PCF cannot really exist, since the value of S (i.e., the height of the 1st peak) is too large. Consequently, it is impossible to construct a model with central interaction and with exactly the same PCF.

Nevertheless, we can try. We can try to build a model with the same target PCF-1 as curve 1 in Figure 1. Let's call this virtual substance para-Pu. For this purpose, we applied the Schommers algorithm [18]. The potential in this algorithm is calculated iteratively during the calculation process. At each iteration with a trial potential, the difference between the obtained PCF and the specified (target) PCF is calculated, and a correction to the current pair potential is determined in the form of a table.

Let us denote by g⁰(r_i) the diffraction ("target") PCF, and by g(r_i) the PCF of the model constructed with some trial potential φ(r_i). Then a new approximation of the potential is calculated, according to the Schommers algorithm [18], by the formula:

φ'(r_i) = φ(r_i) + α ln[g(r_i)/g⁰(r_i)]

(3)

Here α is an empirical coefficient regulating the rate of convergence of the results to the asymptotic limit. Using potential φ'(r_i), a new model of liquid is constructed, new PCFs g(r) are found, etc., until the difference between the diffraction and model PCFs becomes small enough. Since in the molecular dynamics method the PCFs are calculated as histograms, the discrepancy between them can be found numerically as standard deviations ("residuals"):

$${\mathrm{R}}_{\mathrm{g}}=\left\{{\displaystyle \frac{1}{{\mathrm{n}}_2-{\mathrm{n}}_1+1}\sum _{{\mathrm{n}}_1}^{{\mathrm{n}}_2}}\right[{\mathrm{g}}^0\left({\mathrm{r}}_{\mathrm{j}}\right)-\mathrm{g}\left({\mathrm{r}}_{\mathrm{j}}\right){]}^2{\}}^{1/2}$$

(4)

Here g⁰(r_j) is the histogram of the target PCF, g(r_j) is the histogram of the model PCF, n₁ and n₂ are the summation boundaries of the tabular data, and j is the number of the histogram element. As n₁, the number of the first non-zero element of the PCF was chosen, and n₂ was determined by the cutting radius of the interaction. If the value of R_g is several hundredths, then the graphs of the functions g⁰(r) and g(r) are visually almost indistinguishable. Therefore, the iterative process can be carried out until the residual values of the order of 0.01 are reached.

3. Results

The initial model of the Schommers algorithm at 950 K contained 2048 particles in the main cube with an edge length of 36.751 Å (density N/V = 0.04126 at/ų or 16.58 g/cm³). The copper potential was used as a seed [21]. At distances less than r_min = 2.40 Å, the values of Pu potential are not determined in the Schommers procedure and were chosen empirically in the form of an ascending branch:

φ(r), eV = 1.82336 + 0.025(r_min - r)² + 5.35964(r_min - r) + 25.03(r_min - r)^2.20.

(5)

The Pu model (M1) at 950 K was constructed in 140 Schommers iterations using the PCF-1. The details of the modeling are described, for example, in [17]. In this algorithm, the interparticle potential φ₁(r) was calculated in the form of a table with a step of 0.05 Å. The characteristics of the structure of this model are given in Table 1. The deviation of the parameter S of the target PCF-1 from the topological standard led to the same deviation for the PCF of the M1 model. The discrepancy between the PCF of the M1 model and the PCF-1 Pu was R_g = 0.078, which could not be reduced.

And that's not all. It turned out that the models constructed by this method practically do not expand when heated to 2000 K. For example, the model had a density of 16.543 g/cm³ at 950 K, and 16.514 g/cm³ at 2000 K. In this case, the coefficient of volume expansion is 1.67^.10^-6 K^-1, an order of magnitude smaller than that of ordinary metals. This means that the resulting potential is very close to a symmetrical parabola. This solution to the problem and the potential φ₁(r) itself had to be rejected.

The second variant was based on the indication of the authors [8] that the PCF of liquid Pu at 950 K and aluminum at 943 K are very similar if the rn^1/3 values are plotted on the abscissa axis (n = N/V is the density). In the case of Pu at 950 K, the density n = 0.04126 at/ų, and for Al it is slightly higher: n = 0.06029. The structural parameters of aluminum at 943 K are quite acceptable; the height of the 1st peak of the PCF is 2.828 and the coefficient S = 0.9899 is very close to unity. Therefore, we constructed a model of Pu using the PCF of aluminum at 943 K [19]. To do this, it was necessary to take the aluminum PCF, stretch its abscissa axis by 0.06029^1/3 = 0.3921 times and then reduce it by 0.04126^1/3 = 0.3455 times, that is, stretch it by 1.13488 times. The resulting PCF (PCF-2) is shown in Figure 1 and 3.

Further calculations consisted of reconstructing the interparticle potential φ₂(r) of a liquid metal using the Schommers algorithm, which generated the potential φ₂(r) and PCF-2 function. Since the PCF-2 is no longer a function of the natural Pu, although close to it, the resulting interparticle potential will be a certain approximation to the exact potential.

Figure 2. PCF of Pu at 950 K. Series 1 – diffraction PCF-2, calculated using the PCF of aluminum at 943 K. Series 2 – PCF of the M2 model of liquid Pu, constructed at 950 K via Schommers algorithm.

Figure 2. PCF of Pu at 950 K. Series 1 – diffraction PCF-2, calculated using the PCF of aluminum at 943 K. Series 2 – PCF of the M2 model of liquid Pu, constructed at 950 K via Schommers algorithm.

Using the Schommers algorithm, the spherically symmetric pair potential φ₂(r) was calculated, and a liquid metal model M2 was constructed at 950 K. The discrepancy between the PCF M2 and PCF-2 is 0.1074. The density of the model was 16.58 g/cm³ and is exactly equal to the data (1). The pressure of this model is 0.04 GPa and is quite close to zero. The energy of the model (-317.15 kJ/mol) is almost exactly equal to the energy of real Pu (-317.7) [5]. Table 2 shows the values of the pair potential φ₂(r) of the Pu model, which was then used as the pair contribution φ(r) to the EAM potential for modeling Pu. The potential graph is shown in Figure 3.

The characteristics of the PCF-2 structure are given in Table 3. The above-described method of recalculating the PCF of aluminum to the PCF-2 leads to an underestimated S value of the target structure, so that difficulties could arise with constructing a computer model. However, the structure of the M2 model constructed by the Schommers algorithm partially corrects the shortcomings of the calculated PCF-2. The discrepancy between the PCFs of these models is quite large (0.1071), while in the case of good results of the Schommers algorithm it can be 0.01-0.04. As can be seen from Figure 2, this discrepancy is due to a very small shift of the sharp 1st peak of PCF along the abscissa axis.

The Embedded atom model (EAM) potential. In this paper, EAM potential [22] was applied. It has the form:

$${\mathrm{U}}_{\mathrm{pot}}{\displaystyle \sum _{\mathrm{i}}\mathsf{\Phi }\left(\mathsf{\rho }\mathrm{i}\right)}+{\displaystyle \sum _{\mathrm{i}<\mathrm{j}}\phi \left({\mathrm{r}}_{\mathrm{ij}}\right)}$$

(6)

Here U_pot is the potential energy of the system, Φ(ρ_i) is the embedding potential of the i-th atom, depending on the effective electron density ψ at the location of the center of the atom, and the second sum over pairs of atoms - the pair contribution - contains the usual pair potential calculated above by the Schommers algorithm. The effective electron density at the location of the atom is created by the surrounding atoms and is determined by the formula:

$${\mathsf{\rho }}_{\mathrm{i}}={\displaystyle \sum _{\mathrm{j}}\mathsf{\psi }\left({\mathrm{r}}_{\mathrm{ij}}\right)},$$

(7)

where ψ(r_ij)) is the contribution to the effective electron density from the neighbor number j. It is convenient to choose the function ψ(r) in the form:

ψ(r) = p₁ exp(- p₂r)

(8)

Embedding potential is taken as piecewise function of the form:

Φ(ρ) = a₁ + c₁ (ρ - ρ₀)²               at ρ₁ ≤ ρ ≤ ρ₄,
Φ(ρ) = a_i + b_i (ρ - ρ_i-1) + c_i (ρ - ρ_i-1)²       at ρ_i ≤ ρ ≤ ρ_i-1 (i = 2,3),
Φ(ρ) = [a₄ + b₄(ρ - ρ₃) + с₄(ρ - ρ₃)²] [2ρ/ρ₃ - (ρ/ρ₃)²]       at ρ ≤ ρ₃,

(9)

where ρ₀ = 1, and for ρ = ρ_i the function Φ(ρ) itself and its first derivative are continuous. The function φ(ρ) and all coefficients a, b and c are expressed in eV. The coordinates of the division points of the abscissa axis increase in the sequence ρ₃ < ρ₂ < ρ₁ < ρ₀ < ρ₄ . As a result, the EAM potential is determined by the parameters p₁, p₂, a₁, c₁ – c₄, ρ₁ – ρ₄, m, n, which allow, in principle, to fit the properties of the models to the selected experimental data. This fitting is carried out based on the dependence of the density and energy of the metal on the temperature along the binodal, as well as on the data of static and shock compression. Expressions for ρ < ρ₀ are used when modeling states with normal and reduced density, and for ρ > ρ₄ – for compressed states. We don’t consider here the compressed states. The parameter p₂ in (8) is adjustable. The parameter p₁ was determined in such a way as to obtain the average value <ρ> = ρ₀ = 1 for the liquid model in the “standard” state (near the melting point). With this value of <ρ>, the embedding potential has almost no effect on the motion of particles, since dΦ(ρ)/dρ is close to zero at ρ ≈ 1. Coefficient a₁ is determined through the energy of the model in the "standard" state at ρ = 1, and the coefficient c₁ is determined through the bulk modulus in this state. The coefficients c₂ – c₄ are selected based on the dependence of the metal density on temperature along the binodal. The coefficients a₂–a₄, b₂–b₄ are calculated from the continuity condition of the embedding potential and its derivative at the points ρ_i. Note that in the state at ρ = ρ₀ = 1 the contribution from the embedding potential disappears, and the properties of the liquid are determined only by the pair contribution.

The characteristics of the Pu potential are given in Table 4. The coefficients p₁ and p₂ were selected using the condition that in the "standard state", namely near the melting point, the average value of the effective electron density of atoms <ρ> is equal to unity. The coefficients c_i can be determined by dependence of the liquid density on temperature. Considering the values of density and energy of real Pu, the coefficients of the expression for the EAM potential can be determined (see Table 4).

The Pu models with a size of 2048 particles in the main cube were constructed using the MD method with a time step of 0.0022t₀, where t₀ is an internal time unit equal to 1.5904^.10^-13 s. The coefficient 0.0022 was selected from the condition of minimal drift of the total energy of the system in the NVE mode. Particle displacements were calculated using L. Verlet algorithm. The potential cutoff radius was 9.40 Å. The results of the simulation liquid Pu are given in Table 5.

Figure 4 shows the dependence of the density of Pu models on temperature. This dependence is linear at T > T_m, but in the supercooling region, where T < 913 K, upward deviations are observed.

The calculated heat capacities C_v and C_p are given in Table 6. Our computer calculation slightly overestimates the ratio C_p/C_v (by ~0.02), but generally dependence of the heat capacities of models and real Pu on temperature is the same.

According to [7], the sound speed during heating passes through a maximum at 2000 K. At liquid densities above 15 g/cm³ (i.e. at 900 - 2500 K, see Table 5), the bulk modulus K_T of liquid Pu changes little and is close to 24 GPa, and with further heating it decreases with a slope of 3.6 GPa/(g/cm³). We calculated the isothermal modulus K_T using the dependence of the model pressure on the volume at a constant temperature in the NVE mode. The values of the modulus at temperatures up to 5000 K are given in Table 5. Our calculations show that the dependence of the modulus on temperature really changes little at T < 2500 K, and then begins to decrease at heating in accordance with experiment. These features can be explained by the proposition that 5f electrons of Pu play the role of valence electrons because of delocalization, and the increase in the sound speed during heating (maximum at 2500 K) is caused by the hybridization of 5f electrons with 6d7s electrons of the conduction band of Pu [7].

The self-diffusion coefficients D of liquid Pu were measured from the time dependence of the mean square atomic displacements. The values of the coefficient D at temperatures of 913-5000 K are given in Table 7 and shown in Figure 5, where the notation D5 = D*10⁵ is adopted. They are well straightened in logarithmic coordinates and approximated by the equation:

ln D5 = 1.5019lnT - 10.063, or D = 4.263^.10^-10 T^1.5019 cm²/s.

(10)

Equations of this type are well satisfied for liquid metals, but the exponent at temperature is usually close to 2.0 (2.0047 for Rb [23], 2.2109 for Cu [24], 1.7411 for Tl [25], etc.). For example, in the case of liquid uranium, the expression D = 5.17 × 10^–12T^2.1029 cm²/s was obtained [16].

The graph in Figure 5 shows that upon transition to the supercooling region (T < 913 K), the coefficient D begins to deviate downwards from expression (10). The models seem to sense that they have entered the supercooling region. In the theory of simple liquids, sensitivity of this kind is not expected, since the models are actually one-component. Analysis of the structure of the models revealed only a gradual decrease in the coordination number upon heating from 13.33 at 500 K to 10.64 at 5000 K. In this case, the coordinate of the 1st peak shifts from 3.14 to 3.12 Å.

This data can be compared with the diffusion characteristics of hard-sphere liquids. The effective diameter σ of a hard sphere can be found from the condition that the packing fraction of hard-sphere system (HSS) η near the melting point of liquid metals is ~0.45. For Pu at 913 K, the condition η = (N/V) (πσ³/6) = 0.45 yields σ = 2.75 Å. This value is slightly smaller than the usual atomic diameter of Pu (3.02 Å). In the case of liquid Si and Ge, the interparticle potential is indeed close to the hard-sphere potential with a step down [26]. This allows us to check the applicability of the formula for the self-diffusion coefficient of the hard-sphere model [27]:

D_HS = (D₀/n) (1-n/1.09) (1+n²(0.4 – 0.83 n²)),

(11)

where n = (N/V) σ³, N is the number of particles in volume V, σ is the effective diameter of the sphere, m is its mass, D₀ = (3/8) σ(kT/πm)^1/2. Multiplying the numerator and denominator of the fraction in brackets by Avogadro's number, we obtain the following at T = 913 K: σ = 2.75 Å, N/V = 0.041354 at/ų, n = 0.8591, D₀ = (3/8) 2.279×10^–8×(8.314×10⁷×913/π/243.06)^1/2 = 1.028×10^–4 cm²/s, and finally, D_HS = 2.13×10^–5 cm²/s. This value is in poor agreement with the figure of 1.16×10^–5 obtained by the MD method. The reason is that, unlike Si and Ge, the pair contribution to the interparticle potential of Pu (Figure 3) differs significantly from the hard-sphere one. The value of D found by the MEAM method (2.3×10^–6 cm²/s at 913 K [8]) is apparently greatly underestimated.

Table 7 presents the literature data on the self-diffusion of liquid Pu obtained theoretically in comparison with our calculations. Quantum-mechanical calculations [9] were carried out using the VASP package. These data are overestimated in relation to ours by about 1.5 times, but upon heating the difference between them decreases.

It was difficult to calculate the viscosity of liquid Pu using the Green-Kubo equations [28] due to the large scatter of results. However, the steady-state flow method can be used [29]. If an external force of the form:

F_z = F₀sin(2πx/L),

(12)

acting on each atom is included in the molecular dynamics model located in the main cube with periodic boundary conditions, then a steady-state double flow of liquid along the z axis should arise in the model, and the velocity profile should also have the form v_z = v₀ sin(2πx/L). The solution to the hydrodynamic problem has the form [29]:

$$\mathsf{\eta }={\displaystyle \frac{N{L}^{2 }{F}_0}{4{\pi }^{2}V{v}_0}}$$

(13)

Here η is the dynamic viscosity of the liquid, N/V is the particle number density, L is the edge length of the main cube, F₀ and v₀ are defined above. In our case, with a model size of 2048 atoms in the main cube with periodic boundary conditions, a steady-state laminar flow should arise along the z axis at coordinates 0 < x < L/2 and in the opposite direction at L/2 < x < L. For verification, the main cube was conditionally divided into 10 sections with coordinates 0 < x < 0.1L, 0.1L < x < 0.2L, etc., and in each section the average value of the projection of the velocity of atoms <v_z> in the MD run was determined. If the force F₀ had an order of magnitude of 0.01 eV/Å, then the profile of the average velocity <v_z> in z direction had the form of a sinusoid. An example of the distribution of velocities v_z for the case with F₀ = 0.01 eV/Å is given below:

Layer number:	1	2	3	4	5	6	7	8	9	10
1000<vz>, Å/Δt:	2.6	4.3	6.5	5.7	2.8	-2.8	-7.0	-5.9	-4.4	-1.9

Here Δt is the MD time step: Δt = 0.01t₀, where t₀ = 1.5836^.10^-13 s is the internal time unit of the Pu models. Indeed, the dependence of v_z on x has a sinusoidal shape. When the external force was reduced to F₀ = 0.002 eV/Å, the velocity profile became chaotic. The MD program calculated the average value of the particle velocity projection <v_z> with coordinates 0 < x < L/2 and found the viscosity η using formula (13), considering that for a sinusoidal profile the value v₀ = (π/2) <v_z>. Typically, the viscosity was calculated on 10-15 MD runs with 5000 steps each. The standard deviation of the viscosity value from the average was ~9%. With the specified modeling parameters, there is no noticeable heating of the model with normal temperature control, and the pressure remains acceptable in value with fluctuations of about 0.01 GPa.

The steady-flow method was used to calculate the dynamic viscosity of Pu at temperatures up to 5000 K (see Table 7). Our data are in good agreement with measurements [12,13] up to 1200 K, but decrease with further heating more slowly than the calculations in [9,10,11]. The values of Pu viscosity at temperatures up to 5000 K are shown in Figure 6. The graph also shows the viscosity values in the supercooled state (T < 913 K).

At temperatures of 913-5000 K, the viscosity can be described by an expression like the formula for the self-diffusion coefficient (10):

Ln η, mPa^.s = -0.6351lnT + 6.1477, or η, mPa^.s = 467.64*T^-0.6351

(14)

The viscosity of supercooled metal, like self-diffusion, is sensitive to entering the supercooling region.

The MD calculations with the EAM potential can be assessed using the Stokes-Einstein Equation, which relates the self-diffusion coefficient D to the viscosity η. It is usually valid for simple liquids [30,31], including liquid metals [32]:

$$D= {\displaystyle \frac{kT}{4\pi \eta r_a}}$$

(15)

Here r_a denotes the effective hard-sphere radius of the atom. Taking D = 1.16^.10^-5 cm²/s and r_a = σ/2 = 1.375^.10^-8 cm, we find an estimate of the viscosity of Pu at the melting point: η = 6.30 mPa^.s. This value coincides well with the viscosity of Pu of 6.30 mPa^.s, measured near the melting point [10,12,13,14].

If the formula (15) is valid, the relation Dη/T = const should be satisfied. Using the results for the values of D and η obtained in the simulation with the EAM potential, the dependence of the values Dη/T on temperature can be calculated. They are given in Table 7 (the values of Z₁ = 10¹⁰Dη/T in units of g/(cm^.s²K). The values of Dη/T depend rather weakly on temperature, so that the Stokes-Einstein relation (15) is fulfilled quite well for liquid Pu in the entire range of 913 - 5000 K.

In the work [33], the correction to the Stokes-Einstein equation is proposed, in which this equation takes the form DηV^1/3/T = const, where V is the volume [9]. The values of Z₂ = Z₁(V/V₉₁₃)^1/3 are also given in Table 7. For Pu, Z₂ = 7.31 ± 0.52, i.e. the deviations do not exceed 7%, and the Stokes-Einstein equation with this correction is fulfilled very well. In theoretical calculations, an underestimation of viscosity is observed at elevated temperatures [9,11].

Additional information about liquid Pu is provided by calculating its surface properties. It is possible to calculate the excess surface energy of nanoclusters (not surface tension!) by constructing a series of them with different sizes. Following the work [34], we write the energy of a spherical nanoparticle of radius R as:

E = aN + bN^2/3, or E/N = a + bN^-1/3

(16)

where N is the number of atoms in the cluster. Here bN^2/3 is considered generally as the excess surface energy. Dividing it by the surface area, we obtain the specific surface energy σ_E = bN^2/3/4πR². The number of atoms in a spherical drop is (4/3)πR³/v₀, where v₀ = V_mol/N_A (V_mol is the molar volume, N_A is Avogadro's number). Hence, σ_E = (4π/3)^2/3(1/4π) b/v₀^2/3.

To calculate the values σ_E, a series of Mackay icosahedral nanoclusters with a free surface and particle numbers of 13, 55, 147, 309, 561, 923, 1415 and 2869 were constructed and simulated at 950 K. Their energies are shown in Figure 7. The value at N → ∞ (-3.2919 eV/atom) was obtained by conventional modeling of liquid Pu. For all N >13 the graph is perfectly straightened, and b = 5.8155 eV/at^2/3. Considering that at 950 K V_mol = 14.60 cm³/mol, we obtain as a result σ_E = 2.30 J/m². This value is much greater than the real surface tension ~0.55 J/m² [14]. It means that the term bN^2/3 is not common excess surface energy but contains also the contributions from the bulk nanoparticle atoms.

4. Discussion

This paper presents an attempt to construct a model of liquid metal having only its pair correlation function (PCF) and density. Analysis of the topological characteristics of liquid metal allowed us to determine, to what extent this PCF can correspond to the structure of a real liquid.

The procedure for calculating the pair potential using the Schommers algorithm is stable and allows us to confidently reconstruct the interparticle potential via the shape of the PCF. In the case of Pu, a testing potential was previously discovered at which the expansion coefficient of liquid metal turned out to be abnormally low – an order of magnitude less than the usual values. This potential was calculated using the Schommers algorithm and the original PCF Pu obtained in [8]. For example, in calculations with this potential in the NpT mode, the densities of the liquid at temperatures of 950 and 2000 K were equal to 16.543 and 16.514 g/cm³, respectively, i.e. the coefficient of volume expansion α was equal to 1.67^.10^-6 K^-1. For real Pu in this range, the coefficient α = 6.18^.10^-5 K^-1, i.e. 40 times greater.

The Pu potential obtained using the structure of liquid aluminum turned out to be more accurate [8]. Several properties of Pu models - density, heat capacity, compressibility - showed a rather weak dependence on temperature. With heating to 5000 K, only a small increase in C_p to ~30 J/mol/K is observed. This agrees with the data of [5], where the heat capacity of real liquid Pu does not depend on temperature at all. Consequently, the small excess heat capacity of the Pu models (in relation to the classical value 3R) is due to the shape of the potential, and the high heat capacity of real Pu is due to the internal degrees of freedom of the atoms. The contribution from the conduction electrons is apparently small, since their heat capacity would be visible over a wide temperature range. The calculated heat capacity ratios C_p/C_v (Table 6) are in good agreement with the experimental data [7].

Three important properties of liquid metal – density, diffusion rate, and viscosity – deviate from high-temperature trends in the same low-temperature ranges. Using the EAM potential, it was possible to calculate the values of the self-diffusion coefficients in the entire range from 500 to 5000 K. It should be noted that the order of magnitude of Pu self-diffusion coefficients is typical for metals (D ~ 10^-5 cm²/s near the melting point). Theoretical calculations [8,9] do not agree with our results. The temperature dependence D(T) at T > T_m is well approximated by the expression D = aT^b, which follows from the R. Swalin theory [35]. Arrhenius-type dependences do not reflect the physical picture here, since the act of diffusion in simple liquids does not require activation.

The results of calculations of the dynamic viscosity of Pu by the laminar flow method look promising [29]. Here, good agreement with experiment was obtained for temperatures near T_m [12,13]. In the supercooling region, the viscosity increases to ~20 mPa^.s, i.e. an order of magnitude higher than the usual values for liquid metals. Real Pu behaves as a fairly viscous liquid. Liquid Pu models are also rather viscous and, in addition, are prone to strong supercooling, so they easily retain the properties of a liquid even when supercooled by 400 K. At temperatures above T_m, the Stokes-Einstein relation (15) with a correction [33] is quite well satisfied. The behavior of self-diffusion coefficient and viscosity at temperatures below the melting point is noteworthy. These properties seem to indicate the fact of transition to a supercooling state, although the theory of simple liquids does not attach special importance to entering a two-phase region in the actual absence of a second phase. It seems quite probable that characteristic structure changes, such as signs of dimerization, could appear in the supercooling region of real Pu and its models. However, the pair correlation functions of the Pu models at 500-913 K do not show any special deviations from the normal form. As an example, we can cite the PCF of supercooled Pu at 500 K (see Figure 8).

This graph has a form typical for liquid metals, and the sign of supercooling is only the increased height of the 1st peak, equal to 3.63.

Analysis of the properties of real Pu and its models allows us to assume that the models with the EAM potential constructed in this work adequately describe some important properties of real Pu.