Effects of system-size and inner-core 2 p states on melting of dense sodium at high pressure : ab initio molecular-dynamics simulation

We report the effects of system size on melting of dense sodium at about 100 GPa. We have performed systematic ab initio molecular-dynamics (MD) simulations with the number of atoms at most 250 for bcc and 256 for fcc structures based on 3s-valence-electron model and 2p3s-valence-electron model. We have shown that the inner-core 2p states play important roles at extremely high pressures and that the melting temperature estimated by one-phase approach depends on the system size; in particular, the melting temperature for 3s model decreases with increasing system size and deviates from the experimental value, while the melting temperature for the 2p3s model approaches to the experimental value slowly. This result suggests that we have to carry out large-scale ab initio MD simulations with the inner-core 2p states to obtain the melting temperature of dense Na at extremely high pressure.


Introduction
Since Gregoryanz et al. [1] have performed x-ray diffraction experiments on sodium (Na) at high pressures of up to 130 GPa and have shown that there exists a maximum of the melting temperature at 30 GPa and a minimum of the melting temperature, which is unusually low temperatures, reaching 300 K at 118 GPa, some ab initio moleculardynamics (MD) simulations [2][3][4][5][6] have been performed on dense Na.
Hernández and Íñiguez [2] have performed ab initio MD simulations, treating not only the 3s electron but also the 2p electrons as valence electrons to avoid an artifi cial dimerization of Na, and suggested the existence of the melting point maximum from the compressibility, though they did not obtain the melting curve.
Raty et al. [3] have performed ab initio MD simulations treating only the 3s electron as valence electron (hereafter referred to as 3s-valence-electron model or 3s model) to obtain the melting curve of dense Na, and shown that their simulations can explain the unusual melting behaviour in dense Na.
Koči et al. [4] have carried out ab initio MD simulations treating 2p as well as 3s electrons as valence electrons (hereafter referred to as 2p3s-valence-electron model or 2p3s model) and shown that there exists a maximum melting point around 40 GPa, though the obtained melting temperature starts to deviate higher than experimental values at extremely high pressures larger than 50 GPa.
We [5,6] have also performed ab initio MD simulations based on two different models, 3s-valence-electron model and 2p3s-valence-electron model, and have shown that the a e-mail: khoshino@hiroshima-u.ac.jp difference of these models are significan at extremely high pressures, since inner-core 2p states play important roles at extremely high pressures, where the 2p wavefunctions of neighbouring Na atoms overlap due to the shortening of average nearest-neighbour interatomic distance.Our melting curves obtained based on 3s model and 2p3s model are similar to those obtained by Raty et al. [3] and Koči et al. [4], respectively.
In Fig. 1 we show the melting curves obtained experimentally and theoretically using ab initio MD simulations.The melting temperatures obtained by 2p3s-model is higher than those obtained by 3s-model and the experimental values at pressures higher than 50 GPa, while they agree with each other at relatively lower pressures.Physically speaking, 2p3s-model is better than 3s-model, since the effects of inner-core 2p states play important roles at extremely high pressures [5,6].Therefore it is a theoretical mystery why 2p3s-model gives higher melting temperature than those obtained by 3s model and experiment.
There are two possible reasons for this mystery.One of them is the effect of system size on the melting temperature obtained by ab initio MD simulations.In the previous ab initio MD simulations [2-6], the system sizes, i.e. the number of atoms in the simulation cell are at most 108 for fcc and 128 for bcc structures, respectively.So far no systematic study has been done for the effect of system size on melting temperature of Na.
The other possible origin is the method of determining melting temperatures.The melting temperatures are obtained in previous ab initio MD simulations [2][3][4][5][6] from the one-phase approach, i.e by heating the solid (bcc or fcc structures) until it melts, in which the melting temperature is determined from the drastic change of the timedependence of the mean-square-displacements (MSD), the 1.Melting curves of Na obtained from our ab initio MD simulations [5,6] compared with previous ab initio studies [3,4] and experimental results [1,7].The large circles and large triangles show our resutls obtained for 3s and 2p3s model, respectively, where black symbols show solid states and white ones show liquid states.The small triangles and squares show the results by Raty et al. [3] and Koči et al. [4], respectively.The experimental results are shown by small diamonds [1] and small circles [7].
structure factors and the radial distribution functions.It is pointed out recently by Bouchet et al. [8] using ab initio MD simulations that the melting temperature obtained by one-phase approach tend to be overestimated and the deviation from the true melting temperature is significantl large at extremely high pressures.It should be noted that the melting temperature obtained by one-phase approach is an upper limit of the true melting temperature.
In this paper, we report the system-size effects on melting of dense sodium at about 100 GPa.For this purpose we have performed ab initio MD simulations with the number of atoms at most 250 for bcc and 256 for fcc structures, respectively, based on 3s-valence-electron model and 2p3svalence-electron model.

Method of calculation
We carried out ab initio MD simulations [9][10][11] based on the density functional theory with the generalized gradient approximation (GGA) [12] for the exchange-correlation energy and on the projector augmented-wave (PAW) method [13,14] for the electron-ion interaction.In order to investigate systematically the effects of system size on the melting temperature at extremely high pressures, we carry out ab initio MD simulations for 3s and 2p3s models at the density of 3.7 g/cm 3 , which corresponds to the pressure of 100 GPa, for bcc structures with 54, 128 and 250 atoms/cell and for fcc structures with 32, 108 and 256 atoms/cell, for 4,000 time steps (7.74 ps).
We have determined the melting temperature as follows.With increasing temperature, (1)the time dependence of the mean-square-displacement (MSD) is qualitatively different in the solid and liquid states, i.e. below and above the melting temperature, respectively; the MSD oscillates around the constant value in the solid, while the MSD increases in proportion to time in the liquid, where the gradient of MSD is proportional to the diffusion constant, (2)the structure factor S (k) of solid has sharp spiky peaks, while S (k) of liquid show a smooth oscillating behaviour, (3)the radial distribution function g(r) of solid has a fin structure peculiar to crystal structures, while g(r) of liquid is a smooth oscillating curve.Using these criteria, we estimated the melting temperature, which exists between the highest temperature of solid state and the lowest temperature of liquid state.As mentioned in Sect.1, the melting temperature obtained in this way, i.e. one-phase approach, is an upper limit of the true melting temperature.

Results of calculation 3.1 3s model
In Fig. 2 we show the system-size dependence of the structure factors S (k) of Na for 3s model at 100 GPa in the (a) solid and (b) liquid states.The S (k) of the solid states show sharp spiky peaks peculiar to bcc and fcc structures.The S (k) of liquid states are calculated at the temperatures 50 K above the corresponding solid states.As is seen from these figures the melting temperature obtained based on 3s model depends strongly on the system size and decreases with increasing the number of atoms N in the simulation cell.For the fcc structure with N=256 we did not obtain the solid state, since there appears liquid state at the low temperature of 80 K, which is far below the room temperature.
In Fig. 3 we show the system-size dependence of the radial distribution functions g(r) of Na at 100 GPa in the (a) solid and (b) liquid states calculated at the temperatures corresponding to those shown in Fig. 2. As is seen from the figure the g(r) of solid states have some fin structures which are reminiscent of solid structure, while g(r) of liquid states show smooth oscillation.
The system-size dependence of the melting temperature obtained for 3s model is shown in Fig. 4. It is clearly seen that the melting temperature of 3s model decreases with increasing N and deviates from the experimental value substantially for N 250.

2p3s model
In Fig. 5 we show the system-size dependence of the S (k) of Na for 2p3s model at 100 GPa in the (a) solid and (b) liquid states.The S (k) of the solid states show sharp spiky peaks peculiar to bcc and fcc structures.The S (k) of liquid states are calculated at the temperatures 100 K above the corresponding solid states.As is seen from these figures the melting temperature obtained based on 2p3s model depends weakly on the system size and decreases very slowly with increasing N.
In Fig. 6 we show the system-size dependence of the g(r) of Na for 2p3s model at 100 GPa in the (a) solid and 01009-p.2(b) liquid states calculated at the temperatures corresponding to those shown in Fig. 5.As is seen from the figure the g(r) of solid states have some fin structures which are reminiscent of solid structure, while g(r) of liquid states show smooth oscillation.
In Fig. 7 the system-size dependence of the melting temperature obtained for 2p3s model is shown.The melting temperature of 2p3s model decreases very slowly with increasing N and tends to approach the experimental value.

Conclusion
Applying ab initio MD simulations to solid and liquid Na at 100 GPa based on 3s and 2p3s model with at most 256 atoms/cell, we have shown that the inner-core 2p states play important roles at extremely high pressures and that the melting temperature estimated by one-phase approach depends on the system size; in particular, the melting temperature for 3s model decreases with increasing system size and deviates from the experimental value, while the melting temperature for the 2p3s model approaches to the experimental value slowly.This result suggests that N 100 is not large enough to obtain the melting temperature.From these results we conclude that we have to carry out large-scale ab initio MD simulations with the innercore 2p states to obtain the correct melting temperature at extremely high pressures.Therefore we need further study for the melting of Na.

Fig. 2 .
Fig. 2. The system-size dependence of the structure factors S (k) of Na for 3s model at 100 GPa are shown in the (a) solid and (b) liquid states.

Fig. 3 .
Fig. 3.The system-size dependence of the radial distribution functions g(r) of Na for 3s model at 100 GPa are shown in the (a) solid and (b) liquid states.

Fig. 4 .
Fig. 4. The system-size dependence of the melting temperature for 3s model at 100 GPa.

Fig. 5 .
Fig. 5.The system-size dependence of the S (k) of Na for 2p3s model at 100 GPa are shown in the (a) solid and (b) liquid states.

Fig. 6 .
Fig. 6.The system-size dependence of the g(r) of Na for 2p3s model at 100 GPa are shown in the (a) solid and (b) liquid states.

Fig. 7 .
Fig. 7.The system-size dependence of the melting temperature for 2p3s model at 100 GPa.