Molecular dynamics computer simulations of molten zinc chloride

Physica 131B (1985) 227-233 North-Holland, Amsterdam

M O L E C U L A R DYNAMICS C O M P U T E R S I M U L A T I O N S OF MOLTEN ZINC CHLORIDE P.J. G A R D N E R and D.M. H E Y E S Department of Chemistry, Holloway College, University of London, Egham, Surrey, TW20 OEX, UK

Molecular dynamics computer simulations are presented of model ZnC12 at a temperature of 1200K and molar volume of 61 cm3mo1-1. Pair radial distribution functions, coordination numbers, energy, internal pressure and self-diffusion coefficients are used as probes for the effects of model parameters on calculated properties. Pairwise additive rigid ion potentials are successful in reproducing many features of the real melt. A split like-ion first peak in the pair radial distribution function, however, is in marked contrast to recent neutron scattering work which gives approximate coincidence of Zn-Zn and CI-CI first peaks at (3.57-+0.1)A. o

1. Introduction

T h e r e is an increasing technological use for molten ZnCI 2 as a solvent and catalyst in the hydrocracking of coal slurries and heavy crude oil fractions [1]. It is also interesting in other respects. For example close to its melting point (~590 K) it exhibits an unusually high viscosity, low ionic diffusivities and electrical conductivity and ready glass formation on cooling, which are in m a r k e d contrast to the behaviour of many other metal halides [2-5]. Anomalously low enthalpy and volume changes on fusion are consistent with the persistence of a crystalline-like structure in the melt. This has been [6] expressed in terms of extended 'crosslinked' network structures, (ZnCIz) n CI- which gradually collapse to smaller units e.g. Zn2CI62-, ZnC142-, ZnCI 3etc. as the t e m p e r a t u r e increases. Currently, Molecular Dynamics, MD, simulations have produced no definitive evidence for long-lived ionic aggregates such as these [7]. Indeed the role of covalency versus ionicity even in such compounds as SiO 2, U O 2 and CuCI is still hotly debated just as much for the solid state as for the liquid state [8]. Although angle-dependent bending forces are required for certain properties such as the dispersion curve of SiO 2, it is well established that the pair radial distribution function, g(r), of SiO 2 melts, for example, can be

adequately reproduced in terms of a purely ionic model (containing Si 4+ and 0 2- here) [4]. The metal dihalides are of particular interest here. Neutron scattering experiments using isotopic substitution have recently been used to obtain the c o m p o n e n t pair radial distribution functions, g+_, g++ and g__ (via the structure factors) of a series of metal dichlorides. Two features in these pair radial distribution functions provide evidence for deviations from a simple ionic model. First, because the cation has twice the absolute charge of the anion it is to be expected that the first peak in g++ should be significantly further from the origin than that of g__. This is observed both experimentally and by c o m p u t e r simulation for molten SrCI 2 and BaCI z [9]. H o w e v e r for ZnCi 2, MgCI 2 and CaC12 there is near coincidence of g+÷ and g__ within the first coordination shell [9-11]. T h e results from X - R a y diffraction for ZnCI 2 [12] support this near coincidence. Secondly, relationships between the first peaks of g+_ and g__ indicate other anomalous behaviour. A wide range of sizes of cationic (M 2÷) radii have been considered experimentally by neutron scattering: Mg (0.66), Zn (0_.74), Mn (0.80), Ca (0.99), Sr (1.12) and Ba (1.34 ~ ) [9-11]. (These are Pauling ionic radii which are inversely proportional to an effective nuclear charge q* =

0378-4363/85/$03.30 O Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

228

P.J. G a r d n e r a n d D . M . H e y e s / M o l t e n z i n c chloride

q'-s

where q' is the atomic number and s is a screening coefficient for electrons, which are usually evaluated by Slater's empirical rules). These cations are all somewhat smaller than C1(1.81 ~). Integration of the first g+_(r) peak gives 4-fold coordination for all but only for ZnCI 2 is there evidence of a true tetrahedral arrangement of "just-touching" anions about a cation. If r+_ and r_m_ are the positions of the first maxima of g+_ and g__ then tetrahedral packing gives rm 1.63, the value obtained - - / l r +m_ = ~ v / 8 ~ : experimentally for ZnCI 2. For all other metals this ratio lies between 1.3 and 1.5. For example [10] 1.47_+0.03 for MgCI 2 and 1.43+0.03 for MnC12. For CaCI2, SrCi 2 and BaC12 the values of 1.35, 1.35 and 1.3 are derived [10, 11]. These structures are essentially 4-fold coordinated but not close packed. There is also an inconsistent relationship between this ratio for MgCI 2 and ZnCI 2. Magnesium dichloride has a smaller cation than ZnC12 and a smaller r m__./rm+_ value, which is a result expected for a larger cation however. The anomalous behaviour of MnCI 2, MgC12 and ZnC12 is also manifest in the solid phase in which they adopt low coordination layered structures indicative of some form of "covalent" character. SrC12 in contrast crystallises in a fluorite structure and BaCI 2 has the PbCI z structure. The resolution of such inconsistencies and the provision of thermodynamic and transport data for ZnC1 z (for which there is a serious shortage) provided the motivation for this work.

2. T h e m e t h o d

It is usually assumed in MD that the interaction between the particles in the system is pairwise additive. We use the following pair potential form [4]:

dpij(r ) = qiqj e2+ bi/ exp((o-ij - r)/p).

(1)

r

Here r is the distance between ions i and j, and qi is the ionic charge ( - 1 for CI- and +2 for

ZnZ+). The repulsive term is composed of cr~j, the sum of ionic radii and b~j which is a hardness parameter. These parameters for a wide spread of possible potentials are given in table I. In MD the Newtonian equations of motion of N interacting particles are usually solved numerically using Verlet's Leapfrog algorithm [13, 14]. The instantaneous temperature at time t is given by N

T(t) = (3NkB)-I ~'~ A4N~E(t),

(2)

i=1

where M÷ = 65.37 g mol -~ and M_ = 35.453 g mol 1, and where V/ is the velocity of ion i. The internal pressure of the system is given by

(3) where ~b : - r

0~b;i

(4)

Or

is called the pair virial function. The internal energy is given by Table I V a l u e s for the p a r a m e t e r s of the ZnCI2 p o t e n t i a l s c o n s i d e r e d in this w o r k . T h e W A C p o t e n t i a l w a s o b t a i n e d f r o m eq. (3) in ref. 4 t a k i n g z+ = 0, z = - 1 a n d n_ = 8. T h e p o t e n t i a l G H 1 c o r r e s p o n d s to the f o r m of eq. (3) in ref. 4 with z+ = 2, n+ = 8, z_ = 1 a n d n÷ = 8. e is the a b s o l u t e e l e c t r o n i c charge. N T is the n u m b e r of t i m e steps u s e d in the p r o d u c t i o n simulation Parameters

WAC

GH1

GH2

GH3

Units

q÷ q b+÷ b-b+ o'+÷ o'_ _ o-+_ p

+2 -1 0.190 0.143 0.166 2.86 3.356 3.108 0.34

+2 -1 0.285 0.2375 0.26125 2.86 3.356 3.108 0.34

+1.5 -0.75 0.190 0.143 0.166 2.86 3.356 3.108 0.34

+1.25 -0.625 0.190 0.143 0.166 2.86 3.356 3.108 0.34

e e 10-19 J 10-19 J 10 -19 J .~ ,~ ,~ ,~

NT

8000

2000

4000

8000

P.J. Gardner and D.M. Heyes I Molten zinc chloride

(5)

g(r)

= 1 ~,, ~, P=PBg=~(r) p2

The constant volume specific heat per mole, Cv, is

CJ3NAka = 1.5(1-1.5N(STZ)(T)-2) -1 ,

(6)

where N A is Avogadro's Number and where

8 T = T(t) - (T(I)),

229

=

a

(12)

8 4

4

~g++(r) + ~g+_(r) + ~ g _ _ ( r )

(13)

in this case. Then n=a(r) are evaluated in histogram form at each time step of the calculation. The melt structures are also characterised by radial coordination numbers,

(7) r

and the relaxational bulk modulus is

V N - f r2g+-(r) dr'

n+_(r)

(14)

0

K ® - K~o = k - ~

(SP2(t))'

(8)

r

4zr (r)=~(N_-l)f

r2g__(r) dr,

(15)

47r N n + + ( r ) = - ~ ( + - l ) f r2g÷+(r) dr,

(16)

n where 8P(t) = P ( t ) - (P(t)),

0

(9)

The equilibrium melt structure is analysed in terms of pair radial distribution functions. These functions measure the factor by which the local density p(r) at some distance r from an ion in the melt deviates from the bulk density, p. If the average number of ions in a spherical shell of thickness dr at a distance r from a specific partiCle is denoted by n(r) then the total pair radial distribution function g(r) is defined by

4zrr2pg(r) dr = n(r) .

(lO)

In the MX 2 melt three components of g(r) are necessary for a complete description of the pair correlations. There are the distribution functions for pairs of cations: g÷÷(r), pairs of anions g__(r) and pairs of unlike ions, g+_(r). The distribution function g~(r) is defined by the relation

4~rrEpag~a(r) dr = %a(r) ,

(11)

where pa is the density of species fl and n~a(r) is the average number of ions fl in a spherical shell of radius r and thickness dr around an ion of species t~. Therefore the overall distribution function g(r) is

r

0

where N+ and N_ are the numbers of cations and anions, respectively, V is the volume of the MD cell. During the MD the electrostatic contributions to the potential ~b=Land forces cannot be evaluated by direct summation as the terms in the summation converge too slowly. The Ewald method was used to circumvent this problem (15). The Ewald potential, &~l, is qb=l" = qiqj r

e 2 erfc(ar)

2q,a 77"112 "

No reciprocal space terms were included here; a = 5.6/L. The MD time step was 0.65 × 10 -14 S, the molar volume, Vm, was 61 cm3mo1-1 (the experimental value at 1200 K [4]) and the chosen temperature, T, was 1200 K. The number of ions in the MD cell was 324 (216 anions and 108 cations) which was created by expanding an arbitrary array of 12 ions by 3 x 3 x 3. The MD cell was cubic and had a sidelength, L, of 22.2/~. The truncation distance for the inter-

230

P.J. G a r d n e r a n d D.M. Heyes / Molten zinc chloride

actions was L/2 (=11.1~,). This is sufficient so that the potential ~bij always is only a few percent of kBT (viz ~ 80 K) at the cut-off radius.

3. Results and discussion In ionic melts which contain multi-valent cations it is well established that competition for tetrahedral coordination with the more abundant anions is a major feature in determining structural features [16, 17]. In line with this trend

ZnC! 2 also favours a 4-fold coordination. In any M D treatment it is important to establish the effect of changes in the pair potential in influencing this process. Fig. 1 contains the ~b+ , ~b++ and ~b_ for the four potentials considered here and characterised explicitly in table 1. Potential GH1 is more repulsive than the W A C potential [4] due to the enhanced b~. In contrast unlike-ion potentials GH2 and GH3 are shallower than the W A C example because of a reduction in the ionic charges. The electroneutrality of the melt is maintained in each case.

(a)

12-

(b) @H1

I/ ~/,~1++ 8-

uo-%-: .L

+.

,

I

'1

,

I

,

2

/

-4

I

,

\

\

I

/Z + ~

8

i

[

i

I

i

8

/

-4

r(A)

r(~,) [c)

12

(d)

I

12

t-+

-l/i-+ 8

+_

OH2

@H3

8

%4.

z+ @

, I

-4. -4

2

4

6

8

rim

Fig. 1. The inter ionic pair potentials used in this work: + denotes cation, - denotes anion; (a) W A C [4]; (b) G H I ; (c) GH2; (d) GH3. The r -l Coulomb term is replaced by e r f c ( a r ) r -1 as described in the text.

P.J. Gardner and D.M. Heyes / Molten zinc chloride

Fig. 2 presents the component pair radial distribution functions for model molten ZnCl 2 using the Woodcock, Angell and Cheeseman (WAC) potential. Also shown on this figure are the experimental g~a(r) [9, 11]. There are a number of significant discrepancies. The first peak in g+_(r) is at 2.3A from experiment when compared with 2.4 ~ for the simulation. This could possibly be due to the difference in temperature between the experimental and simulation temperatures; 600K and 1200K respectively (the latter was chosen to aid rapid equilibration and to correspond to a state point treated in the original simulation with the WAC potential [4]). The heights of these peaks are nevertheless in good agreement. The major difference which is, however, unlikely to be due to a temperature change, is the position of the first maximum in g++ for the simulation which occurs at 4.6 ~ in contrast to 3,8/~, by experiment [11]. The g__ for experiment and simulation are at 3.71 and 3.65/~ respectively which are hardly distinguishable. Other preliminary MD simulations of ZnCI 2 have also failed to generate a realistic Zn-Zn component of the pair radial distribution function and is indicative of a major flaw in the simple ionic model [18-20].

231

Potential GH1 has more repulsive short range terms than WAC (see fig~ 3). Consequently the first peak in g+_ is at 2.7 A as opposed to 2.4 for the WAC simulation, The first peak in g__ is little changed however at 3.7 ~k. The first peak in g++ again shows more sensitivity to the potential form, appearing at 4.85 ~ for GH1 but 4.6/~ for WAC. The q+ and q_ of potential WAC of table I were changed to 1.5 and -0.75 respectively in GH2 (see fig. 4). These maintain the electroneutrality of the melt but weaken the like-ion repulsions. Unfortunately there is still quite a large separation between g÷÷ and g__ first peaks (at 4.9 and 3.7 ~ respectively). The major effect is on g÷_ which has a broader and smaller first peak. This is accompanied by a reduction in total

n,_(r)

GHI

gij(r)

__ ++ f'\ /

\

/

(a)

O| M.D.

n+_[r)

I

2

N.S. expt. oe°oo

J

z,

6 8 r(A] r .......

I0_~

q- --

Fig. 3. A s for fig. 2 except that the GH1 potential was used.

{b)

0 5

,

~

,

I/

2. •

I

3

,

I

4

-.-

-

+-r

[a)

4 gij[ r)

-'I

3

:i

2

],l

0

'

:

--

+,

WAG

.... 1 2[~3

, 4

/;, 2

a

6

8

10

OH2

++

gij[r) 21

r{A) Fig. 2. (a) unlike ion radial coordination numbers, n+_(r) for the M D simulation with the W A C ionic potentials; (b) the corresponding c o m p o n e n t pair radial distribution functions + - (), - - ( - - - - ) and + + ( ) are compared with the experimental functions [11].

0

i

i

2

i

4

6

8

10

r(A) Fig. 4. A s for fig. 2 except that the G H 2 potential was used.

l

P.J. Gardner and D.M. Heyes / Molten zinc chloride

232

energy from -2132 to -1060 kJmol i (see table II). This trend is continued further by adopting q+ = 1.25 and q_ = - 0.625 for potential G H 3 (see fig. 5). Interestingly the relative positions of the g++ and g__ first peaks are much the same here and change little down this sequence. This is unexpected as the difference between like-ion repulsions is dropping markedly down this sequence. In fact, the positions of all the first peaks in g+_, g++ and g are all very close for GH1, G H 2 and GH3. The potentials GH1, G H 2 and G H 3 generate less cohesively bound liquids so that the positions and spreads of the g,,o peaks are generally increased (the g__ however is least affected by these changes: r_m_ = 3.7 in each case whereas rm ~ 2 . 7 and r m The n+ (r) also reflect ++ + ~4.9). this relaxation of strongly bonded units. These changes also influence the thermodynamic and transport properties of these systems, summarised in table II. The total energies decrease along the series W A C - G H 3 , all are at least 500kJmol -] smaller than experiment however [4]. The pressure increases also. The selfdiffusion coefficients (obtained from the mean square displacements of the ions [4]) minimise at GH1 and then increase slightly through G H 2 to GH3. These values are all within experimental uncertainties however.

4. C o n c l u s i o n s

It is tempting to speculate that the origin of these differences (indeed any unexpected lea-

(a) n+_(r)4

(b)

/-*~I 0/

gutr)

I

_._L

I

I

I

Z////II

I

OH3

I

+--

.....

t /

;I

2

/

I

I

J

I

4

6

Potential WAC GH1 GH2 GH3 Experiment

-2132 -1993 -1060 -672 -2604 ")

a) Ref. [4].

I

I

t

10

r(A) Fig. 5. As for fig. 2 except that the OH3 potential was used.

tures!) lies in chemical bonding of a nonpairwise additive nature. This is a necessary condition because the strength of a covalent and purely ionic " b o n d " can be comparable and it would therefore be difficult to distinguish them if viewed on a pairwise additive basis only. An estimate of the fractional ionic character, f, as defined by Pauling is [21] faB = 1 - exp(--0.250¢A -- XB)2) where gi is the electronegativity of species i ( = A or B). Now, )¢zn2+= 1.65 and Xo = 3.16 (22) so f = 0.43, i.e. the bond is nominally only 43% ionic. An ion in a tetrahedral site develops an induced octupole moment which can be viewed as equivalent to a covalent (viz. angle dependent) bond [23], but this would necessitate a large increase in computing time and would involve additional uncertainties with parameters. (e.g. an

Table II Derived quantities from the simulations. 1 GPa = 10 kbar Total energy (kJ mol i)

i

8

P (kbar)

D+ D (10 5cm2s-I)

3.6 9.4 14 24

2.4 1.7 4.2 4.8 1-9

3.0 2.9 8.1 9.9 1-15 a)

C

J~

(J moV 1 K i)

(GPa)

93 81 69

25 18 9

K0 s

P.J. Gardner and D.M. Heyes / Molten zinc chloride

octupole polarizability for Zn2+). An alternative route would be to parameterise near HartreeFock-limit quantum mechanical calculations of the ZnnCl m potential surface, But this is a computationally demanding task [24]. At present empirical modifications to the pair wise additive potential have been made in order to gain insights into possible adjustments which would give coincidence of r~+ and r_m_.It appears that ~b÷+ must decay more rapidly with distance than the examples considered. This would probably also have the associated effect of improving the agreement between the experimental and simulation total energy, U. The calculations were confined to the pairwise additive potential framework and demonstrate that the melt structure is a delicate balance between strong manybody interactions. Small changes in the form of the potential produce unexpected results because the ~b~a are so large when compared with kBT. This makes the task before us imposing indeed.

Acknowledgements D.M.H. gratefully acknowledges The Royal Society for a Research Fellowship. References [1] R.T. Struck and C.W. Zielke, Fuel 60 (1981) 795. C.W. Zielke and W.A. Rosenhover, U.S. Patent No. 4424111, Jan 3, 1984. [2] C.A. Angell and J. Wong, J. Chem. Phys. 53 (1970) 2053. G.J. Janz, F.W. Dampier, G.R. Lakshminarayarian, P.K. Lorenz and R.P.T. Tomkins, Molten Salts, NSRDS-NBS 15 (Washington D.C., 1968). [3] C.A. Angell, E. Williams, K.J. Rao and J.C. Tucker, J. Phys. Chem. 81 (1977) 238. [4] L.V. Woodcock, C.A. Angell and P. Cheesemen, J. Chem. Phys. 65 (1976) 1565.

233

S.K. Mitra, Phil. Mag. B 47 (1983) L63. M.-L. Saboungi, A. Rahman and M. Blander, J. Chem. Phys. 80 (1984) 2141. [5} C.A. Angell, I.M, Hodge and P.A. Cheeseman, Proc. Int. Symp. Molten Salts (Electrochem. Soc., Princeton, N.J. Electrochemical Soc., 1976) p. 138. [6] H. Bloom, The Chemistry of Molten Salts (Benjamin, New York, 1967). W.R. Busing, Trans. Amer. Cryst. Assoc. 6 (1970) 57. [7] S.W. de Leeuw, Mol. Phys. 36 (1978) 103, 37 (1979) 489. B.R. Sundheim, J. Chem. Phys. 73 (1980) 2474. [8] N.H. March and M.P. Tosi, Coulomb Liquids (Academic Press, New York, 1983). See also M. Leslie and A.C. Lasaga, presented at the NATO Adv. Res. Workshop on Computer Simulation of Condensed Matter (Norwich, 1984). [9] J.E. Enderby, Contemp. Phys. 24 (1983) 561. R.L. McGreevy and E.W.J. Mitchell, J. Phys. C. 15 (1982) 4635. [10] S. Biggin, M. Gay and J.E. Enderby, J. Phys. C: Solid State Phys. 17 (1984) 977. S. Biggin and J.E. Enderby, J. Phys. C.: Solid State Phys. 14 (1981) 3577. [11] S. Biggin and J.E. Enderby, J. Phys. C.: Solid State Phys., 14 (1981) 3129. [12] R. Triolo and A.H. Narten, J. Chem. Phys. 74 (1981) 703. [13] L. Verlet, Phys. Rev. 159 (1967) 98. [14] D.M. Heyes, Chem. Phys. 82 (1983) 285. [15] D.M. Heyes, J. Chem. Phys. 74 (1981) i924. [16] T.F. Soules, J. Chem. Phys. 71 (1979) 4570. [17] S.H. Garofalini, J. Chem. Phys. 78 (1983) 2069. [18] J.A.E. Desa, A.C. Wright, J. Wong and R.N. Sinclair J. Non-Cryst. Solids 51 (1982) 57. [19] A. Thompson, Ph.D. Thesis, University of Sheffield, U.K., 1982. [20] P.A. Cheeseman, Ph.D. Thesis, Purdue University, U.S.A., 1980. [21] C.R.A. Catlow and A.M. Stoneham J. Phys. C: Solid State Phys. 16 (1983) 4321. [22] J.E. Huheey, Inorganic Chemistry (Harper and Row). L. Pauling, Nature of the Chemical Bond, Corneli University press., 3rd ed., 1960. [23] G.D. Mahan, Chem. Phys. Lett. 76 (1980) 183. [24] D.G. Bounds (personal communication). See also D.G. Bounds and P.J. Bounds, Mol. Phys. 50 (1983) 25. D.G. Bounds, Chem. Phys. Lett. 86 (1982) 1.