- Email: [email protected]
The structural and dynamical properties of ZnCl2 melt—a molecular dynamics simulation study
Journal of Molecular Liquids 115 (2004) 81 – 88 www.elsevier.com/locate/molliq
The structural and dynamical properties of ZnCl2 melt—a molecular dynamics simulation study Shiping Huang a,*, F. Yoshida b, Wenchuan Wang a a
Laboratory of Molecular and Materials Simulation, College of Chemical Engineering, Beijing University of Chemical Technology, P.O. Box 100, Beijing 100029, PR China b Department of Physics, Shiga University of Medical Science, Seta Ohtsu, Shiga 520-2192, Japan Received 8 January 2004; accepted 25 February 2004 Available online 26 June 2004
Abstract The shell-model computer simulation is performed to study polarization effects of ions in structural properties of molten ZnCl2. Experimental results of partial structure factors as well as self-diffusion coefficients are well reproduced by a suitable choice of potential parameters. By using the same parameters, time-correlation functions of charge and mass current are evaluated. Discussion is given about broad bands in the infrared absorption spectra, and longitudinal (or transverse) and acoustic (or optic) modes in excitation spectra, with varying core charges of the Cl ion. D 2004 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics simulation; ZnCl2 melt; Polarization effects
1. Introduction Molten ZnCl2 has some uncommon physical properties, such as high viscosity, low conductivity and high tendency to form a glass. For this reason, the properties of ZnCl2 in melt and glass states have attracted special attention, with a number of theoretical and experimental works [1,11]. Biggin and Enderby [2] have used the technique of the 35Cl/37Cl isotopic substitution in neutron diffraction experiments to determine three partial distribution functions, whereby assisted the interpretation for characteristic structures of molten ZnCl2. Wong and Lytle [3] have used the extended X-ray absorption fine structure to study both glass and liquid states. Recently, Allen et al. [4] have revealed by using the technique of pulsed neutron diffractions that the basic ZnCl4 structural units are extremely well defined, and the coordination number decreases from 3.9 to 3.7 with temperature on heating from 330 to 660 jC. The properties of molten ZnCl2 in the temperature range of 623 –853 K are also investigated by using high energy-photon diffractions. It is found that * Corresponding author. Tel.: +86-010-64444905; fax: +86-01064427616. E-mail address: [email protected] (S. Huang). 0167-7322/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2004.02.045
the first peak position of Zn – Cl in partial distribution functions does not change with increasing temperature, but the first shell Zn – Zn distance shows some decrease [5]. All these experimental facts lead us to the microscopic understanding that the small Zn ion occupies tetrahedral coordination sites formed by the close Cl ions, and the tetrahedral coordination remains well-defined above the melting point. Another characteristic feature of molten ZnCl2 is a first sharp diffraction peak (FSDP) at scattering ˚ 1, which is interpreted as intermediate vector about 1 A range order (IRO) in the melt. Molecular dynamics simulation has been widely used as complementary but powerful method to study the structure and dynamic properties of liquids and glasses. For molten ZnCl2, the interaction between two ions is complicated due to the polarization of ions or covalent bond to some degree. The development of a potential model consistent with experimental results is therefore a challenging task. In the first approach to ZnCl2 melt, Woodcock et al. [6] have constructed a potential the Born-Mayer form (called the WAC potential model). Their result for the distance of Zn – Zn is however too large compared to experimental results. The KDR potential model introduced by Kumba et al. [7] including a dispersion term between zinc ions seems to generate the Zn – Zn correlation functions, but there is no
82
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
physical meaning for the zinc ion radius. The authors have used the compressed potential model for simulating molten ZnCl2, but could not predict accurately the Zn –Zn correlation function. Pusztai and McGreevy [8] have applied the reverse Monte Carlo (RMC) simulation method to analyze the structure of the ZnCl2 molten and glass state. Their results show that the structural units (ZnCl4 tetrahedral) are connected at corners. Bassen et al. [9] have used the Monte Carlo and the RMC method to perform a simulation of molten zinc chloride on the basis of the X-diffraction experiments at 623 and 873 K. They have also modified the Coulomb interaction with a dielectric screening function, and found that the simulation using their potential model is capable of reproducing experimental results for correlation functions, the total structure factor and prepeaks. Salmon [10] has systematically analyzed partial structure factors of the molten and glassy MX2 binary system, and found that the first diffraction peaks or pre-peaks have mainly their origin in the number –number partial structure factor SNN(k) and number– concentration partial structure factor SNC(k). Recently, Ribeiro et al. [11] and Gray-Weale et al. [12] have developed a dynamically polarizable model which incorporates a polarization force by combining the Cl dipole polarization with the WAC potential for MX2 melts. They have shown reasonable consistency of the model with experiments for static structures, and further discussed dynamical features in excitation or absorption spectra. In this paper, we use the shell-model potential in the molecular dynamics simulation to study the structural properties of ZnCl2 molten salt. This method is found useful to take into account the polarizability of ions compared to the rigid ion model. The present study has two objectives. The first is to see how our approach using the polarizable potential model is useful by comparing with experimental results of static structures as well as self-diffusion constants. The second, more interesting, is to get an insight into the dynamical structure of this network liquid. It is reported by the instantaneous normal mode analysis that in dynamics the system could not be simply regarded neither as being composed of weakly coupled tetrahedral units with intramolecular type modes, nor of strongly coupled individual ions allowing extended phonons [11]. We discuss polarization effects in dynamical properties by calculating spectra of mass and charge current correlation functions.
is complicated because of its directional character in some degree. Therefore, the polarizable effects or higher body forces must be taken into account. The dynamical shell model recently devised by Fincham and Mitchell [13] incorporates polarizability into the molecular dynamics simulation. In the shell model, the ion is divided into two components, a core with the ion mass and a shell with no mass, which are connected by a harmonic spring. The two parts carry different charges, and the electrostatic interaction between the core and the shell does not include the self interaction (the core and shell of the same ion) [14]. In the simulation of molten zinc chloride, the chloride ion is assumed to have a core and a shell connected by a harmonic spring with the spring constant k. They carry the different charges, but the sum of their charges is equal to 1. The interaction potential /ij(r) between ions i and j is given by ! qi qj Cij r þ Aij exp /ij ðrÞ ¼ 6 qij r r It consists of the long-range Coulomb part and the shortrange repulsion part, with r as the inter-ionic distance, qij the effective hardness and Cij the dispersion parameter. The parameters adopted in the simulation are listed in Table 1: core charges are in units of the elementary charge e, and the value for the Cl ion is taken as 0.984 e from Ref. [15] fitting the crystal properties of ZnCl2. The simulation system consists of 1500 ions (500 Zn and 1000 Cl ions), and the initial position of ions is chosen as CaF2 type in a cubic box. The Ewald summation is used for calculating the electrostatic interaction. The time step is taken as 0.001 ps in the simulation. The temperature is initially controlled to be 2000 K up to 10 000 steps, and it is then reduced to the set temperature. In order to get the density at the set temperature, the simulated system is taken under constant pressure (0.11 MPa) in conjunction with temperature control by using a loose-coupling method, and then the average density is calculated for a constant temperature running. The calculated density at 600 and 873 K is 2.601 and 2.50 g cm 3, respectively. The experimental values at temperatures corresponding to them are 2.52 and 2.40 g cm 3 [6]. It takes 20 ps for the simulation system to achieve an equilibrium state (constant temperature, volume and number of particle), where the energy fluctuation is less than Table 1 Potential parameters for the molecular dynamics simulation [15]
2. Computational detail In the molecular dynamics simulation, each ion is treated as a point particle and the motion of an ion is based on Newton’s equation. The thermodynamic quantities and transport coefficients are calculated by the time average over the ensemble of particles. For the molten and glass state of zinc chloride, the interaction between Zn and Cl ions
Zn Cl
Zn – Zn Zn – Cl Cl – Cl
Q (core) (e)
q (shell) (e)
˚ 2) Spring constant k (eV/A
2.000 0.984
0.000 1.984
17.250
Aij (eV)
˚) U (A
˚ 6) Cij (eV/ A
0.0000 9704.8900 3296.5700
0.0000 0.2320 0.3289
0.0000 0.0000 107.2000
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
83
0.2%. The structure and dynamical quantities are calculated by taking the average over more than 25 ps. All the simulations are performed by using DL_POLY code [16].
3. Results and discussion The static structure of molten systems is usually described by the partial radial distribution function (RDF). It represents the probability of finding another particle from a particle at origin as a function of the radial distance. Fig. 1 shows the partial RDF at T = 1000, 873 and 600 K. Table 2 presents the summary of results for the partial RDF, in comparison with the experiments. The results of the molecular dynamics simulation are found to be in good agreement with experiments and the RMC simulation. The ratio rClCl/ rZnCl for first peak positions from the partial RDF is about 1.61, close to the value of 1.63 for the regular tetrahedral [2,9,17]. The partial structure factors sij(k) for i and j species (i, j = Zn, Cl) are obtained by the Fourier transformation of the partial RDF gij(r) [18] Z l 4p sij ðkÞ ¼ dij þ pffiffiffiffiffiffiffiffi ½gij ðrÞ 1sinðkrÞrdr k ni nj 0 with ni being the number density of i species. The total structure factor is then calculated from sij(k) according to the formula [4] FðkÞ ¼ c2i b2i ½sii ðkÞ 1 þ c2j b2j ½sjj ðkÞ 1 þ 2ci cj bi bj ½sij ðkÞ 1 with bi being the scattering length and ci the concentration of particle i. It is experimentally determined by neutron or X-ray diffractions. The calculated results for sij(k) and F(k) are shown in Figs. 2 and 3, respectively. It is noted that there exists the lower limit of k = 2p/L owing to the finite size L of the simulated system in calculating partial structure factors ˚ 1). We can observe clearly the pre-peak at (kmin = 0.20 A 1 ˚ ), with good agreement between our simulated k f 1(A results and neutron diffraction experiments in the total structure factor. These results imply a validity of our potential model incorporating polarizability to investigate the structure of ZnCl2 molten salt. In the present potential model, the electric field for unlike charges in an ion (Cl ion) due to surrounding ions tends to make the ion to separate the shell from the core, and it gives rise to a polarization dipole. The polarizability of the particle in the equilibrium depends on the core charge. We have therefore changed the core charge from the set value of 0.984 e in Table 1 down to 0 to see polarization effects. Fig. 4 shows the calculated result for the dependence of the partial RDF gZn – Zn on the core charge. It is found that both the position and the height of the first peak remarkably decrease as the core charge (and thus polarizability)
Fig. 1. Partial radial distribution function : (a) 1000 K; (b) 873 K; (c) 600 K; Line: gZn – Cl(r) Dot: gCl – Cl(r) ; Dash: gZn – Zn(r).
increases. From micro-structural viewpoint, these facts mainly reflect the appearance of the structure of the [ZnCl4] tetrahedral, whose units are connected by the chlorine bridge or edge [9]. The Coulomb repulsion between doubly charged cations appears, in effect, to be significantly reduced by the polarization effect. The connection of the
84
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
Table 2 Summary of static results Temperature (K)
Method
600 873
MD MD Neutron2 RMC9 RMC9 X-ray17
623 873
First peak position (nm)
Coordination
Zn – Zn
Zn – Cl
Cl – Cl
Zn – Zn
Zn – Cl
Cl – Cl
0.39 0.39 0.38 0.385 0.385 0.366
0.23 0.23 0.229 0.220 0.225 0.229
0.37 0.37 0.371 0.355 0.360 0.385
4.271 4.359 4.7
40.080 4.106 4.3
10.377 10.124 8.6
[ZnCl4] tetrahedral, or the tendency for forming a network structure, would get stronger with increasing polarizability. The diffusion process of ions can be investigated through its mean-square displacement as a function of time. The self-diffusion constant is evaluated by the formula < Arðt þ sÞ rðsÞA2 > t!l 6t where t is the time interval, and the angular bracket stands for the ensemble average. It is evaluated by using fitting procedures during 6 ps after 20 ps required to be achieved the equilibrium in the molecular dynamics simulation. The simulated result obtained by using the values of the parameters in Table 1 (the core charge is 0.984 e) is shown as a function of temperature in Fig. 5. It is observed that the value of the self-diffusion coefficient for Zn is always smaller than that for Cl in this temperature range, and the behavior appears similar for both Zn and Cl in the dependence on temperature. This is because the formation of the [ZnCl4] unit or units connected to form more complex structure restricts the mobility of ions. The experimental results for the self-diffusion coefficient are 0.010 and 0.013 (unit is 1 10 5 cm2/s) for Zn and Cl ions at 600 K [19], respectively. The corresponding values of our simulation are 0.04241 and 0.08069 (unit is 1 10 5 cm2/s) for Zn and Cl ions, respectively, which corresponds with the same magnitude of experimental value. D ¼ lim
Fig. 2. Partial structure factors at 600 K. Solid line: sZn – Cl(k), Dash-Dot: sZn – Zn(k), Dot: sCl – Cl(k).
4
4
12
The diffusive motion of ions is not independent of polarization effects, which is responsible for the formation of the tetrahedral unit. Fig. 6 shows the dependence of diffusion coefficients on the core charge at T = 873 K. It is found that diffusion coefficients of both ions show similar behavior of increase with the core charge. It depends on temperature, but the result may be interpreted in the way that the increase of the core charge facilitates the formation of many tetrahedral units, but these units would be rather loosely connected with each other. As a result, ions in the molecular unit appear to diffuse more easily, compared to the case of the individual motion at zero core charge. The electrical conductivity is related to the correlation function of the total charge current correlation function J(t) via its Laplace transform a(x), defined by the following equations [19] X jðtÞ ¼ qa via i;a
< jðt þ sÞjðsÞ > < jðsÞjðsÞ > Z l aðxÞ ¼ expðixtÞJ ðtÞ
J ðtÞ ¼
0
where j(t) is the total charge current, qa is the charge of a species and via is the velocity of particle i of a species.
Fig. 3. Total structure factor at 600 K, with the experimental data from Ref. [4].
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
Fig. 4. Partial radial distribution function of Zn – Zn as a function of core charge of Cl ion at 873 K.
85
Fig. 6. Self-diffusion coefficient as a function of core charge at 873 K.
The results of the molecular dynamics simulation for a(x) at representative temperatures ranging from 600 to 1200 K are shown in Fig. 7. The total charge current correlation function reflects a collective property of the system, and according to the fluctuation theory, the infrared absorption can be obtained from the Fourier transformation of the dipole correlation function, which is also related with a(x). There are two broad bands in Fig. 7; one band is around 3.141 THz ( f 104.77 cm 1), and the other is around 10.050 THz ( f 335.23 cm 1). As the temperature increases, the high frequency band makes a shift towards the low frequency region with broader width. Angell and Wong [20] have measured the infrared spectra of ZnCl2 molten salt, and observed two broad bands at 100– 115 and at 255 cm 1. The position of the high frequency band in our results is found to be rather higher in frequency compared with the experimental result. A similar situation was reported for the dynamically polarizable model in comparison with experiments [12]. The peak positions in the Fourier transform of the total current correlation function represent the transverse optic
frequency. The mass and charge currents are defined by the equations [21]
Fig. 5. Self-diffusion constant as a function of temperature.
Fig. 7. Spectrum of the total charge current correlation function.
Jm ðk; tÞ ¼
N X
expðik ri ðtÞÞmi vi ðtÞ
i¼1
Jq ðk; tÞ ¼
N X
expðik ri ðtÞÞqi vi ðtÞ
i¼1
Here, mi and qi are the mass and charge of ion i, respectively. The transverse and longitudinal correlation functions for the mass and charge currents are given by T CXX ðk; tÞ ¼ ð1=2N Þ < k JX ðk; t þ sÞ k JX ðk; sÞ > L ðk; tÞ ¼ ð1=N Þ < k JX ðk; t þ sÞk JX ðk; sÞ > CXX
with X = m for the mass current and X = q for the charge current, where the symbol denotes the vector product. When the spectra or Fourier transforms of the transverse and longitudinal mass correlation function exhibit sharp peaks, the peak position is identified as the transverse and longi-
86
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
˚ 1), Dot: k = 1.0 (A ˚ 1). (a) TA, Fig. 8. Spectrum of transverse and longitudinal mass and charge correlation functions at T = 873 K. Solid line: k = 0.45 (A (b) LA, (c) TO, (d) LO.
tudinal acoustic (TA and LA) frequency, respectively. Similarly, peaks in the spectrum for the charge current correlation function imply the transverse and longitudinal optic mode’s frequencies (TO and LO). We have calculated these spectra, and the results at ˚ 1) and k = 1.0 (A ˚ 1), T = 873 K are shown at k = 0.45 (A in Fig. 8. It is noted that the latter is the wave number at prepeaks indicating the specific structure of the system. In the spectrum of acoustic modes, there are many fine structures for the LA mode, compared to the TA mode. It may correspond to the result by the INM that normal modes are not interpreted by single plane waves [12], or the density fluctuation with a definite wave number is not a welldefined variable. Both the TA and LA modes make a shift to higher frequency as the wave number k increases, which implies that these modes have a positive dispersion. As for optic modes representing the relative motion of ions, the peaks of TO and LO modes correspond to the two peaks in Fig. 7. It is found that these modes tend to slightly show a negative dispersion. These are in accordance with the results for similar systems [12]. We have further investigated the polarization effects on the excitation spectra in Fig. 9. The calculated results at T = 873 K for the transverse (or longitudinal) mass (or ˚ 1), charge) correlation function are shown at k = 0.28 (A
with representative values of the Cl ion core charge, i.e., 0, 0.5 e and 0.984 e. The complicated dependence on the core charge is seen for these spectra. For the acoustic spectra, the position of the main peak for 0.984 e seems to make a shift to lower frequency, compared to the result for 0.5 e. It implies an effective softening of the coupling constant for individual ions, presumably averaged over tight coupling in the tetrahedral unit, and relatively weak coupling between the units. For the optic spectra, it is found by comparing the cases of 0.984 e and 0.5 e that the two large peaks tend to separate, indicating a kind of the level repulsion, the results for the zero core charge show different behavior for the acoustic and optic spectra. It may be presumably correlated to the rapid change of the diffusion coefficient accompanied by the increase of the core charge from 0 to 0.5 e, as seen in Fig. 6.
4. Conclusion We have discussed the polarization effects in the static and dynamical properties of molten ZnCl2 by performing the shell-model computer simulation. It is found that experimental results of the total structure factor and diffusion constants are fairly well reproduced by appropri-
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
87
˚ 1). The Fig. 9. Dependence on the core charge for spectra of transverse and longitudinal mass and charge correlation functions with T = 783 K, at k = 0.28 (A Cl ion core charge is: Dash dot: 0.00 e; Dot: 0.500 e; Solid line: 0.984 e.
ately choosing parameters in the effective interionic potential and core charge for the ion (only for Cl ion). In addition, by investigating the dependence of the partial structure factor of Zn –Zn on the core charge, we have explicitly shown a significant contribution of the polarizability of ions: the ratio of the first peak positions for partial structure factors is found to be consistent with the formation of the tetrahedral unit [ZnCl4]. We have moreover found that the self-diffusion coefficients of Zn and Cl show the similar dependence on the core charge, which also might indicate the molecular structure mentioned above in the network liquid. As for dynamical properties, we have calculated the spectrum of the total mass and charge current correlation function. It is found by adopting a few values for the core charge that excitation spectrum for acoustic and optic modes are qualitatively understandable by the polarization effects. On the whole, our results are in accordance with other studies such as the instantaneous mode analysis. In comparison with observed infrared absorption spectra, however, some discrepancies are found for the two broad bands. This and similar results reported in the dynamically polarizable model suggest partly an inadequacy of polarization models in treating the bonding.
The significant dependence of the diffusion coefficient on the core charge of the Cl ion shown in Fig. 6 indicates that polarization effects coupled with temperature would improve the situation. For this and a more detailed understanding of the system, it is interesting to make an ab initio calculation for the structure, which will be done in the future.
Acknowledgements This work is supported by the National Basic Research Program of China, Grant No. 2003CB615807, and the Nation Natural Sciences Foundation of China (No. 20236010).
References [1] [2] [3] [4] [5]
M. Rover, M.P. Tosi, Rep. Prog. Phys. 49 (1986) 10001. S. Biggin, J.E. Enderby, J. Phys. C 14 (1981) 3129. J. Wong, F.W. Lytle, J. Non-Cryst. Solids 37 (1980) 273. D.A. Allen, R.A. Howe, N.D. Wood, J. Chem. Phys. 94 (1991) 5071. J. Neuefeind, K. Todheide, A. Lemke, H. Bertagnolli, J. Non-Cryst. Solids 224 (1998) 205.
88
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
[6] L.V. Woodcock, C.A. Angell, P. Cheeseman, J. Chem. Phys. 4 (1976) 1565. [7] P.N. Kumta, P.A. Deymier, S.H. Risbud, Physica. B + C 153 (1986) 85. [8] L. Pusztai, R.L. McGreevy, J. Non-Cryst. Solids 117/118 (1990) 627. [9] A. Bassen, A. Lemke, H. Bertagnolli, Phys. Chem., Chem. Phys. 2 (2000) 1445. [10] P. Salmon, Proc. R. Soc. Lond., A 437 (1992) 591. [11] M.C.C. Ribeiro, M. Wilson, P.A. Madden, J. Chem. Phys. 109 (1998) 9859. [12] A. Gray-Weale, P. Madden, M. Wilson, J. Chem. Phys. 113 (2000) 6782.
[13] [14] [15] [16] [17] [18]
D. Fincham, P.J. Mitchell, J. Phys., Condens. Matter 5 (1993) 1031. P.J.D. Linda, M.J. Gillan, J. Phys., Condens. Matter 5 (1993) 1019. D.J. Binks, PhD thesis, University of Surrey, 1994. W. Smith, T.R. Forester, The DL_POLY user manual (1999) 55. R. Triolo, A.H. Narten, J. Chem. Phys. 14 (1981) 703. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publication, New York, 1996. [19] J.P. Hansen, R. McDonald, Theory of Simple Liquids, 2nd ed., Academic Press, London, 1986. [20] C.A. Angell, J. Wong, J. Chem. Phys. 53 (1970) 2053. [21] B.J. Berne, R. Pecora, Dynamic Light Scattering, Wiley, New York, 1976, Sec 15.3.
The structural and dynamical properties of ZnCl2 melt—a molecular dynamics simulation study Shiping Huang a,*, F. Yoshida b, Wenchuan Wang a a
Laboratory of Molecular and Materials Simulation, College of Chemical Engineering, Beijing University of Chemical Technology, P.O. Box 100, Beijing 100029, PR China b Department of Physics, Shiga University of Medical Science, Seta Ohtsu, Shiga 520-2192, Japan Received 8 January 2004; accepted 25 February 2004 Available online 26 June 2004
Abstract The shell-model computer simulation is performed to study polarization effects of ions in structural properties of molten ZnCl2. Experimental results of partial structure factors as well as self-diffusion coefficients are well reproduced by a suitable choice of potential parameters. By using the same parameters, time-correlation functions of charge and mass current are evaluated. Discussion is given about broad bands in the infrared absorption spectra, and longitudinal (or transverse) and acoustic (or optic) modes in excitation spectra, with varying core charges of the Cl ion. D 2004 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics simulation; ZnCl2 melt; Polarization effects
1. Introduction Molten ZnCl2 has some uncommon physical properties, such as high viscosity, low conductivity and high tendency to form a glass. For this reason, the properties of ZnCl2 in melt and glass states have attracted special attention, with a number of theoretical and experimental works [1,11]. Biggin and Enderby [2] have used the technique of the 35Cl/37Cl isotopic substitution in neutron diffraction experiments to determine three partial distribution functions, whereby assisted the interpretation for characteristic structures of molten ZnCl2. Wong and Lytle [3] have used the extended X-ray absorption fine structure to study both glass and liquid states. Recently, Allen et al. [4] have revealed by using the technique of pulsed neutron diffractions that the basic ZnCl4 structural units are extremely well defined, and the coordination number decreases from 3.9 to 3.7 with temperature on heating from 330 to 660 jC. The properties of molten ZnCl2 in the temperature range of 623 –853 K are also investigated by using high energy-photon diffractions. It is found that * Corresponding author. Tel.: +86-010-64444905; fax: +86-01064427616. E-mail address: [email protected] (S. Huang). 0167-7322/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2004.02.045
the first peak position of Zn – Cl in partial distribution functions does not change with increasing temperature, but the first shell Zn – Zn distance shows some decrease [5]. All these experimental facts lead us to the microscopic understanding that the small Zn ion occupies tetrahedral coordination sites formed by the close Cl ions, and the tetrahedral coordination remains well-defined above the melting point. Another characteristic feature of molten ZnCl2 is a first sharp diffraction peak (FSDP) at scattering ˚ 1, which is interpreted as intermediate vector about 1 A range order (IRO) in the melt. Molecular dynamics simulation has been widely used as complementary but powerful method to study the structure and dynamic properties of liquids and glasses. For molten ZnCl2, the interaction between two ions is complicated due to the polarization of ions or covalent bond to some degree. The development of a potential model consistent with experimental results is therefore a challenging task. In the first approach to ZnCl2 melt, Woodcock et al. [6] have constructed a potential the Born-Mayer form (called the WAC potential model). Their result for the distance of Zn – Zn is however too large compared to experimental results. The KDR potential model introduced by Kumba et al. [7] including a dispersion term between zinc ions seems to generate the Zn – Zn correlation functions, but there is no
82
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
physical meaning for the zinc ion radius. The authors have used the compressed potential model for simulating molten ZnCl2, but could not predict accurately the Zn –Zn correlation function. Pusztai and McGreevy [8] have applied the reverse Monte Carlo (RMC) simulation method to analyze the structure of the ZnCl2 molten and glass state. Their results show that the structural units (ZnCl4 tetrahedral) are connected at corners. Bassen et al. [9] have used the Monte Carlo and the RMC method to perform a simulation of molten zinc chloride on the basis of the X-diffraction experiments at 623 and 873 K. They have also modified the Coulomb interaction with a dielectric screening function, and found that the simulation using their potential model is capable of reproducing experimental results for correlation functions, the total structure factor and prepeaks. Salmon [10] has systematically analyzed partial structure factors of the molten and glassy MX2 binary system, and found that the first diffraction peaks or pre-peaks have mainly their origin in the number –number partial structure factor SNN(k) and number– concentration partial structure factor SNC(k). Recently, Ribeiro et al. [11] and Gray-Weale et al. [12] have developed a dynamically polarizable model which incorporates a polarization force by combining the Cl dipole polarization with the WAC potential for MX2 melts. They have shown reasonable consistency of the model with experiments for static structures, and further discussed dynamical features in excitation or absorption spectra. In this paper, we use the shell-model potential in the molecular dynamics simulation to study the structural properties of ZnCl2 molten salt. This method is found useful to take into account the polarizability of ions compared to the rigid ion model. The present study has two objectives. The first is to see how our approach using the polarizable potential model is useful by comparing with experimental results of static structures as well as self-diffusion constants. The second, more interesting, is to get an insight into the dynamical structure of this network liquid. It is reported by the instantaneous normal mode analysis that in dynamics the system could not be simply regarded neither as being composed of weakly coupled tetrahedral units with intramolecular type modes, nor of strongly coupled individual ions allowing extended phonons [11]. We discuss polarization effects in dynamical properties by calculating spectra of mass and charge current correlation functions.
is complicated because of its directional character in some degree. Therefore, the polarizable effects or higher body forces must be taken into account. The dynamical shell model recently devised by Fincham and Mitchell [13] incorporates polarizability into the molecular dynamics simulation. In the shell model, the ion is divided into two components, a core with the ion mass and a shell with no mass, which are connected by a harmonic spring. The two parts carry different charges, and the electrostatic interaction between the core and the shell does not include the self interaction (the core and shell of the same ion) [14]. In the simulation of molten zinc chloride, the chloride ion is assumed to have a core and a shell connected by a harmonic spring with the spring constant k. They carry the different charges, but the sum of their charges is equal to 1. The interaction potential /ij(r) between ions i and j is given by ! qi qj Cij r þ Aij exp /ij ðrÞ ¼ 6 qij r r It consists of the long-range Coulomb part and the shortrange repulsion part, with r as the inter-ionic distance, qij the effective hardness and Cij the dispersion parameter. The parameters adopted in the simulation are listed in Table 1: core charges are in units of the elementary charge e, and the value for the Cl ion is taken as 0.984 e from Ref. [15] fitting the crystal properties of ZnCl2. The simulation system consists of 1500 ions (500 Zn and 1000 Cl ions), and the initial position of ions is chosen as CaF2 type in a cubic box. The Ewald summation is used for calculating the electrostatic interaction. The time step is taken as 0.001 ps in the simulation. The temperature is initially controlled to be 2000 K up to 10 000 steps, and it is then reduced to the set temperature. In order to get the density at the set temperature, the simulated system is taken under constant pressure (0.11 MPa) in conjunction with temperature control by using a loose-coupling method, and then the average density is calculated for a constant temperature running. The calculated density at 600 and 873 K is 2.601 and 2.50 g cm 3, respectively. The experimental values at temperatures corresponding to them are 2.52 and 2.40 g cm 3 [6]. It takes 20 ps for the simulation system to achieve an equilibrium state (constant temperature, volume and number of particle), where the energy fluctuation is less than Table 1 Potential parameters for the molecular dynamics simulation [15]
2. Computational detail In the molecular dynamics simulation, each ion is treated as a point particle and the motion of an ion is based on Newton’s equation. The thermodynamic quantities and transport coefficients are calculated by the time average over the ensemble of particles. For the molten and glass state of zinc chloride, the interaction between Zn and Cl ions
Zn Cl
Zn – Zn Zn – Cl Cl – Cl
Q (core) (e)
q (shell) (e)
˚ 2) Spring constant k (eV/A
2.000 0.984
0.000 1.984
17.250
Aij (eV)
˚) U (A
˚ 6) Cij (eV/ A
0.0000 9704.8900 3296.5700
0.0000 0.2320 0.3289
0.0000 0.0000 107.2000
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
83
0.2%. The structure and dynamical quantities are calculated by taking the average over more than 25 ps. All the simulations are performed by using DL_POLY code [16].
3. Results and discussion The static structure of molten systems is usually described by the partial radial distribution function (RDF). It represents the probability of finding another particle from a particle at origin as a function of the radial distance. Fig. 1 shows the partial RDF at T = 1000, 873 and 600 K. Table 2 presents the summary of results for the partial RDF, in comparison with the experiments. The results of the molecular dynamics simulation are found to be in good agreement with experiments and the RMC simulation. The ratio rClCl/ rZnCl for first peak positions from the partial RDF is about 1.61, close to the value of 1.63 for the regular tetrahedral [2,9,17]. The partial structure factors sij(k) for i and j species (i, j = Zn, Cl) are obtained by the Fourier transformation of the partial RDF gij(r) [18] Z l 4p sij ðkÞ ¼ dij þ pffiffiffiffiffiffiffiffi ½gij ðrÞ 1sinðkrÞrdr k ni nj 0 with ni being the number density of i species. The total structure factor is then calculated from sij(k) according to the formula [4] FðkÞ ¼ c2i b2i ½sii ðkÞ 1 þ c2j b2j ½sjj ðkÞ 1 þ 2ci cj bi bj ½sij ðkÞ 1 with bi being the scattering length and ci the concentration of particle i. It is experimentally determined by neutron or X-ray diffractions. The calculated results for sij(k) and F(k) are shown in Figs. 2 and 3, respectively. It is noted that there exists the lower limit of k = 2p/L owing to the finite size L of the simulated system in calculating partial structure factors ˚ 1). We can observe clearly the pre-peak at (kmin = 0.20 A 1 ˚ ), with good agreement between our simulated k f 1(A results and neutron diffraction experiments in the total structure factor. These results imply a validity of our potential model incorporating polarizability to investigate the structure of ZnCl2 molten salt. In the present potential model, the electric field for unlike charges in an ion (Cl ion) due to surrounding ions tends to make the ion to separate the shell from the core, and it gives rise to a polarization dipole. The polarizability of the particle in the equilibrium depends on the core charge. We have therefore changed the core charge from the set value of 0.984 e in Table 1 down to 0 to see polarization effects. Fig. 4 shows the calculated result for the dependence of the partial RDF gZn – Zn on the core charge. It is found that both the position and the height of the first peak remarkably decrease as the core charge (and thus polarizability)
Fig. 1. Partial radial distribution function : (a) 1000 K; (b) 873 K; (c) 600 K; Line: gZn – Cl(r) Dot: gCl – Cl(r) ; Dash: gZn – Zn(r).
increases. From micro-structural viewpoint, these facts mainly reflect the appearance of the structure of the [ZnCl4] tetrahedral, whose units are connected by the chlorine bridge or edge [9]. The Coulomb repulsion between doubly charged cations appears, in effect, to be significantly reduced by the polarization effect. The connection of the
84
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
Table 2 Summary of static results Temperature (K)
Method
600 873
MD MD Neutron2 RMC9 RMC9 X-ray17
623 873
First peak position (nm)
Coordination
Zn – Zn
Zn – Cl
Cl – Cl
Zn – Zn
Zn – Cl
Cl – Cl
0.39 0.39 0.38 0.385 0.385 0.366
0.23 0.23 0.229 0.220 0.225 0.229
0.37 0.37 0.371 0.355 0.360 0.385
4.271 4.359 4.7
40.080 4.106 4.3
10.377 10.124 8.6
[ZnCl4] tetrahedral, or the tendency for forming a network structure, would get stronger with increasing polarizability. The diffusion process of ions can be investigated through its mean-square displacement as a function of time. The self-diffusion constant is evaluated by the formula < Arðt þ sÞ rðsÞA2 > t!l 6t where t is the time interval, and the angular bracket stands for the ensemble average. It is evaluated by using fitting procedures during 6 ps after 20 ps required to be achieved the equilibrium in the molecular dynamics simulation. The simulated result obtained by using the values of the parameters in Table 1 (the core charge is 0.984 e) is shown as a function of temperature in Fig. 5. It is observed that the value of the self-diffusion coefficient for Zn is always smaller than that for Cl in this temperature range, and the behavior appears similar for both Zn and Cl in the dependence on temperature. This is because the formation of the [ZnCl4] unit or units connected to form more complex structure restricts the mobility of ions. The experimental results for the self-diffusion coefficient are 0.010 and 0.013 (unit is 1 10 5 cm2/s) for Zn and Cl ions at 600 K [19], respectively. The corresponding values of our simulation are 0.04241 and 0.08069 (unit is 1 10 5 cm2/s) for Zn and Cl ions, respectively, which corresponds with the same magnitude of experimental value. D ¼ lim
Fig. 2. Partial structure factors at 600 K. Solid line: sZn – Cl(k), Dash-Dot: sZn – Zn(k), Dot: sCl – Cl(k).
4
4
12
The diffusive motion of ions is not independent of polarization effects, which is responsible for the formation of the tetrahedral unit. Fig. 6 shows the dependence of diffusion coefficients on the core charge at T = 873 K. It is found that diffusion coefficients of both ions show similar behavior of increase with the core charge. It depends on temperature, but the result may be interpreted in the way that the increase of the core charge facilitates the formation of many tetrahedral units, but these units would be rather loosely connected with each other. As a result, ions in the molecular unit appear to diffuse more easily, compared to the case of the individual motion at zero core charge. The electrical conductivity is related to the correlation function of the total charge current correlation function J(t) via its Laplace transform a(x), defined by the following equations [19] X jðtÞ ¼ qa via i;a
< jðt þ sÞjðsÞ > < jðsÞjðsÞ > Z l aðxÞ ¼ expðixtÞJ ðtÞ
J ðtÞ ¼
0
where j(t) is the total charge current, qa is the charge of a species and via is the velocity of particle i of a species.
Fig. 3. Total structure factor at 600 K, with the experimental data from Ref. [4].
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
Fig. 4. Partial radial distribution function of Zn – Zn as a function of core charge of Cl ion at 873 K.
85
Fig. 6. Self-diffusion coefficient as a function of core charge at 873 K.
The results of the molecular dynamics simulation for a(x) at representative temperatures ranging from 600 to 1200 K are shown in Fig. 7. The total charge current correlation function reflects a collective property of the system, and according to the fluctuation theory, the infrared absorption can be obtained from the Fourier transformation of the dipole correlation function, which is also related with a(x). There are two broad bands in Fig. 7; one band is around 3.141 THz ( f 104.77 cm 1), and the other is around 10.050 THz ( f 335.23 cm 1). As the temperature increases, the high frequency band makes a shift towards the low frequency region with broader width. Angell and Wong [20] have measured the infrared spectra of ZnCl2 molten salt, and observed two broad bands at 100– 115 and at 255 cm 1. The position of the high frequency band in our results is found to be rather higher in frequency compared with the experimental result. A similar situation was reported for the dynamically polarizable model in comparison with experiments [12]. The peak positions in the Fourier transform of the total current correlation function represent the transverse optic
frequency. The mass and charge currents are defined by the equations [21]
Fig. 5. Self-diffusion constant as a function of temperature.
Fig. 7. Spectrum of the total charge current correlation function.
Jm ðk; tÞ ¼
N X
expðik ri ðtÞÞmi vi ðtÞ
i¼1
Jq ðk; tÞ ¼
N X
expðik ri ðtÞÞqi vi ðtÞ
i¼1
Here, mi and qi are the mass and charge of ion i, respectively. The transverse and longitudinal correlation functions for the mass and charge currents are given by T CXX ðk; tÞ ¼ ð1=2N Þ < k JX ðk; t þ sÞ k JX ðk; sÞ > L ðk; tÞ ¼ ð1=N Þ < k JX ðk; t þ sÞk JX ðk; sÞ > CXX
with X = m for the mass current and X = q for the charge current, where the symbol denotes the vector product. When the spectra or Fourier transforms of the transverse and longitudinal mass correlation function exhibit sharp peaks, the peak position is identified as the transverse and longi-
86
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
˚ 1), Dot: k = 1.0 (A ˚ 1). (a) TA, Fig. 8. Spectrum of transverse and longitudinal mass and charge correlation functions at T = 873 K. Solid line: k = 0.45 (A (b) LA, (c) TO, (d) LO.
tudinal acoustic (TA and LA) frequency, respectively. Similarly, peaks in the spectrum for the charge current correlation function imply the transverse and longitudinal optic mode’s frequencies (TO and LO). We have calculated these spectra, and the results at ˚ 1) and k = 1.0 (A ˚ 1), T = 873 K are shown at k = 0.45 (A in Fig. 8. It is noted that the latter is the wave number at prepeaks indicating the specific structure of the system. In the spectrum of acoustic modes, there are many fine structures for the LA mode, compared to the TA mode. It may correspond to the result by the INM that normal modes are not interpreted by single plane waves [12], or the density fluctuation with a definite wave number is not a welldefined variable. Both the TA and LA modes make a shift to higher frequency as the wave number k increases, which implies that these modes have a positive dispersion. As for optic modes representing the relative motion of ions, the peaks of TO and LO modes correspond to the two peaks in Fig. 7. It is found that these modes tend to slightly show a negative dispersion. These are in accordance with the results for similar systems [12]. We have further investigated the polarization effects on the excitation spectra in Fig. 9. The calculated results at T = 873 K for the transverse (or longitudinal) mass (or ˚ 1), charge) correlation function are shown at k = 0.28 (A
with representative values of the Cl ion core charge, i.e., 0, 0.5 e and 0.984 e. The complicated dependence on the core charge is seen for these spectra. For the acoustic spectra, the position of the main peak for 0.984 e seems to make a shift to lower frequency, compared to the result for 0.5 e. It implies an effective softening of the coupling constant for individual ions, presumably averaged over tight coupling in the tetrahedral unit, and relatively weak coupling between the units. For the optic spectra, it is found by comparing the cases of 0.984 e and 0.5 e that the two large peaks tend to separate, indicating a kind of the level repulsion, the results for the zero core charge show different behavior for the acoustic and optic spectra. It may be presumably correlated to the rapid change of the diffusion coefficient accompanied by the increase of the core charge from 0 to 0.5 e, as seen in Fig. 6.
4. Conclusion We have discussed the polarization effects in the static and dynamical properties of molten ZnCl2 by performing the shell-model computer simulation. It is found that experimental results of the total structure factor and diffusion constants are fairly well reproduced by appropri-
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
87
˚ 1). The Fig. 9. Dependence on the core charge for spectra of transverse and longitudinal mass and charge correlation functions with T = 783 K, at k = 0.28 (A Cl ion core charge is: Dash dot: 0.00 e; Dot: 0.500 e; Solid line: 0.984 e.
ately choosing parameters in the effective interionic potential and core charge for the ion (only for Cl ion). In addition, by investigating the dependence of the partial structure factor of Zn –Zn on the core charge, we have explicitly shown a significant contribution of the polarizability of ions: the ratio of the first peak positions for partial structure factors is found to be consistent with the formation of the tetrahedral unit [ZnCl4]. We have moreover found that the self-diffusion coefficients of Zn and Cl show the similar dependence on the core charge, which also might indicate the molecular structure mentioned above in the network liquid. As for dynamical properties, we have calculated the spectrum of the total mass and charge current correlation function. It is found by adopting a few values for the core charge that excitation spectrum for acoustic and optic modes are qualitatively understandable by the polarization effects. On the whole, our results are in accordance with other studies such as the instantaneous mode analysis. In comparison with observed infrared absorption spectra, however, some discrepancies are found for the two broad bands. This and similar results reported in the dynamically polarizable model suggest partly an inadequacy of polarization models in treating the bonding.
The significant dependence of the diffusion coefficient on the core charge of the Cl ion shown in Fig. 6 indicates that polarization effects coupled with temperature would improve the situation. For this and a more detailed understanding of the system, it is interesting to make an ab initio calculation for the structure, which will be done in the future.
Acknowledgements This work is supported by the National Basic Research Program of China, Grant No. 2003CB615807, and the Nation Natural Sciences Foundation of China (No. 20236010).
References [1] [2] [3] [4] [5]
M. Rover, M.P. Tosi, Rep. Prog. Phys. 49 (1986) 10001. S. Biggin, J.E. Enderby, J. Phys. C 14 (1981) 3129. J. Wong, F.W. Lytle, J. Non-Cryst. Solids 37 (1980) 273. D.A. Allen, R.A. Howe, N.D. Wood, J. Chem. Phys. 94 (1991) 5071. J. Neuefeind, K. Todheide, A. Lemke, H. Bertagnolli, J. Non-Cryst. Solids 224 (1998) 205.
88
S. Huang et al. / Journal of Molecular Liquids 115 (2004) 81–88
[6] L.V. Woodcock, C.A. Angell, P. Cheeseman, J. Chem. Phys. 4 (1976) 1565. [7] P.N. Kumta, P.A. Deymier, S.H. Risbud, Physica. B + C 153 (1986) 85. [8] L. Pusztai, R.L. McGreevy, J. Non-Cryst. Solids 117/118 (1990) 627. [9] A. Bassen, A. Lemke, H. Bertagnolli, Phys. Chem., Chem. Phys. 2 (2000) 1445. [10] P. Salmon, Proc. R. Soc. Lond., A 437 (1992) 591. [11] M.C.C. Ribeiro, M. Wilson, P.A. Madden, J. Chem. Phys. 109 (1998) 9859. [12] A. Gray-Weale, P. Madden, M. Wilson, J. Chem. Phys. 113 (2000) 6782.
[13] [14] [15] [16] [17] [18]
D. Fincham, P.J. Mitchell, J. Phys., Condens. Matter 5 (1993) 1031. P.J.D. Linda, M.J. Gillan, J. Phys., Condens. Matter 5 (1993) 1019. D.J. Binks, PhD thesis, University of Surrey, 1994. W. Smith, T.R. Forester, The DL_POLY user manual (1999) 55. R. Triolo, A.H. Narten, J. Chem. Phys. 14 (1981) 703. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publication, New York, 1996. [19] J.P. Hansen, R. McDonald, Theory of Simple Liquids, 2nd ed., Academic Press, London, 1986. [20] C.A. Angell, J. Wong, J. Chem. Phys. 53 (1970) 2053. [21] B.J. Berne, R. Pecora, Dynamic Light Scattering, Wiley, New York, 1976, Sec 15.3.