Inherent structure of a molten salt
Randall A. LaViolette and Joanne L. Budzien
Idaho National Engineering and Environmental Laboratory
P.O. Box 1625, Idaho Falls, ID 834152208
Frank H. Stillinger
Bell Laboratories
600 Mountain Avenue, Murray Hill, NJ 079740636
Abstract
We calculated the inherent structure of a model melt of zinc (II) bromide over a wide range of densities. Stable, metastable and unstable branches were obtained for the zero temperature pressurevolume isotherm of the inherent structure. The pressurevolume isotherm, the void distribution and the structure factor were used to identify the spinodal, independent of any model equation of state.
I. Introduction
The concept of the inherent structure of a liquid was introduced in order to divide the structure of a liquid into the part due to thermal fluctuations, and the part due to the packing in the absence of fluctuations. One can calculate the inherent structure by extracting, from the simulated liquid, a sequence of statistically independent instantaneous configurations, and subjecting each to an instantaneous steepestdescent quench to the nearest local minimum on the potential energy hypersurface. The inherent structure is the ensemble of the resulting nonequilibrium, zerotemperature, and usually amorphous, configurations, often represented by the radial pair distribution function (PDF) calculated for the ensemble. The inherent structure of strongly associated liquids (e.g., water) substantially reproduces the structure of the equilibrium liquid, i.e., the peaks and valleys of the PDF differ in width and magnitude but line up with those of the liquid, and no new peaks or valleys are found in the inherent structure. For every other type of liquid studied to date (e.g., one or twocomponent van der Waals atoms^{,}, metals and metal alloys, silicon, sulfur, and purely theoretical systems^{,}), the inherent structure contains more short and mediumrange order than the corresponding equilibrium liquid. In particular, the second peak of the liquid PDF is usually split in the inherent structure PDF, although the first peak is usually changed only in its width and magnitude.
The inherent structure of ionic melts has not been examined as far as we know. Therefore we studied the inherent structure of a model ionic melt. We chose a system with a substantial size and charge mismatch between the cation and the anion, rather than the canonical 1:1 system with identical radii, because the former are usually better glass formers, and because the latter corresponds closely only to potassium fluoride. To be specific, we chose a simple model of the interactions between the ions that corresponds to zinc (II) bromide (ZnBr_{2}). The ZnCl_{2} melt has been much more thoroughly studied than ZnBr_{2}, both theoretically and experimentally^{,}, is known to be a good (but relatively fragile with respect to transport^{12c,d,e}) glass former with a low glasstransition temperature^{,} and good empirical force fields have been developed for that system. Nevertheless, the ZnBr_{2} system has more immediate interest for us because of its utility as a gammaray shielding material and issues arising from its degradation and disposal, although we shall have nothing more to say about such applications here.
We conclude this Introduction with a plan of this paper. Section II. describes the details of the model and the calculations. Section III. contains the results of the calculations, including the pair distribution functions and the pressurevolume isotherm (at zero absolute temperature) of the inherent structure of the ionic melt. We conclude this paper in Section IV. with a discussion of the results.
II. Methodology
We follow the approach employed previously for simple monatomic liquids, combining traditional molecular dynamics with instantaneous quenching to obtain the pressurevolume isotherm (or "equation of state") for the system. Molecular dynamics has been employed to study ionic melts for over two decades^{,,}. The primary input to a molecular dynamics calculation is a model for the forces, or the energy from which the forces may be deduced. Our model for the pair potential energy between the ions is a slight modification of the models by Wilson and Madden and is listed below in Table 1. Although we expect only a qualitative correspondence between this model and the ZnBr_{2}, the model contains the essential features that we seek: a doublycharged cation, a singlycharged anion, and a ratio of radii between 1:2 and 1:3.
Table 1: Ion pair interaction potentials
The potential energy U is given in kcal/mol and the pair distances r in Å.
Zn^{2+}Zn^{2+} 

Zn^{2+}Br^{} 

Br^{}Br^{} 
Our constantenergy molecular dynamics simulation of the melt employed 486 Zn^{2+} and 972 Br^{} ions in a cubic box with periodic boundary conditions, for twelve densities, from 4.174 to 1.941 g/cm^{3}. For reference, the equilibrium solid at room temperature is 4.2 g/cm^{3}. The Coulombic interactions were calculated using the method of Ewald sums. The cutoff parameters for the sum were determined using the method of Karasawa with an accuracy of 0.001 kcal/mol/atom. The dynamics were advanced with the Verlet leapfrog algorithm with a time step of 1 fs. For each density, the velocities were adjusted at the beginning of the run to keep the temperature of the melt at 2000 ºC, on average. Although this temperature is much higher than that of the experimental boiling temperature (650 ºC) at atmospheric pressure, even at atmospheric pressure the simulated melt showed no signs of phase separation. The high temperature is employed in order to make the sampling of configuration space (and of the "basins" that partition it) efficient. The density series was generated by starting the melt simulation from the ending configuration of the next highest density. The box was expanded to obtain the desired density. The system was then subjected to about 30 ps of microcanonical dynamics to allow the melt to equilibrate. Subsequently, twenty instantaneous configurations were extracted at 12 ps intervals in order to represent the inherent structure. The interval for each density was constant; however, some of the higher densities had smaller intervals (1 ps) than the lower densities (2 ps).
The minimizations (instantaneous quenches) were performed by employing the extracted configurations of the melt as starting points for a constantvolume steepestdescent trajectory to the nearest potential energy minimum of the ions. A literal implementation of steepestdescent algorithm requires orders of magnitude more time than more sophisticated methods that essentially reproduce the steepestdescent path. A combination of steepestdescent and quasiNewton procedures has been used successfully in previous work, but in this study, we found that the conjugate gradient technique was even more efficient and just as accurate. We employed the PolackRibiere algorithm with a line search precision of 0.1 so that the final gradient value is less than 10% of the original gradient value. For energy convergence, the energy difference between successive steps must be less than 0.001 kcal/mol; the rootmeansquared (r.m.s.) displacement of the atoms, 0.003 Å; and the r.m.s. force, 0.100 kcal/mol/Å. The minimization was refined by requenching using more demanding convergence criteria (energy difference < 10^{4} kcal/mol, r.m.s. displacement < 10^{5} Å, maximum force < 5 ´ 10^{3} kcal/mol/Å, maximum displacement < 5 ´ 10^{5} Å) than that used initially in the quench.
We calculated the average pressure, energy, ion diameter, relative volume accessible to Br^{}, total static structure factor S(k), and the atomic pair distribution functions (PDF) for the twenty quenched configurations that we employed to represent the inherent structure for each density. Neither the pressure nor the energy contain contributions from the kinetic energy, which is zero for the inherent structure. We calculated the PDFs in the usual way by binning all pair distances less than half the box length, after the application of periodic boundary conditions, with a bin width of 0.05Å. Although the PDF is formally the Fourier transform of S(k) – 1, the structure factor was computed independently^{ }as the average of exp(ik •r) on a grid in kspace. The average Br^{} diameter was determined by locating the first maximum of the Br^{}Br^{} PDF of the inherent structure. Because almost no contacts between the strongly repulsive zinc ions were found, the average Zn^{2+} diameter was determined indirectly by subtracting the Br^{} diameter from twice the diameter of the Zn^{2+}Br^{} pair (also determined from the location of the first maximum of the PDF),. The relative available volume was determined on a grid system. The value of the length of the box was divided by 0.415 Å (about half the zinc ion radius) and rounded the result to the nearest integer to obtain the spacing. A test sphere of the same radius as the ion of interest was placed at each grid point. If there were no overlaps with ions actually present in the cell, the point was counted as a available volume point. The ratio of the number of such available volume points to the total number of points is the available volume fraction for the ion.
III. Results
A. Mechanical and thermodynamic properties
We simulated the melt, and calculated the inherent structure, for twelve densities, from 4.174 to 1.941 g/cm^{3}, or in terms of specific volume, from 0.2396 to 0.5152 cm^{3}/g. The minimum energy of the inherent structure occurs between 0.2862 and 0.3039 cm^{3}/g (3.494 to 3.291 g/cm^{3}). Zero pressure occurs between 0.3096 and 0.3220 cm^{3}/g (3.291 to 3.106 g/cm^{3}) and the minimum pressure (slightly below –0.4 GPa) occurs between 0.3220 and 0.3621 cm3/g (3.106 to 2.762 g/cm^{3}). The former interval contains the boundary between the thermodynamically stable (positive pressure and compressibility) and the metastable (negative pressure but positive compressibility) branches of the inherent structure equationofstate; the latter interval contains the "spinodal", i.e., the boundary between the metastable and the unstable (negative pressure and compressibility) branches. The ionic diameters are essentially constant throughout the density range: the Zn^{2+} diameter varies between 2.0 and 2.2 Å; the Br^{} diameter varies between 3.6 and 3.8 Å. For comparison, the empirical mean diameter in a crystal is 1.48 Å and 1.66 Å for the four and sixcoordinated Zn^{2+}, respectively, and 3.92 Å for the Br^{}. For another comparison, the model gives a Zn^{2+}Br^{} distance of 2.56 Å for the smallest neutral gas phase configuration of a zinc ion between two bromides (at the antipodes of the zinc ion); as is usual for the fused salts, this distance is shorter than that (2.7  2.8 Å) predicted from the condensed phase diameters. The available volume (or void) around the the Zn^{2+} ion increases smoothly with increasing specific volume for all branches, but the available volume around the Br^{} ion begins to increase significantly (and abruptly) only in the metastable branch, corresponding to the growth of cavities in the metastable and unstable branches. Figure 1 displays both the pressure vs. specific volume, and the volume available to the Br^{} vs. specific volume.
Figure 1. Pressurevolume isotherm of, and available volume to the bromide, in the inherent structure of molten ZnBr_{2}
The specific volume is given in cm^{3}/g, the pressure is in GPa. Closed circles mark the pressure of the inherent structure. The solid curve marks the % volume available to the bromide in the inherent structure. Block arrows point to the boundaries of the stable and unstable regions, respectively, of the inherent structure. The doublebar doublearrow marks the metastable region; the smaller solid doublearrow marks the region that contains the spinodal (its width is our uncertainty). Open circles mark the pressure of the calculated melt. The dashed line marks zero pressure or per cent for reference.
For comparison we include the pressure calculated for the model melt. Although we have not found an experimental pressurevolume isotherm for molten ZnBr_{2}, we have found an empirical fit of the variation of the density of the zeropressure melt with temperature, for 400 ºC < T < 620 ºC. Extrapolating the fit well beyond its intended range gives a metastable (superheated) melt density of 1.93 g/cm^{3} (specific volume 0.52 cm^{3}/g) at 2000 ºC, which coincides closely with our estimate as shown in Figure 1. This suggests that the model interactions between the ions, while simple and not intended to provide the acme of realism, constitute an adequate approximation.
B. Structure
We calculated the static structure factor S(k) first in order to see if we could obtain independent verification that the minimum in the inherent structures pressurevolume isotherm actually corresponds to a spinodal. For densities above the minimum, S(k) should approach zero for k® 0, because S(0) is proportional to the product T c, the isothermal compressibility c is finite, and the temperature is absolute zero for the inherent structure. However, at the spinodal, S(k) should begin to diverge for small k, even at zero temperature, because of the rapid divergence of the compressibility there. Figure 2 shows the static structure factor for some representative densities. The divergence of S(k) for small k is apparent for densities corresponding to the unstable branch.
Figure 2. The total structure factor for the inherent structure of liquid ZnBr_{2}
The wavenumber k is given in Å^{1}. The specific volumes for each structure factor are listed (in cm^{3}/g) in the legend. The curves are splines through the data.
For the melt, the PDFs for all three interactions show only minor differences as function of density (Figure 3). The primary differences are the peak shift so that the neighbors are slightly closer together at higher density. The second peak of the BrBr PDF is split at high densities, and gradually coalesces at lower densities.
Figure 3. Pair distribution functions for the melt
The top, middle and bottom panels correspond to the zinczinc, bromidebromide, and zincbromide PDFs, respectively. The pair distance r is given in Angstroms. The legend labels the curves according to density in g/cm^{3} (0.24 and 0.52 in cm^{3}/g respectively).
Figure 4. Zn^{2+}Zn^{2+} inherent structure PDF compared with the melt.
The density is 3.707 g/cm^{3} (specific volume 0.27 cm^{3}/g). The distance r is in Angstroms. The first peak is beginning to split. The inset is the region around the first peak of the inherent structure for various specific volumes.
For all densities, the PDFs of the inherent structure show (as was indicated implicitly by S(k)) that the inherent structure is amorphous. The Zn^{2+}Zn^{2+} PDFs of the melt and the inherent structure (Figure 4 ) are essentially the same at high densities, except for the enhancement of the peaks and valleys of the latter. The fact that there is no peak in the vicinity of the diameter of the zinc ion shows that there are almost no pairs of zinc ions that are in contact. The location of the first peak occurs at about the right distance (5.6 Å, compared with 5.45.6 Å expected from the mean ion diameter) for the Zn^{2+}Zn^{2+} pair to be separated by a bromide at low density, with the zinc ions at the antipodes (as part of the complex in which a bromidebromide pair is separated by a zinc ion). As the density is lowered to 3.707 g/cm^{3 }(0.27 cm^{3}/g), the first peak begins to split in two, (in contrast to the melt PDF, which never splits). The split persists for all lower densities out to 1.914 g/cm^{3} (0.52 cm^{3}/g). The distance separating the first and second peaks of the split is about 1 Å. Apparently the zinc ions at lowtointermediate densities penetrate two kinds of empty spaces between the bromides, where the mean Zn^{2+}Br^{}Zn^{2+} angle is at about 110º, as well as the highdensity case of about 180º.
The first peak of the Br^{}Br^{} PDF (Figure 5) at high densities is at the expected distance of about the diameter of a bromide in a crystal. As the density is lowered to 3.291g/cm^{3 }(0.30 cm^{3}/g), the first peak also begins to split in two, (again in contrast to the melt PDF, which never splits). The pressure of the inherent structure is still positive where the split occurs. The split persists for all lower densities out to 1.914 g/cm^{3} (0.52 cm^{3}/g).
Figure 5. The Br^{}Br^{} inherent structure PDF compared with the melt.
The specific volume is 0.27 cm^{3}/g; axes as in Figure 4. The inset is the region around the first peak of the inherent structure for various specific volumes.
The Zn^{2+}Br^{} PDF (Figure 8) is unremarkable at all densities, for both the melt and the inherent structure. The first peak has its maximum about where one would expect based on the empirical ion diameters.
Figure 6. The Zn^{2+}Br^{} inherent structure PDF compared with the melt
The specific volume is 0.27 cm^{3}/g; axes are as in Figure 4.
IV. Discussion and Conclusion
We have presented the first calculation that we know of for the inherent structure of any ionic melt. We chose a system with both a size and charge mismatch. In some respects, the inherent structure appears little different from the inherent structure of liquids of spherical atoms with only shortranged forces; the Br^{}Br^{} PDF follows this sort of behavior for high densities. However, the absence of new structure in the Zn^{2+}Zn^{2+} PDF (compared with the melt) at high densities reminds us of the results obtained for all of the pair distribution functions for strongly associated liquids. Although the explanation lies with the strong attraction of the zinc to the bromide ion, coupled with the even stronger repulsion between the zinc ions, at the same time that sort of behavior relies on a delicate balance of forces whose outcome would be difficult to guess in advance of the calculations.
For lower densities (starting at positive pressure but persisting for negative pressures), the Zn^{2+}Zn^{2+} PDF does uncover new structure, i.e., a split in the first peak, which is not evident in the melt PDF. This split in the first peak is not seen in the inherent structure of other substances examined so far (see Section I.); for van der Waals substances, for example, it is the second peak that splits. Furthermore, this split occurs at about the same density as the underlying solidsolid phase transition for crystalline zinc bromide. It would be tempting to conclude that the local structure of the ultimately amorphous inherent structure imitates this transition, especially since there is a precedent for such behavior in a model nonionic system. There are at least two obstacles to drawing this conclusion. First, we have not been able to identify the structure corresponding to the split in the Zn^{2+}Zn^{2+} PDF to any part of a known structure for zinc bromide. Second, we have not determined if the split of the first peak of the Zn^{2+}Zn^{2+} and the subsequent split of the first peak of the Br^{}Br^{} inherent structure PDFs at even lower density are both the result of one single structural change or if they are each the result of two distinct structural changes. Therefore the question of the impact (if any) of the solidsolid transition of the crystalline phases upon the inherent structure remains open for this system. Now the zinc chloride has similar phase behavior, and because more realistic models are apparently available for that system, any further investigation of this question might be better carried out for models of the chloride.
As in other studies, we obtained a pressurevolume isotherm (at zero temperature) with stable, metastable, and unstable branches, respectively, with decreasing density. It is interesting to see a van der Waalslike isotherm that is obtained without resort to a model pressurevolume isotherm: The separation into metastable and unstable branches is due not to some contrived mechanical constraint (e.g., a small or confined system) or a defect in a model equation of state, as is the case with the original van der Waals model, but instead is due only to supression of fluctuations by the quenching process. That the transition between the convex and concave part of the equationofstate actually corresponds to the traditional concept of a spinodal is verified by examination of the static structure factor, which begins to diverge for small (but finite) wavenumber as the density is lowered from the metastable to the unstable branch of the pressurevolume isotherm. Our identification of the spinodal is also supported by the rapid growth of volume available to the bromide, suggestive of an arrested nucleation of the vapor in the condensed phase.
We conclude with some suggestions for further work on this subject, in addition to the exploration of the consequences of underlying crystalline transitions on the inherent structure. For example, we regret the lack of experimentally determined critical parameters of zinc bromide (and for that matter, the whole equation of state), because we might have otherwise been able to test the recent proposal that the ratio of the magnitude of the minimum pressure in the inherent structure is 2030 times the critical pressure^{15b}. Our estimate of the minimum pressure at –0.4 GPa would give, by this proposal, an estimate of the critical pressure of 1320 MPa, below that of water (P_{c} = 22 MPa ). For another example, it is not known how the theories that have been employed to describe equilibrium ionic melts^{,} will perform on the inherent structures. In particular, it is not known if the corresponding states theory for equilibrium molten salts will hold for their inherent structures. For yet another example, it is not known if or how the phase separation predicted for the highly dissymmetric binary ionic mixtures (e.g., H^{+}He^{++}, encountered in astrophysics and fusion studies) would appear in the inherent structure. Finally, there remains the elucidation of the StillingerLovett sum rules for the 0^{th} and 2^{nd} moments of the PDFs of the inherent structure. These rules are among the few exact results for equilibrium ionic systems. Although the first of these rules is expected to hold for all densities in view of its derivation from the chargeneutrality condition, we do not know if the second sum rule holds in the inherent structure. We have not been able to reliably calculate either of the moments because the volume term (4pr^{2}) in the integrand conspires with the peaksharpening effects of the quench itself to amplify the oscillations in the pair distribution functions, so much so that these do not decay in the 1458atom system that we employed. In the absence of a reliable estimate of the asymptotic oscillations, between ten and 100 times as many atoms would be needed to calculate the moments of the inherent structure PDFs; such a task is within the capabilities of the most powerful supercomputers, but beyond our reach for the moment.
Acknowledgements
Prof. C. Austen Angell (Arizona State) provided helpful discussions concerning the experimental literature of ZnBr_{2}. Brittany D. M. Hodges (INEEL) assisted with the preparation of the Figures.
Table and Figure Captions
Table 1: Ion pair interaction potentials
The potential energy U is given in kcal/mol and the pair distances r in Å.
Figure 1. Pressurevolume isotherm of, and available volume to the bromide, in the inherent structure of molten ZnBr_{2}
The specific volume is given in cm^{3}/g, the pressure is in GPa. Closed circles mark the pressure of the inherent structure. The solid curve marks the % volume available to the bromide in the inherent structure. Block arrows point to the boundaries of the stable and unstable regions, respectively, of the inherent structure. The doublebar doublearrow marks the metastable region; the smaller solid doublearrow marks the region that contains the spinodal (its width is our uncertainty). Open circles mark the pressure of the calculated melt. The dashed line marks zero pressure or per cent for reference.
Figure 2. The total structure factor for the inherent structure of liquid ZnBr_{2}
The wavenumber k is given in Å^{1}. The specific volumes for each structure factor are listed (in cm^{3}/g) in the legend. The curves are splines through the data.
Figure 3. Pair distribution functions for the melt
The top, middle and bottom panels correspond to the zinczinc, bromidebromide, and zincbromide PDFs, respectively. The pair distance r is given in Angstroms. The legend labels the curves according to density (g/cm^{3}).
Figure 4. Zn^{2+}Zn^{2+} inherent structure PDF compared with the melt.
The specific volume is 0.27 cm^{3}/g (or 3.707 g/cm^{3}). The distance r is in Angstroms. The first peak is beginning to split. The inset is the region around the first peak of the inherent structure for various specific volumes.
Figure 5. The Br^{}Br^{} inherent structure PDF compared with the melt.
The specific volume is 0.27 cm^{3}/g; axes as in Figure 4. The inset is the region around the first peak of the inherent structure for various specific volumes.
Figure 6. The Zn^{2+}Br^{} inherent structure PDF compared with the melt
The specific volume is 0.27 cm^{3}/g; axes are as in Figure 4.
References
1
2 For an exception see F. H. Stillinger and R. A. LaViolette, Phys. Rev. B 34 5136 (1986).
3 J.P. Hansen and I. R. MacDonald, Theory of Simple Liquids, 2^{nd} ed. (Academic Press, Orlando FL, 1986).
4 A. Pohorille, L. R. Pratt, R. A. LaViolette, M. A. Wilson, and R. D. MacElroy, J. Chem. Phys. 87 6070 (1987).
5
6
7
8 F. H. Stillinger, T. A. Weber and R. A. LaViolette, J. Chem. Phys. 85 6460 (1986).
9 R. A. LaViolette and F. H. Stillinger, J. Chem. Phys. 82 3335 (1985).
10 R. H. Doermus, Glass Science, 2d ed., (Wiley, New York, 1994).
11
12
13 M. Wilson and P. A. Madden, J. Phys.: Condens. Matter 11 A237 (1999).
14 W. B. Doe, Argonne National Laboratory Report ANL4879 (1952).
15
16 M. J. L. Sangster and M. Dixon, Adv. Phys. 25 211 (1976).
17 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987).
18
19 CRC Handbook of Chemistry and Physics, 73^{rd} ed., edited by D. R. Lide (CRC Press, Boca Raton FL, 1992).
20 N. Karasawa and W. A. Goddard III, J. Phys. Chem. 93 7320 (1989).
21R. A. LaViolette and F. H. Stillinger, J. Chem. Phys. 85 6027 (1986).
22 W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C (Cambridge University Press, New York, 1992).
23 An alternative and more detailed treatment of the available volume or void distributions would employ Voronoi polyhedra, e.g.,
24 P. G. Debenedetti, Metastable Liquids (Princeton University Press, Princeton, 1996).
25 H. B. Callen, Thermodynamics and an Introduction to Thermostatics, 2nd ed. (Wiley, New York, 1985).
26 G. J. Janz, F. W. Dampier, G. R. Lakshminarayanan, P. K. Lorenz and R. P. T. Tomkins, Molten Salts: Volume 1, Electrical Conductance, Density, and Viscosity Data (National Bureau of Standards, Washington D. C., 1968).
27 E. M. Levin, C. R. Robbins and H. F. McMurdie, Phase Diagrams for Ceramists: 1969 Supplement Volume 2, (American Ceramic Society, Columbus OH, 1969) Fig. 2949. We graphically extrapolated the transition to zero absolute temperature to find a transition pressure of about 1.8 GPa, which according to our Figure 1 is at about 0.28 cm^{3}/g.
28 R. A. LaViolette and D. M. Stump, Phys. Rev. B 50 5988 (1994).
29 The existence of the the van der Waals "loop" and the spinodal in unconstrained onecomponent systems has been controversial and denied by many, e.g.,
For a variety of alternatives, of which some are also discussed in (e) above, see e.g.,
and Refs. 15 and 24.
30
31
32
33