Thermodynamic Properties and Self-diffusion of Molten Sodium Chloride A Molecular Dynamics Study BY JOHNW. E. LEWIS*-/-SINGERAND KONRAD Department of Chemistry, Royal Holloway College, Egham Hill, Egham, Surrey Received 15th May, 1974 Molecular dynamics calculations are reported for liquid sodium chloride simulated by Born- Huggins-Mayer type pair potentials. Thermodynamic properties and self diffusion constants have been obtained for 22 VT points ranging between 1000 and 2000 K and between -1 and 14 kbars. The agreement with experimental data is, on the whole, satisfactory. Algebraic expressions have been fitted to the calculated data for the pressure, internal energy, the (variation of the) Helmholtz free energy, and the ionic self diffusion constants.The variation of the pair correlation functions and of the velocity autocorrelation functions with V and T are examined. Some features of computer generated films showing the ionic motion are briefly discussed. Computer simulation of liquid systems by the Monte Carlo (MC) and molecular dynamics (MD) methods have led to a considerably better insight into the structure of simple 1iquids.l The application of these methods to hard spheres and Lennard- Jones atoms has been well exploited and the results of these calculations have had a considerable impact on the development of theories of liquids.2* Woodcock and Singer have applied the MC method to molten potassium chloride and more recently MC calculations for a range of liquid and solid alkali metal halides have been carried Limited calculations by a modified MD method have been reported by Woodcock.6a Rahman and co-workers have recently published some MD calculations for beryllium fluoride and lithium fluoride, which show among other things that a pair potential with purely electrostatic cohesive terms correctly predicts the tetrahedral coordination and high viscosity of liquid beryllium fluoride.' Lantelme et al.have reported limited MD calculations for liquid sodium chloride.8 Following the success of the MC calculations on molten potassium ~hloride,~ it was decided to attempt a molecular dynamics study in order to see whether this tech- nique could give estimates for the transport properties of a molten salt system. It seemed reasonable to first try and calculate the self diffusion constants and to compare these to experiment.Molten sodium chloride was chosen for study for two reasons, first the availability of experimental data and second because the self diffusion constants are different for sodium and chloride ion in the melt in contrast to the case of potassium chloride. Thus such a study should provide a fairly searching test of the assumed pair potential and indicate whether it is worthwhile to attempt the computation of other transport coefficients. In the present paper we report a molecular dynamics study of molten sodium chloride over a fairly wide range of temperature and density with the aim of studying 'f present address : Applications Software Group, Atlas Computer Laboratory, Chilton, Didcot, Berks.41 THERMODYNAMICS OF MOLTEN SALTS both the thermodynamic properties and the diffusion constants as function of these state parameters. THEORY AND COMPUTATION The method of molecular dynamics essentially involves the setting-up of the coordinates of a relatively small number of atoms or ions, assuming the form of the interactions between these particles, and solving the classical equations of m~tion.~ In the present calculations the following basic assumptions have been made. (a)The interaction potential is pairwise additive. (b) The form of the potential function was taken to be that suggested by Tosi and Fumi lo in a study of solid alkali metal halides : where the parameters used were obtained from the paper of Tosi and Fumi.l0 They are J~/10-19 B/lO1'Jm-l (or+oj)/lO-'O m C/10-79 J m6 D/10-99J mS ++ 0.4225 3.15 2.34 -1.68 -0.80+-0.3380 3.15 2.75 -11.20 -13.90 --0.2535 3.15 3.17 -116.00 -233.00 The first .term in (1) represents the Coulomb interaction, the second the Born-Huggins exponential repulsion with parameters obtained by Tosi and Fumi while the third and fourth terms represent the dipole-dipole and dipole-quadrupole dispersion energies (with para- meters obtained by Mayer ll).(c) The number of particles chosen for the computer experiments was 216, i.e. 108 sodium ions and 108 chloride ions contained in a cube of side L. The reason for choosing this number was computational economy. Some runs with 512 particles indicate that compara- tively little would have been gained in the present investigation by the use of a larger sample.If initially we assume a random set of velocities assigned to the particles such that : 216c miui = 0 (2)i= 1 then after evaluation of the total energy 0and the forces acting on the particles from eqn (3) and (4), the Newtonian equations of motion can be solved numerically for a particular time step giving rise to a new configuration : This time step was chosen sufficiently small to maintain conservation of total energy, the value being 8.0~ s. Although eqn (3) and (4) are straightforward to use, the presence of the Coulomb term in the pair potential causes problems. As is well-known, direct summation of Coulomb potentials in crystal lattices is only slowly convergent for some lattices and divergent for others, and indirect techniques such as Ewald's method, have to be em- ployed.12 This method has been employed in previous MC and MC calculations for ionic systems and of plasma^.^^ l3 The Coulomb energy sum is transformed into two convergent series J.W. E. LEWIS AND K. SINGER and the forces are calculated by taking the gradient of these. In the computations was truncated at Y = L/2 and at lnI2 = 1, with aL = 5.714. After most of the calculations had been completed it was realised that this truncation, while fairly satisfactory for the calculation of the energy sum, is too severe for the calculation of the l5 Fortunately, however, the rather large errors in the individual inter- particle forces have in this case very little effect on the computed average properties : the values for the internal energy, pressure, and diffusion constants reported here are in excellent agreement with those obtained in later MD calculations in which the forces were calculated with an accuracy of -0.3 % (see the Appendix).The cut-off at Y = L/2was also employed for the non-Coulombic forces and energies and a long-range correction was applied to the r6and Y-~energy terms in the usual manner by the replacement of the distant particles by a continuum.16 Initially the system was equilibrated around the desired temperature for about loo0 time intervals (.v10-l1 s). After this, during the production stage of 1800 time intervals, the microcanonical ensemble averages of the potential energy and the virial were calculated as time averages and the instantaneous coordinates and velocities and pair-distance histograms were output to magnetic tape, together with the running averages of energy, virial coefficients, and temperature.During the production run the total energy E was conserved with 0.1 %, From the primary data one obtains the thermodynamic properties internal energy : (E = total energy of the system, No = Avogadro's number, kg = Boltzmann's constant) IN (Y) = -4 1Fij rij. i#j The brackets ( ) indicate ensemble averages. The self-diffusion constants axe calculated either from the mean square displacement of single particles according to the Einstein formula or from the integral of the velocity autocorrelation function D = lim (ui(0) .vi(s))(l -s/t) ds.t-+m s'0 Eqn (7) and (8) are, of course, equivalent and the estimates of the diffusion constants calcu- lated by both formulae are mutually consistent. RESULTS AND DISCUSSION THERMODYNAMIC PROPERTIES In table 1 the molecular dynamics data for internal energy, enthalpy, pressure and self diffusion constants are given. The variation of the pressure along isochores is linear. The slopes of the p against T isochores yield the thermal pressure coefficients Parabolae can be fitted to the interpolated values of p at constant T. One of the roots of each of these parabolae corresponds to the molar volume at zero pressure. THERMODYNAMICS OF MOLTEN SALTS The volumes so obtained, as well as other calculated thermodynamic properties at p = 0, are listed together with experimental data in table 2.The isothermal compressibility and thermal expansion coefficients have been evaluated in two ways. (i) The p against V-isotherms were differentiated near p = 0 to yield the isothermal compressibility KT = -i/v(avjap)T. (10) The thermal expansion coefficients up = i/~(av/a~), (1 1) were obtained by fitting straight lines to the density against temperature curves at p = 0. (ii) The relationship was used. Parabolae were fitted to the values of pv along isotherms ; the enthalpies are easily calculated as U+pY and are found to lie on linear isochores. From these '0 2 4 6 8 10 12 r/10-8 cm FIG.1.-Pair distribution functions for sodium chloride just above the melting point : -, g+-, ---, 9+-+ b+++9-42.1qI ,q '0 2 4 6 8 10 12 r/1W8cm FIG.2.-Pair distribution functions for sodium chloride just above the melting point : -, g++; Y 9--. J. W. E. LEWIS AND K. SINGER isotherms H(V) are constructed; parabolae can be fitted to these, which lead to (dH/aV)and through eqn (12) to KT. The thermal expansion coefficient is then obtained from the chain relation This method is inferior for the calculation of a, due to accumulation of errors when using eqn (13). The values of a, shown in table 2 have been evaluated by method (i). The values obtained by the two routes agree within about 15 %. TABLE1 D, 10-9 D-110-9 V/10-6m3 T/K Plkbar UjkJ mol-1 H/kJ mo1-I m2 s-l m2 s-1 37.51 1056 1.23 -694.15 -689.54 7.4 7.1 1375 5.16 -678.74 -659.38 11.8 10.3 1622 8.46 -665.66 -633.93 14.9 13.1 1762 9.96 -659.52 -622.16 17.8 16.7 2060 14.00 -644.74 -592.23 24.2 18.5 39.47 1123 0.45 -688.21 -686.43 8.8 8.1 1377 3.03 -675.17 -663.21 13.3 12.5 1600 6.31 -661.41 -636.50 17.5 15.7 2056 10.12 -643.17 -603.23 25.0 23.6 41.35 1073 -1.47 -688.69 -694.77 8.2 7.9 1381 1.53 -672.81 -666.48 14.5 12.3 1723 4.34 -656.33 -638.34 17.4 16.8 2143 8.51 -637.20 -602.01 26.4 25.4 43.29 1262 -1.13 -676.41 -681.29 12.7 12.0 1330 -0.39 -672.95 -674.65 13.9 12.6 1406 0.19 -669.29 -668.45 15.5 14.2 1485 1.01 -665.31 -660.93 16.0 15.3 1562 1.44 -661.45 -655.22 18.9 18.5 47.35 1459 -0.09 -661.96 -666.65 17.4 15.1 1566 -0.63 -657.01 -659.99 17.7 16.6 1678 0.17 -651.72 -650.92 22.2 20.3 1784 0.98 -646.72 -642.08 26.5 23.6 TABLE2 V/m3 mol-1 cr~/lO-~K-' K~110-6bar-' pv/bar K-1 Cv/J K-i mol-l TlK calc.exp. calc. exp. calc. exp. calc. exp. calc. exp. lo00 38.15 3.04 31.1 11.9 49.1 48.9 1073(m.p.) 39.14 37.70 3.12 3.07 32.4 28.7 11.1 10.7 48.9 48.6 1100 39.51 38.01 3.15 3.09 33.0 30.1 10.8 10.3 48.8 47.6 1200 40.88 39.23 3.26 3.19 35.3 35.5 9.8 9.0 48.6 47.2 1300 42.27 40.54 3.37 3.30 38.3 41.7 8.8 7.8 48.3 1400 43.70 41.90 3.49 3.41 42.3 8.O 48.1 1500 45.17 43.40 3.60 3.53 47.7 7.2 47.8 1600 46.71 44.96 3.73 3.66 55.8 6.5 47.7 1738(b.p .) 49.00 46.88 3.91 3.58 74.6 5.7 47.4 THERMODYNAMICS OF MOLTEN SALTS The molar heat capacity at constant volume C,,is determined by the slopes of the U(T),isochores.Cpcan be similarly obtained from the gradients of the H(T),plots. The variation of Cpwith pressure is shown in table 3. From the values of Cpand Cv and KT,it is possible to re-compute a, by means of VTaic,-cy = -. KT The agreement with the previously calculated values is satisfactory. The calculated and experimental thermodynamic properties at zero pressure are compared in table 2. The magnitude of the discrepancies is, for U(T)-k0.5%, for V(T)-+5 %, for Cv-0.3 %, for C,-2 %, for up,KT, BB-&5-10 %. The agreement is surprisingly good, particularly if it is remembered that the parameters in the pair potential were obtained by fitting the properties of the crystal at 298 K.1° The most obvious discrepancy is the over-estimation of the volume by -5 %.The statistical error in the calculated pressures is -kO.3 kbar, which is equivalent to an error of k0.3 cm3 in the molar volume at 1100 K. This clearly does not account for the systematic discrepancy of +2-3 cm3. The results are consistent with those obtained by Monte Carlo calculations for liquid KCl and for other alkali metal halides, and must be considered as a general feature of the Tosi-Fumi pair potentials. The calculated densities conform to : p(T, p = O)/g ~m-~ = 1.993 -4.662 x T/K (14) for the experimental values p (T,p = 0) = 2.061 -4.759 x T/K.I7 For H(T,p = 0) one obtains the equation (-761.5+0.066 T/K) kJ mol-' which is to be compared with the -764.5+0.067 T/K for the experimental values 18* 19* lo between the melting point and 1300 K; the agreement is better than 0.5 % The agreement between observed and calculated heat capacities is very good.The linear dependence of the enthalpy on temperature, which is characteristic of molten alkali metal halides, applies to the calculated values over the entire temperature range. The value of CVis very close to R = 3NkB,showing that the law of Du-Long and Petit holds. Since this law holds for systems of classical harmonic oscillators, one is led to the conclusion that the ions move on average in parabolic potential wells. The discrepancy between the calculated and experimental values 179 2o of the pv,KT, up is probably no greater than would be expected from the statistical errors attached to the MD data, combined with the margins of error in the experimental measurements. The thermal expansion coefficients agree with the experimental values within 3 % (except at the highest temperature).For the isothermal compres- sibility the discrepancy is of the order of 10 % and there is an indication that the increase with temperature is too small. The calculated thermal pressure coefficients are by 5-10 % larger than the observed values. In summary, one can conclude that the Born-Huggins-Mayer-Tosi-Fumi potential simulates the thermodynamic properties of liquid NaCl quite well. It must be remembered, however, that, owing to the absence of experimental data, the com- parisons are restricted to the behaviour at low pressure.EQUATIONS OF STATE AND FOR THE FREE ENERGY Over the range covered by the MD data, the pressure can be expressed as plkbar = (42.2V-2-0.898 V-l+ 0.006 44)T/K- 287.0 Y-1-4.480 (Y in cm3). (1 5) Although this equation is valid within 50.3 kbar, the values of kVand KTderived from it are not in close agreement with those reported in table 3. The discrepancy arises J. W. E. LEWIS AND K. SINGER from the additional curve-fitting step required to obtain the coefficients in eqn (15). The internal energy can be expressed as U/kJ mol-l = 48.44 x T/K-2655V-l/~m-~-675.5. (16) The coefficient of Tis the average value of Cv; its small variation with V is neglected.TABLE3 ~~ plkbar calc. exp. 0.0 65.9 67.0 1.o 65.1 2.0 64.4 3.O 63.8 4.0 63.2 5.0 62.8 6.0 62.4 7.0 62.1 8.O 61.8 9.0 61.5 10.0 61.2 Although eqn (1 5)and (1 6) have a reasonable form, they are not derivable from the same expression for the Helmholtz free energy A( V, T),and are thus thermodynamic- ally not compatible. This can be seen by calculating (aU/aY), (i) from eqn (16), (ii) by means of (aU/aY), = T(i?p/~V)~-p,from (15). Thermodynamically compatible equations can be obtained in the form U = C,T+b,V2+b2V+b3 (17) p = (c, V-2+c2V-1+~3)T+a,Y+a, = P,T+alV+a,. (18) The coefficients cl, c,, c3 and Cvhave the same values as in (1 5) and (1 6) respectively The remaining coefficients can be determined by least squares curve fitting, but the primary data are not sufficiently accurate to ensure the compatibility of the equations so obtained.Instead the coefficients al and a2 were obtained by the curve fitting, and bl, b2 from the conditions of compatibility i.e. b, = -a1/2, b, = -a2 (19) and b, finally by mean difference from (17). The resulting expression U = (48.44 x T/K-7.989 x V2/cm6+ 1.8106V/cm3-801.17) kJ mol-l (20) fits the data within +O.l %, while the equation p = (Pv(V)T+0.1598 V/cm3-18.1 1) kbar (21) holds with an accuracy of better than k0.5 kbar. Integration of (18) with respect to volume leads to A(V, T) = -sp dV = (c, Y-l -c2 In V -c3V)T-+a,Y2-a2V -const(T). (22) Substitution of (20) into A(V, T) = -T dT’ leads to A( Y,T) = -CVTIn T + (b, V2+b, V+b3)+T const’(Y). THERMODYNAMICS OF MOLTEN SALTS Eqn (22) and (23) are compatible if const(T) = C,Tln T-b3 and const'(V) = c1Y-l-c2In V-c3V, so that finally A( V, T) = -C,T In T+ (c,V-' -c2In V- c3 V)T-$al V2-a2VSconst.(24) Eqn (24) could serve as a starting point for the derivation of equations of state for other liquid alkali metal halides by corresponding states arguments.2 TABLE4 A B C cations 1.0207~ -47.92 1537.6 anions 9.047~ -39.05 1239.5 D = (ATV+B+C/T') x cm2s-' TABLE5 ~+/10-5crn2 s-~ ~-/10-5cm2 s-l T/K v/cm3 mol-' calc. exp. calc. exp. 1073 39.14 8.1 7.3 7.5 5.8 1100 39.51 8.5 8.0 7.9 6.3 1200 40.88 10.2 10.5 9.4 8.4 1300 42.27 12.2 13.2 11.3 10.7 1400 43.70 14.6 16.1 13.5 13.1 1500 45.17 17.4 19.1 16.1 15.7 1600 46.71 20.6 19.1 1738 49.00 26.0 24.0 RADIAL DISTRIBUTION FUNCTIONS The pair correlation functions g++(r),g+-(r) and g--(r), as well as the mean or total correlation function gm(r),have been calculated.Some of these, obtained for systems with 512 ions, are shown in fig. 1 and 2 and the characteristic features at different VT points are summarised in table 6. The results are similar to those obtained in previous computer simulations 4-6* : the maxima of g++and g--coincide and fall on the minima of g+-and vice versa. The first peak of g+-(and of gm)occurs at a distance smaller by 0.15-0.2A than obtained by diffraction experiments. There is some over- lap of the g++ and g--curves and the first peak of gm,showing that the like ion pairs contribute >5 % to the first shell coordination number; this also has the conse- quence that the first minimum of gm(t)occurs at smaller distances than that of g+-(r), and the "number of nearest neighbours " estimated by integration of g,(r) to its first minimum is much lower than the coordination number obtained by integration of g+-(r) to its first minimum.These and other features have been more fully discussed el~ewhere.~' Table 6 shows that (rmax)+-is fairly independent of temperature and 22 volume. This is in contrast with the results for liquid KCl where there is a small but definite decrease of (rmaX)+-with increasing v01ume.~ DYNAMIC PROPERTIES The most readily computed quantities are the diffusion constant D and the velocity autocorrelation function.The diffusion constant can be computed from MD data either by the Einstein formula J. W. E. LEWIS AND K. SINGER from the linear increase of the mean square displacement of single particles, or by integration of the velocity autocorrelation function D = 3 lim (~~(0)-ui(s))( 1-sit) ds. t-rw s'0 The two methods are, of course, equivalent and the estimates of D obtained by them are mutually consistent. Following Verlet and Leve~que,~~ an equation of the form D = ATV2+B+CV-' (25) has been fitted to the primary data, and found to be valid within 5-10 %, which is the same magnitude as the computational error attached to the data. The coefficients are given in table 4.Eqn (25) together with eqn (14) have been used to calculate values of D(T, V) at zero pressure. These values are compared with experimental data 24 in table 5. In view of the margins of error OF the quantities so compared, the agreement is not unsatisfactory. The calculated values particularly for the anion are by up to 20 % too high near the melting point ; at higher temperatures the agreement is better because the calculated D values vary more slowly with Tthan the experimental ones. TABLE6 molar vol- cm3ume/ mol-' TIK rmax g++ rmin h rmax 9-rmin h * rmax g+-rmin h n rmax gr? rmin h n 37.51 1056 4.0 6.2 0.89 14.6 4.0 6.2 0.94 14.5 2.7 4.2 3.76 5.2 2.7 3.3 1.88 3.8 1375 4.0 6.1 0.84 14.2 4.0 6.0 0.88 13.8 2.6 4.1 3.44 5.0 2.6 3.3 1.72 4.0 1622 4.0 6.2 0.80 14.6 4.0 6.1 0.86 14.2 2.6 4.3 3.24 5.4 2.6 3.2 1.63 3.5 1761 4.0 6.2 0.78 14.6 4.0 6.2 0.84 14.6 2.6 4.2 3.17 5.2 2.6 3.3 1.59 4.0 2061 4.0 6.0 0.75 13.3 4.0 6.1 0.80 14.6 2.6 4.2 3.03 5.3 2.6 3.2 1.52 3.5 39.47 1123 4.1 6.1 0.89 13.3 4.1 6.1 0.90 13.3 2.6 4.1 3.78 4.8 2.6 3.3 1.89 3.6 1660 4.1 6.3 0.79 14.5 4.1 6.1 0.82 13.7 2.6 4.2 3.25 5.0 2.6 3.2 1.63 3.2 2056 4.1 6.3 0.75 14.6 4.0 6.1 0.78 13.3 2.5 4.1 3.05 4.7 2.5 3.2 1.53 3.2 41.35 1073 4.2 6.2 0.87 13.2 4.2 6.2 0.91 13.2 2.6 4.0 3.91 4.5 2.6 3.3 1.95 3.6 1381 4.2 6.2 0.82 13.2 4.2 6.2 0.86 13.2 2.6 4.2 3.58 4.8 2.6 3.2 1.79 3.3 1723 4.2 6.1 0.76 12.8 4.2 6.2 0.83 13.2 2.6 4.2 3.29 4.8 2.6 3.2 1.65 3.3 2143 4.2 6.1 0.74 12.8 4.0 6.0 0.76 12.8 2.6 4.1 3.03 4.6 2.6 3.2 1.52 3.2 43.29 1265 4.2 6.3 0.88 12.7 4.2 6.3 0.91 12.7 2.8 4.1 3.80 4.5 2.8 3.3 1.90 3.0 1329 4.1 6.1 0.87 12.7 4.1 6.1 0.90 12.7 2.7 4.0 3.76 4.5 2.7 3.2 1.88 3.0 1406 4.1 6.1 0.87 12.7 4.1 6.1 0.88 12.7 2.7 4.0 3.68 4.5 2.7 3.2 1.84 3.0 1485 4.1 6.1 0.84 12.7 4.1 6.1 0.84 12.7 2.7 4.0 3.58 4.5 2.7 3.2 1.79 3.0 1561 4.1 6.1 0.83 12.3 4.1 6.1 0.87 12.7 2.7 4.0 3.51 4.4 2.7 3.2 1.76 3.0 47.35 1457 4.3 6.4 0.84 12.9 4.3 6.4 0.86 12.9 2.6 4.1 3.80 4.3 2.6 3.3 1.90 3.0 1566 4.3 6.4 0.82 12.9 4.3 6.4 0.85 12.9 2.5 4.1 3.69 4.2 2.5 3.3 1.84 3.0 1678 4.4 6.4 0.80 12.9 4.3 6.4 0.84 12.9 2.5 4.1 3.62 4.2 2.5 3.3 1.82 3.0 1782 4.3 6.4 0.79 12.9 4.3 6.3 0.84 12.0 2.5 4.1 3.58 4.2 2.5 3.3 1.79 3.0 r,,, = position of first peak (A, ; rmin= position of first minimum (A) ; h = height of first peak ; n = coordination number = s:"4nrzg(r)dr.Velocity autocorrelation functions have been obtained for a number of V, T points. Fig. 3 shows normalised velocity autocorrelation functions, defined by P(t>= Q<t>/<lui(0)12> ; Q(t>= <vi(O)*ui(t)> (26) for sodium and chloride ions in a melt at the V,Tpoint (43.29 cm3, 1406 K).Some generalisations concerning the V,Tdependenceof the diffusion constants and that of the velocity autocorrelation functions can be made. It should be borne in THERMODYNAMICS OF MOLTEN SALTS mind, however, that although the MD data cover a fairly large temperature range, this is small compared with the entire liquid range and the system is always far from the critical region. L time (intervals of 8 x 10-l5s) FIG.3.-Velocity autocorrelation functions for the sodium and chloride ions in molten sodium chloride at 43.29cm3 and 1406K. The diffusion constant can be written in the form CQ D = +(1vJ2) 1 P(t)dt = 1; P(t)dt. 0 Expressing D by the Einstein relation in terms of the friction coefficient D = k,T/my, one obtains l/y = T = Jrm P(t)dt.(28)0 Thus if y is independent of temperature then at constant density, D will be a linear function of T. This is seen to be nearly the case [eqn (25)]. The density dependence, on the other hand lies entirely in y. The temperature dependence of D at constant density can therefore be ascribed to the change of the mean square velocity, whereas the integral over the normalised velocity autocorrelation function is almost temperature independent. As the temperature increases, the negative portion of P(t) is diminished and the crossing point from positive to negative values shifts to large times. Since the integral remains virtually constant, one may conclude that the more rapid loss of correlation with time occurs in such a manner that the decrease in area of the negative portion is compensated by an equal decrease in the area of the positive portion.These features are similar to those observed by Levesque and Verlet 23 for the velocity autocorrelation function of the Lennard-Jones liquid. Detailed examination shows, however, that in the case of the ionic liquid there are oscillations of the order of 2 x s [leading to a second maximum and minimum in P(t)]. It seems likely that these are remnants of a strongly damped oscillation, which, in the absence of damping by short range forces, is a dominant feature in the velocity autocorrelation function of plasmas. J. W. E. LEWIS AND K. SINGER In addition to the computation of thermodynamic and transport properties, computer simulation can also provide an intuitive insight into the structure and dynamics of molten salts.To this end computer generated films illustrating the motion of the ions have been made. A fairly general program with full hidden line removal for use in three dimensional systems has been written with graphical output on a SD4020 microfilm recorder, making it possible to construct film sequences. Still frames from such sequences are shown in fig. 4 (a), (b) and (c). The presence of (c) FIG.4.-Still frames of molten sodium chloride simulation display at 100 time interval separations, Ionic radii of 0.75 and 1.00x lo-* cm for Na+ and C1-respectively for the purposes of clarity. large voids within the liquid alkali metal halides has been disc~ssed.~.Similar22 voids have also been observed by Fehder et aZ.25in a two-dimensional MD simulation of Lennard-Jones disks. The moving sequences show the dynamic behaviour of these voids : they appear to have a fairly short life time (10-l2 s) and they tend to collapse rather than being jumped into by a single ion. The word "void " in the present context is preferred, since "hole "would suggest the presence of holes in the sequence of cell-hole models, which have been shown to be invalid for liquid alkali metal halides.6 The presence of voids may be simply a consequence of the small coordination number (-5 nearest neighbours as compared with -10 in Lennard-Jones liquids). It is less easy to understand why the voids are inaccessible to single ions since they seem to be THERMODYNAMICS OF MOLTEN SALTS large enough to accommodate the diameter of an ion and since the electrostatic poten- tial inside them must be either positive or negative. The last argument, however, is probably fallacious : because of the long range of the electrostatic forces, it is insuffi- cient to consider the potential inside the void only and to disregard the changes elsewhere, resulting from the displacement of a single ion.It would be instructive to map the electrostatic field of relevant configurations in three dimensions. Another feature, which is possibly statistically insignificant, can be seen from the stills : tetrahedral arrangements of ions of opposite charge around a central ion are not infrequent.CONCLUSIONS It is initially surprising that a potential of the form given by equation (1) and obtained from solid state data should give such good results for the thermodynamic properties of sodium chloride over a large part of the liquid range. They are by no means perfect, the zero pressure volumes being the most seriously in error, although the other thermodynamic properties appear to agree quite well with experimental values. Experimental measurements of the self-diffusion constants of molten sodium chloride have been confined to zero pressure and little is known about th3 temperature and density dependence of the diffusion process. The reported calculations suggest that the temperature dependence of the self diffusion constant is approximately linear, but as might be expected the density dependance is complicated.However, it appears to be approximately represented by eqn (26). The results of this study indicate that the pair potential (1) is a fairly good represent- ation of the effective pair interaction in molten sodium chloride over a considerable range of temperature, and this potential could be used as a starting point for further studies. Finally the method of molecular dynamics applied to molten salts besides yielding values of the thermodynamic properties of the system and values of self diffusion constant is further capable of yielding information about both dynamic and time averaged microscopic structure of such systems. APPENDIX EFFECT OF THE SEVERE TRUNCATION OF THE EWALDSUMS The magnitude of the error in the MD data calculated on the basis of the severe trunca- tions Y = L/2 (in @I) rzLaX= 1 (in @IT) used in the present work is examined by comparisonof values for the energy, pressure and diffusion constants obtained by MD calculations TABLE7 U/kJ mol-I plkbar D+/IO-scm2s-1 ~-/10-5cm*s-1 V/cm3 T/K this work a.E s.* this work a.E.s.* this work a.E.s.* this work a.E.s.38.4 1165 -687.0 -684.1 ' 9.7 9.8 a 8.8 8.5' 39.5 1225 -682.8 -679.9' 10.5 10.8 a 9.6 9.6a 41.7 1341 -674.6 -673.6' 12.7 14.1 a 11.8 12.2a 39.1 1091 -689.7 -691.8 0.24 0.29 , 8.4 8.6b 7.7 7.7 39.1 1281 -680.5 -679.7 2.34 2.59 11.4 10.6b 10.4 9.5 47.35 1817 -645.5 -644.5 0.91 1.65 26.1 23.2 24.0 21.4b * aEs = accurate Ewald summation ; a ref.(8) ; b Gosling and Singer, unpublished data, truncations in the Ewald sums : this work : rc = L/2, niLx = 1 : a rc = L, nL = 6; b rc = L/2, nt, = 27. J. W. E. LEWIS AND K. SINGER involving a much more accurate evaluation of the electrostatic forces. Such calculations for a few V, T points for liquid NaCl have been carried by Lantelme et aL8 with the trunc- ation Y = L, niax= 6 (aL = 3.5) and with Y = L/2, n&x = 27 (aL= 5.6).25 In these computations the electrostatic forces are calculated with an accuracy of 20.3 %. The figures in table 7 show that the discrepancies are for the most part within the margin of statistical error.Only at the highest temperature is this margin exceeded for the pressure (-0.7 kbar) and for the diffusion constant (-15 %). It would appear that the error in the calculation of the electrostatic forces, even if large for single particles, has a very small effect on the calculated averages. It is not suggested, however, that this favourable situation would necessarily apply in the MD simulation of other molten salts. One of us (J. W. E. L.) is grateful to the S.R.C. for a Research Assistantship. We also thank the staff of the Atlas Computer Laboratory for a generous allocation of computer time and for assistance in running the programs on the IBM 360/195 at the Rutherford Laboratory. K. Singer and I. R. McDonald, Quart. Rev. Chem. Soc., 1970, 24, 238.I. R.McDonald and K. Singer, Ann. Rep. Chem. Soc., 1970, 45. J. A. Barker and D. Henderson, Ann. Rev. Phys. Chem., 1972, 23,439. L. V. Woodcock and K. Singer, Trans. Faraday Soc., 1971, 67, 12. D. J. Adams and I. R. McDonald, J. Phys. C, in press. J. W. E. Lewis, L. V. Woodcock and K. Singer, to be published. 6a L. V. Woodcock, Chem. Phys. Letters, 1971, 10, 257. 'A. Rahman, R. H. Fowler and A. H. Narten, J. Gem. Phys., 1972,57, 3010. 13 F. Lantelme, P. Turq, B. Quentrec and J. W. E. Lewis, to be published. A. Rahman, Phys. Rev. A, 1964, 136,405. M. P. Tosi and F. G. Fumi, J. Phys. Chem. Solids, 1964, 25, 31. J. E. Mayer, J. Chem. Phys., 1933, 1, 270. "P. P. Ewald, Ann. Phys., 1921, 21, 1087. l3 A. A. Barker, Austral.J. Phys., 1965,18,119 ; S. G. Brush, H. L. Sahlin and E. Teller, J. Chem. Phys., 1966, 45, 2102. l4 M. Dixon and M. J. Sangster, unpublished work. E. M. Gosling and K. Singer, unpublished work. l6 W. W. Wood in Physics of Simple Liquids, ed. H. N. V. Temperley, J. S. Rowlinson and G. S. Rushbrook (North Holland, Amsterdam, 1968).A. D. Kirchenbaum, J. A. Cahill, P. J. McGonigal and A. V. Grosse, J. Inorg. Nuclear Chem., 1962, 24, 1287. l8 K. K. Kelley, US.Bur. Mines Bull., 1960, 584. l9 A. S. Dworkin and M. A. Bredig, J. Phys. Chem., 1960, 64, 260. 2o J. O'M. Bockris and N. E. Richards, Proc. Roy. Soc. A, 1957, 241,44. 21 M. Blander, Adu. Chem. Phys., 1967, 11, 83. 22 L. V. Woodcock, Proc. Roy. Soc. A, 1972, 328, 83, 23 D. Levesque and L. Verlet, Phys. Rev. A, 1970, 2,2514. 24 J. O'M. Bockris, S, R. Richards and L. Nanis, J. Phys. Chem., 1965, 69, 627. 25 E. M. Gosling and K. Singer, unpublished work.