Structural studies of liquid alkaline-earth metals: A molecular dynamics approach

Structural studies of liquid alkaline-earth metals -A molecular dynamics approach

J. K. Baria^I; A. R. Janib^II

^IV. P. & R. P. T. P. Science College, Vallabh Vidyanagar – 388 120, Gujarat, India

^IIDepartment of Physics, Sardar Patel University, Vallabh Vidyanagar – 388 120, Gujarat, India

ABSTRACT

In the present research article we have studied various properties like binding energy, the pair distribution function g(r), the structure factor S(q), specific heat at constant volume, and coordination number of alkaline-earth metals (Be, Mg, Ca, Sr, and Ba) near melting point temperature using molecular dynamics (MD) simulation technique with a pseudopotential proposed by us. Good agreement with the experiment is achieved for the binding energy, pair distribution function, and structure factor and these results compare favorably with the results obtained by other such calculations, showing the transferability of the pseudopotential used from solid to liquid environment in the case of alkaline-earth metals.

Keywords: Molecular dynamics; Pair distribution function; Structure factor; Binding energy

1. INTRODUCTION

The alkaline-earth metals (IIA column) can be considered as simple metals in regard of the electronic aspect, they are certainly the less studied elements among the metals. From the experimental point of view, only a few properties have been determined to date. Indeed, we can mention the electrical resistivity [1], the absolute thermopower [2], the static structure [3], the sound velocity [4], and the density [5]. This is due to the high chemical reactivity and to the gas adsorption ability [6] that further increase with temperature. These difficulties have induced a low interest for technological purposes and a disinterest for theoretical work on alkaline-earth metals, while alkali and polyvalent ones were often preferred. However, alkaline-earth metals have a central place. On one hand, they can be considered as simple and free-electron-like, since they have only s-like conduction electrons. Beryllium and magnesium, essentially free-electron-like, are as simple as sodium [7]. On the other hand, in the case of Ca, Sr, and Ba the existence of an empty d-band above the Fermi level, which get nearer as one goes down the (IIA column) [8] has a great influence on the electronic properties. This explains the high electrical resistivity of liquid barium [1,9].

Recent works [8,10-12] have been devoted to this class of elements to see trends through the Periodic Table and to study the transition from alkali to more complex metals in relation to the influence of the electron gas. Alemany et al. [12] also computed the velocity auto-correlation function VACF, its memory function and the self-diffusion coefficient.

In this paper, we study the properties like binding energy, the pair distribution function g(r), the structure factor S(q), specific heat at constant volume, and coordination number using molecular dynamics (MD) simulation technique with a pseudopotential of Baria and Jani [13-15] for alkaline earth metals, namely beryllium, magnesium, calcium, strontium, and barium. We present the results of the static properties i.e., Binding energy, the structure factor, and the pair distribution function and compared with both experimental and theoretical previous works, as well as the VACF and its spectral density. Because of the lack of experimental data for VACF and its spectral density, we are only able to compare our results with the recent ones of Alemany et al. [12]. Our calculations have been performed with the molecular dynamics simulation technique with the pseudopotential of Baria and Jani [13-15] in the second-order perturbation theory.

2. Theory

Simple metals are usually depicted as an assembly of ions of well-defined electric charge immersed in the bath of the conduction electrons, the global system being electrically neutral. Since the number of conduction electrons per atom has an integer value (Z=2 for alkaline-earth metals), the direct ion-ion repulsion, essentially Coulombic, is easy to describe contrary to the transition metals. In contrast, the description of the effective ion-ion interaction is drastically affected by the electron-electron and the electron-ion interactions.

The electron-electron interaction relies on the electronic charge density, which takes also a well-defined value for simple metals. Its description is performed in terms of screening with two functions, the Lindhard-Hartree dielectric function ε_H (q) and the local-field correction G(q), which account for electrostatic and exchange-correlation effects, respectively. Wax and Bretonnet [16] have shown that the important characteristic of G(q) at least for the static structure factor of liquid metals is the low-q limit parameter depending on the correlation energy of the electron gas. The exchange-correlation effect is now a days well-predicted by Monte Carlo techniques, and the Ichimaru and Utsumi [17] IU expression of the local-field correction used in this work is certainly one of the most reliable. Since alkaline-earth metals are simple like metals, the main difficulty in the calculation of the effective ion-ion interaction lies in the description of the electron-ion interaction. Thus, the effective pair potential can be obtained from the perturbation theory applied to second order to the determination of the configurational energy of the material,

F_N (q) being the normalized energy-wave number characteristic,

Where, , being the atomic volume, the electron-gas dielectric function ε(q), and the local-?eld correction G(q) that takes into account the exchange-correlation effects [17] between conduction electrons and W_B (q) is the pseudopotential proposed by Baria and Jani [13-15], (in Ryd. Unit) is given by,

where r_c is the parameter of the potential and e the base of natural logarithm. We have calculated the parameter of the potential rc using standard zero pressure technique. The corresponding parameters are compiled in Table I. One further variable aspect of the pseudopotential formalism is the choice of the exchange-correlation function, G(q). We employ the expression exempt from free parameters developed by Ichimaru and Utsumi [17].

In the framework of the second-order perturbation theory of the pseudopotential, the binding energy is given by the energy of an atom as,

The configurational energy E_con is obtainable directly by molecular dynamics or from its definition in terms of pair distribution function g(r),

The volume-dependent contribution to the binding energy, E_vol, is given according to the prescription of Hasegawa and Watabe [18] by

In this expression, E_eg represents the ground state energy of the electron gas, for which the following Nozieres-Pines [19] interpolation formula is used (within atomic units):

r_s being the radius of the electronic sphere defined by The first term of Eq. is the kinetic energy of the free electron gas, the second term is the attractive exchange energy due to the parallel-spin electrons separated by Pauli's exclusion principle, and the third term is the correlation energy that gives an additional lowering in energy.

In the right-hand side of Eq. , the second term, usually noted Beg, corresponds to a rearrangement of various energetic contributions for the zero-wave vector and the third term noted as Φ, represents the self-energy between an ion and its surrounding cloud of charge. In this context, the local pseudopotential assumed for the ion-electron interaction and the local-field correction G(q) considered for the exchange and correlation effects within the electron gas are crucial for obtaining the volume-dependent energy.

3. MOLECULAR DYNAMICS SIMULATION

The molecular dynamics (MD) simulations are performed in the NVE microcanonical ensemble. The system is a cubic box containing N atoms with periodic boundary conditions and the box sizes are fitted to the desired density, During the thermalization stage, the system relaxes and the velocities are rescaled each 50 time steps to the expected temperature. Once the system is thermalized, the production stage is launched and positions and velocities of the particles are taped each 5 (for dynamical properties) or 10 (for static ones) time steps. During this stage, the temperature is no more renormalized. Because of the statistical fluctuations, for which standard deviations of the temperature are a few tenths of degrees, the effective temperature during the production differs from the expected one.

In this paper, we always indicate the effective value of the temperature at which the results are obtained, since, as it will be seen later, the self-diffusion coefficient is highly sensitive to the temperature. It is worth to mention that the isothermal MD simulation, which requires the renormalization of the velocities, is prescribed in observing the evolution of a tagged particle in order to compute the VACF. Depending on the physical properties under study, two different sizes for the simulation box is used; 256 particles for dynamical properties and 2048 particles for static properties. The bigger is the system, the faster is the thermalization and the smaller are the fluctuations of the temperature. Knowing a great number of successive configurations, the pair distribution function g(r) can be determined from the relationship,

Where, C is the number of configurations taped and n(r) is the number of atoms counted at a distance between r and r+δ r from another atom taken as origin. As we said, for computing the static structure the box contains N=2048 atoms in order to ensure a large enough r extension of g(r). So, the structure factor S(q) is obtained straightforwardly by Fourier transform as,

Another function 4 π ρ r² g (r) obtained from g (r) is used in the discussion of the structure of non-crystalline systems. This has been called the radial distribution function (RDF). This function corresponds to the number of atoms in the spherical shell between r and r+dr, the coordination number, and is obtained from relation [20],

Where, r₀ is the left hand edge of the first peak and r_m corresponds to the first minimum on the right hand side of the first peak in RDF.

4. RESULTS AND DISCUSSION

To determine the binding energy by Eq. , we first calculate the configurational internal energy E_con by MD simulation. We also take advantage of the fluctuations in this energy during the course of the calculation to evaluate the specific heat at constant volume using the standard expression [21],

〈E_kin〉 being the average kinetic energy and δE_con, the fluctuation in the configurational energy, Table 1 presents potential parameter rc, and E_con with the other contribution to the volume-dependent energy, as well as the experimental binding energy [22]. Note that the main contribution Φ, to the binding energy comes from the electrostatic interaction between an ion and its own screening cloud of electrons. Our present values of binding energies are in good agreement with the experimental values within percentage deviation of 1.11% to 4.24%.

As a test of transferability of the pseudopotential under study, we can compare the binding energy with experiment since it corresponds to the sum of the ionization energies of the valence electrons, E_ion, and the cohesive energy, E_coh. Indeed, we recall that the conventional binding energy can be written [23] as

Our results for the pair distribution function g(r) are presented in Fig. 1 for Be, Mg, Ca, Sr, and Ba near their melting point temperature in comparison to the experiments [20]. To the best of our knowledge, no experimental results have ever been published for Be. For Mg, Ca, Sr, and Ba there is a slight overestimation of the height of the first peak with the experimental values of Waseda [20] except these, our present findings are excellent. The corresponding static structure factors S(q) are displayed in Fig. 2. Though our calculations are performed with 2048 particles, we can see that the low-q behavior of S(q) is poorly described because of the smallness of the simulation box. The mesh in the q space being governed by the extension of g(r), it may happen that the curves of S(q) are not smooth in some crucial ranges, like the position of the main peak, due to the lack of points. Thus, for Mg, Ca, and Ba, the height of the first peak of S(q) seems to be slightly underestimated. Except these small differences, the agreement between our calculations and experiments [20] is overall very satisfactory. The first and second peak position and its relative magnitude of g (r), and coordination number are shown in table 2 while first and second peak position and its relative magnitude of S(q) are shown in table 3. Presently calculated values of first and second peak position of g (r), coordination number and S(q) are excellently agrees with the experimental findings of Waseda [20] for Mg, Ca, Sr and Ba while for Be, no experimental data are available for comparison.

5. CONCLUSION

It is rather difficult to estimate the error for the dynamical properties in computer simulations because many factors influence the results of the calculations. These are the relatively small number of particles, the periodic boundary conditions, the number of configurations used in calculating various averages and the shortcomings of the pair potential, and hence the pseudopotential.

The properties of crystalline and liquid phases are found to be adequately modeled by the pseudopotential of Baria and Jani [13-15]. Since the self-diffusion coefficient is very sensitive to the temperature and it is doubled at least at each 400K, it might be a very accurate test of the potential. The static structure does not appear to be very sensitive to it as we can see from the results of g(r) and S(q), both being in very good agreement with experiment. So, for this kind of metals, we have to consider a much more sensitive property, and the self-diffusion coefficient could be this one. Unfortunately, the lack of experimental value prevents definitive conclusions.

To conclude, we have performed molecular dynamics study of liquid alkaline-earth metals with pseudopotential of Baria and Jani [13-15] in order to check the transferability of the model. Considering that this pseudopotential is exempt from adjustable parameters, the good quality of our results of the structure factor is a strong argument in favor of the transferability of the pseudopotential from the solid state to the liquid state for alkaline-earth metals as far as the pair potential and the ionic structure is concerned. The study of the dynamic properties and above the entire self-diffusion coefficient points out the relevance of this latter quantity to deal with the transferability of the pseudopotential. The existence of experimental data would be helpful.

Acknowledgement

One of the authors (JKB) is thankful to the University Grants Commission (UGC-WRO, Pune) for providing financial support under minor research project grant (F. No.: 47-719/08).

(Received on 14 January, 2010)

• [1] J.B. van Zytveld, J.E. Enderby, and E.W. Collings, J. Phys. F: Met. Phys. 2 73 (1972).
• [2] J.B. van Zytveld, J.E. Enderby, and E.W. Collings, J. Phys. F: Met. Phys. 3 1819 (1973).
• [3] Y. Waseda, K. Yokoyama, and K. Suzuki, Philos. Mag. 30 1195 (1974).
• [4] S.P. McAlister, E.D. Crozier, and J.F. Cochran, Can. J. Phys. 52 1847 (1974).
• [5] S. Hiemstra, D. Prins, G. Gabrielse, and J.B. van Zytveld, Phys. Chem. Liq. 6 271 (1977).
• [6] J.G. Cook and M.J. Laubitz, Can. J. Phys. 54 928 (1976).
• [7] G.A. de Wijs, G. Pastore, A. Selloni, and W. van der Lugt, Phys. Rev. Lett. 75 4480 (1995).
• [8] W. Jank and J. Hafner, Phys. Rev. B 42 6926 (1990).
• [9] V.K. Ratti and R. Evans, J. Phys. F: Met. Phys. 3 L238 (1973).
• [10] L.E. Gonzalez, A. Meyer, M.P. Iniguez, D.J. Gonzalez, and M. Silbert, Phys. Rev. E 47 4120 (1993).
• [11] Hong Seok Kang, Phys. Rev. B 60 6362 (1999).
• [12] M.M.G. Alemany, J. Casas, C. Rey, L.E. Gonzalez, and L.J. Cal-lego, Phys. Rev. E 56 6818 (1997).
• [13] J. K. Baria, A. R. Jani, Physica B 328 317 (2003).
• [14] J. K. Baria, A. R. Jani, Pramana J. Phys. 60 1235 (2003).
• [15] J. K. Baria, A. R. Jani, Physica B 404 2401 (2009).
• [16] J.-F. Wax and J.-L. Bretonnet, J. Non-Cryst. Solids 250-252 30 (1999).
• [17] S. Ichimaru and K. Utsumi, Phys. Rev. B 24 7385 (1981).
• [18] M. Hasegawa and M. Watabe, J. Phys. Soc. Jpn. 32 14 (1972).
• [19] D. Pines and P. Nozie'res, Quantum Liquids Benjamin, (New York, 1966).
• [20] Y. Waseda, The Structure of Non Crystalline Materials (McGraw-Hill, New York, (1980).
• [21] J.M. Haile, Molecular Dynamics Simulation: Elementary Methods (Wiley, New York, 1992).
• [22] K.A. Gschneider, Solid State Phys. 16 275 (1964).
• [23] M. Shimoji, Liquid Metals (Academic, London, 1977).

Publication Dates

• Publication in this collection
  23 June 2010
• Date of issue
  June 2010

History

• Received
  14 Jan 2010