Ind. Eng. Chem. Res. 2011, 50, 227–233
227
Nonadditive Hard-Sphere Reference Model for Ionic Liquids Leslie V. Woodcock* Colburn Laboratory, Center for Molecular and Engineering Thermodynamics, UniVersity of Delaware, 150 Academy Street, Newark, Delaware 19716, United States
A nonadditive hard-sphere (NAHS) reference model for ionic liquids with truncated ion-ion pair interaction potentials represented by collision diameters defined as σij ) σ(1 + δijε) (using the Kronecker δ)
is investigated. In a reference 1:1 ionic salt (e.g., AB ) NaCl), the nonadditivity parameter (ε) is positive; hence, like ions interact with a larger diameter than unlike ions; that is, σAA ) σBB ) (1 + ε)σAB
The NAHS model is found to show remarkable similarity of structure to the full-Hamiltonian for electrostatic charged ion models of ionic liquids, their crystalline phases, and freezing properties. Molecular dynamics computations are used to show that the structures of ionic liquids and crystals are determined mainly by short-range geometric and symmetry effects rather than the long-range electrostatic potentials. A preliminary phase diagram for nonadditive hard spheres has been determined for the range of positive nonadditivity for stable cubic lattices CsCl (body-centered), rock salt, NaCl (simple cubic), and sphelerite ZnS (tetrahedral). Model NAHS ionic liquids obey a scaling law for the freezing temperature Tf ) T0(1 + ε)3, where T0 is the HS freezing temperature. NAHS is an effective reference model, analogous to the HS model for molecular and atomic liquids, as a starting point for describing the structures and phase diagrams of ionic liquids and crystals. 1. Introduction The hard-sphere pair potential defines a model Hamiltonian for atomic or molecular fluids that have played a central role in the statistical theory of liquids. It has only one parameter to model a steep positive repulsion when atomic electron clouds overlap at short range. The potential energy of interaction between two hard spheres of diameter σ is φij ) 0 φij ) ∞
for rij > σij
for rij < σij
Unlike simple liquids, ionic liquids, although thermodynamically a single component, are essentially two-component in the statistical mechanical sense that three pair potentials are required to specify the Hamiltonian for a two-body ion-ion interaction model. It would at first appear that this complexity, together with the asymmetry of size, and the long-range Coulombic potentials, that anything resembling a hard-sphere reference model could not reproduce the essential structure. Surprisingly, ionic liquids may be more amenable to a perturbation treatment based upon a hard-sphere reference model than even molecular liquids. One reason is that the Coulombic potential energy per ion is very much greater than the characteristic thermal energy kBT, on the order of 50 times kBT at freezing. Consequently, like ions do not get close enough to interact or “collide” at contact range. The scaling consequence of this is that the thermodynamic properties of simple ionic liquids, such as molten alkali halides, scale with a single distance parameter (r0) that is the minimum in the unlike ion-ion pair potential, as also do thermodynamic * To whom correspondence should be addressed. E-mail:
[email protected].
properties of molecular organic ionic liquids.2 Another observation is that the “effective” pair potentials for ionic liquids are essentially short-range3 due to an effect of local electroneutrality, or Coulomb cancellation. Tosi and co-workers4 successfully calculated partial pair distribution functions and thermodynamic properties of molten alkali chlorides starting with symmetric ionic pair potentials. They found that the interplay of Coulomb attractions and repulsions in the molten salt require only the inclusion of excluded-volume effects between the various ion pairs. The hard-sphere (HS) model pair potential above is easily generalized by a single dimensionless parameter (ε) to the nonadditive hard-sphere (NAHS) model by using the Kronecker δ. σij ) σ(1 + δijε)
(1)
and the collision diameters for a binary A-B system of symmetric 1:1 spheres are σAA ) σBB ) σAB(1 + ε)
(2)
The nonadditivity ε varies from -1 to +infinity; the onecomponent pure hard-sphere fluid corresponds to ε ) 0. Evidence that the nonadditive reference model may represent the structure-generating function part of the full Hamiltonian for ionic liquids can be gleaned from inspection of the widely used Coulombic pair potentials for alkali halides (Figure 1). In the spirit of the WCA perturbation theory for simple liquids,5 the unlike ion-ion Coulomb potential can be truncated at the minimum in the pair potential. The like ion-ion pair potential can be truncated at the point when the repulsive energy at the distance balances the attraction of the unlike ion pair potential; that is, to obtain a WCA-like n-m ionic pair potential (Figure 1).
10.1021/ie101601v 2011 American Chemical Society Published on Web 10/15/2010
228
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
Figure 1. Top: n-m (8-1) ionic liquid ion-ion pair potentials (LHS) and corresponding pair forces (RHS) for NaCl. Bottom: truncated pair potentials in the spirit of Weeks-Chandler-Anderson (WCA; ref 5) model Hamiltonian (LHS), and corresponding truncated ion-ion pair forces (RHS) showing how the NAHS (nonadditive hard-sphere) model is obtained and how the nonadditivity parameter (ε ∏) is determined. Pair energies and forces are reduced in units of minimum cation-anion dissociation energy at distance r0. Red, Na+-Cl-; green, Cl--Cl-; blue, Na+-Na+.
eter is dimensionless, the properties of NA hard spheres obey just the same scaling laws as the hard-sphere model itself. Temperature and pressure cannot be varied independently, and the single reduced thermodynamic state function of reduced number density (F* ) NσAB3/V) and nonadditivity parameter (ε), either 3 P* ) PσAB /kBT(F*, ε)
(3)
T* ) 1/P*
(4)
or reciprocally
Figure 2. Radial distribution function for a WCA reference model of a 1:1 ionic liquid. In this figure, T* is kBTe2/r0 and V* is V/Nr03. r0 is the distance of unlike ion minimum pair potential, T and V being the temperature and molar volume approximately of molten NaCl at freezing.
When this reference Hamiltonian is used to calculate the structure of, say, liquid sodium chloride at its melting point, the radial distribution functions thus obtained (Figure 2) are remarkably similar to those computed from the full-Hamiltonian simulation of model long-range ion-ion pair potentials for sodium chloride, which, in turn, are in good agreement with experiment.4 An inspection of the forces obtained for the truncation of the WCA model in Figure 1 suggests that the properties of the nonadditive hard-sphere (NAHS) system might serve as a reference model for ionic liquids and ionic crystals generally; that is, playing the same role that the additive hard-sphere model plays in the general theory of simple liquids1 and in the WCA perturbation theory in particular.6 The NAHS model then has no other distance scale except the hard-sphere diameter, and because the nonadditivity param-
incorporates all there is to know about the thermodynamic (p-V-T) equation of state. This simple model is consistent with what we already know about corresponding states behavior of the molten alkali halides; that is, that all conformal ionic liquids or crystals with the same σAB, which can be identified with r0 in the pair potential (Figure 1), will obey the same law of corresponding states and have the same freezing point, for example.7 To determine to what degree the nonadditivity ε determines the crystal structure, we need to investigate how different crystal structures are determined by the nonadditivity. 2. Crystalline State Phase Boundaries Before proceeding to determine the properties of the NAHS model at or above freezing, we can immediately specify the most stable crystal structures at low temperatures, that is, when the temperature is zero or the pressure is infinite. In this limit, the crystal structure of maximum close packed density is thermodynamically the most stable. The density of maximum packing of binary nonadditive sphere crystal structures can be calculated analytically as a function of the nonadditivity parameter, ε. Results are given in Table 1 for the fcc structure
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011 Table 1. Maximum Densities (G0* ) NσAB /V0, where Vo is the Minimum close-packed volume per sphere) of Crystal Phases for Non-Additive Spheres for the Various Cubic Structures: Face-Centered Cubic (fcc), CsCl (body-centered cubic), NaCl (simple cubic), and ZnS (sphelerite or tetrahedral) 3
fcc CsCl NaCl ZnS
ε < 0 (2 phases) ε 2/3 ε < 2 ε > 2 ε < 22/3 ε > 22/3
Fo* ) 2 Fo*(ε) ) 2/(1 + ε)3 Fo* ) (33)/4 Fo*(ε) ) 2/(1 + ε)3 Fo* ) 1 Fo*(ε) ) 22/(1 + ε)3 Fo* ) 33/8 Fo*(ε) ) 22/(1 + ε)3
of the hard-sphere crystal and the three crystal structures with cubic symmetry found in nature for ionic crystals, respectively. From these close packed limits, the most stable phases and the value of the nonadditivity parameter at which the phase transitions between these alternative crystal structures occur in the limit of low temperature (T f 0 K) can be calculated. At the phase boundary, the maximum density expression for any two phases is set equal, and the resulting equation can be solved for ε. The maximum densities of the crystal structures are shown in Figure 3. It can be deduced from a principle of maximum packing that at all negative values of the nonadditivity (-1 < ε < 0), the most stable state comprises the two separated fcc phase A and fcc phase B. For low negative nonadditivity, there will be a region of solid solution, that is, close to zero nonadditivity, whereupon the two components are indistinguishable and become one component. The fcc structure is the most stable crystal structure of one-component, pure hard spheres. As the nonadditivity is increased from zero, in the limit of low temperature at constant pressure, or equivalently, limit of high pressure at constant temperature, inspection of Figure 2 shows that as the nonadditivity increases positively from 0, the fcc with coordination number 12 becomes unstable with respect to CsCl, with coordination number 8, at the point where the lines of maximum density cross. At higher values, the CsCl becomes unstable relative to 6-coordinated NaCl, which then competes
229
with the ZnCl at a density at which these two structures become the same. Nonadditivity values at the three transition points are fcc to CsCl
ε(t1) ) 2(5/6) /3(1/2) - 1 ) 0.02875
CsCl to NaCl
ε(t2) ) 21/3 - 1 ) 0.25992
NaCl to ZnS
ε(t3) ) (16/3)1/3 · (2/3)1/6 - 1 ) 0.63299
In Figure 3, the density of maximum close packing of the crystalline phases of the model as a function of nonadditivity using two alternative distance scales. One can define a set of reduced properties either in dimensions of σAB (denoted by superscript *), F* ) NσAB3 /V
(5)
where V is the volume occupied by N spheres, or one can equally define the reduced properties in terms of the distance scale, σAA, (denoted by a superscript #), whereupon 3 F# ) NσAA /V
(6)
Figure 1 further shows that there are two different scaling regimes, for low and negative nonadditivity. The freezing transitions are scaling as σAB, but for higher values of nonadditivity, the maximum density scales as σAA. When ε becomes very large, the two components A and B crystallize independently of the presence of the other. At this point, the maximum packing properties of hard spheres are recovered, but the density of the system has increased by a factor of 2. In this range, the model system has either NaCl or ZnS structures, which eventually show properties identical to the original hard-sphere model. For these systems at very high nonadditivity, there are two hard-sphere systems effectively coexisting in the same volume and oblivious to the presence of the other. We will see below, that this interesting limiting structure also manifests itself in the phase transition properties of the model for large ε at finite temperatures or pressures. 3. Thermodynamic Equations of State: Molecular Dynamics Simulations To determine thermodynamic equation-of-state data for model NAHS systems, we use a standard hard-sphere collision dynamics program using the Alder-Wainwright method with the model generalized to incorporate a binary system of 8000 spheres with the nonadditive hard-sphere potentials. The program calculates the equation-of-state pressure for fluid phases or for any of the four cubic crystal structures considered in section 2. The pressure is calculated both from the virial theorem and the collision rate expression. When the system is in equilibrium and the averaging is on the order millions of collisions, these two methods give the same pressures. The total pressure, P, is resolved into partial pressures from the three different types of collision, P ) PAA + PAB + PBB
(7)
for a symmetric 1:1 salt with size ratio σAA/σBB ) 1, then PAA ) PBB and Figure 3. The infinite-pressure, or zero-temperature, maximum packing density of a system of binary symmetric nonadditive hard spheres as a function of the nonadditivity parameter (ε) when the reduced density is in units of σAB (above) and when the reduced density is in units of σAA(below).
P* ) F* + (2(1 + ε)〈τAA〉 + 〈τAB〉) √π/3
(8)
where 〈τΑΑ〉 and 〈τΑΒ〉 are the mean collision rates, per sphere per unit time, of types AA and AB, respectively, in dimension-
230
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
Figure 4. Radial distribution functions for solid and liquid phases of the NAHS model at reduced densities corresponding to sodium chloride at its melting point. Red, Na+-Cl- unlike ion; blue, Cl--Cl- and Na+-Na+ like ion pair distributions. The respective configuration snapshots show how the NAHS model, with no dissimilarity of A and B, and no long-range effects, reproduces “ionic” liquid and “ionic” crystal structures.
less HS units of frequency (kBT/mσAB2)1/2, where m is the mass of a sphere, and kB is Boltzmann’s constant. To compare with experiment and to find the extent these simple scaling laws relate to real experimental systems, we consider to begin with a model of sodium chloride with a nonadditivity of 1.5, as shown in Figure 1. The radial distribution functions for NAHS-NaCl at the experimental melting and freezing densites of NaCl, of crystal and fluid phases, respectively, are shown in Figure 4, together with sample equilibrium configurations. The crystal RDFs are essentially the same as those obtained from Monte Carlo calculation for full ion-ion Hamiltonian models. Likewise, the liquid NAHS model RDF at the melting point also shows all the characteristic features beyond r0 of liquid NaCl at the melting point, with the same peak height at, except for the sharp cusps, contact distances σAB () r0) and σAA () 1.5r0). This is a consequence of the hardsphere discontinuity. Presumably, this artifict could be alleviated by bringing the reference model into line with real experiment simply by replacing the nonadditive hard-sphere with nonadditive slightly soft spheres. The equation-of-state pressure data obtained for the 1:1 NAHS fluid for ε ) 0.2, that is, within the CsCl stability range in Figure 2 and for the CsCl crystal phase, obtained as a function of density are shown in Figure 5. Just as with the one-component hard-sphere fluid, the fluid branch extends beyond the freezing transition into the metastable supercooled region, and the crystal branch can be extended until the crystal becomes mechanically unstable at the reduced density around 0.7. Accurate prediction of the coexistence properties requires a detailed computation of the chemical potential along a reversible path connecting the two phases. For the present purposes, to obtain a preliminary insight into the phase diagram, however, it suffices to use approximate statistical theories of melting to estimate the phase transition.8-10 In the case of the two-component system, however, one can see that the main difference between the entropy of mixing is
Figure 5. Equation-of-state data for nonadditive spheres with a low nonadditivity in the CsCl stable crystal-structure region; Liquid phase, red; crystal phase, blue. The horizontal line is estimated coexistence pressure.
NkB loge 2 where loge 2 is minus the ideal entropy of mixing for an equimolar binary ideal gas mixture () x1 ln x1 + x2 ln x2). This arises because the entropy of mixing in the fluid is loge 2, whereas there is no entropy of mixing for the crystal, because all lattice sites are allocated. For a one-molar mixture, the crystal and liquid then have equal Gibbs chemical potential when the pressure is approximately NkB loge 2 /∆V where kB is Boltzmann’s constant. This simple first approximation was suggested by Frank van Swol;8 as we will see below, it is quite an accurate preliminary estimate of the liquid-crystal coexistence pressure of these binary systems. Another way to obtain a reasonably accurate estimate of the coexistence densities is to use the Alder-Ross rule,9 that is, to
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
231
Figure 7. Equation-of-state data for nonadditive spheres with a low nonadditivity in the ZnS region showing liquid phase (red circles) and of the alternative crystal phases, NaCl phase (blue circles) and ZnS phase (green circles).
Figure 6. Equation-of-state data for nonadditive spheres with nonadditivity in the NaCl crystal-structure region. Red circles, liquid; blue circles, NaCl crystal; green circles correspond to ZnS structure. The horizontal line is the estimated coexistence pressure.
find the density at which the crystal becomes mechanically unstable, and to equate that with the liquidus density at coexistence. This seems to be sufficiently accurate for all these hard-core systems, and indeed, it is the Lindemann law of melting, which becomes an exact scaling law for the general inverse power hard or soft sphere model.10 Interestingly, the root-2 entropy prediction and the Alder-Ross rule give approximately the same coexistence line for the ε ) 0.2, CsCl crystal. Next, we examine in more detail a higher value of ε in the region of stability of the NaCl simple cubic rock salt structure. The equation-of-state data and the estimated transition are still within the low temperature stability region of NaCl, but at finite temperatures, we see from Figure 2 that the ZnS 4-coordination tetrahedral structure is also possible for the crystal phase. At this nonadditivity, the data in Figure 6 shows that the sodium chloride 6-coordinated structure has a slightly lower pressure than the tetrahedral (ZnS) alternative. This means that at coexistence, at this ε, the NaCl structure is the more stable for the coexisting crystal in addition to being the more stable at 0 K, as seen in Figure 2. When ε is increased to 0.75, however, into the region for which the maximum densities of 6-coordinated NaCl and 4-coordinated ZnS crystal structures become equal, the pressures are very close together, suggesting that the more stable is a close call. Figure 6 suggests the ZnS structure may become marginally the less stable because it has a slightly higher pressure in the vicinity of melting. MD data for the respective crystal structures in the two-phase region show the ZnS pressure now to be slightly lower than the NaCl case. It seems likely that the freeenergy difference between these two phases for higher nonadditivity (i.e., ε > 0.633) may result in a more complex phase diagram, with regions of stability of both phases being possible. Increasing the value of ε further, it can be predicted analytically that the end result of very large ε, that is, in the limit σAA/σAB f ∞, we recover a system that has exactly twice the number density of the hard-sphere fluid, but that has exactly the same value of reduced pressure or temperature as the pure one-component HS model for its intensive thermodynamic state
Figure 8. Equation-of-state data for nonadditive spheres with a high nonadditivity in the ZnS/NaCl crystal-structure region; the solid black fluid, crystal, and dashed coexistence lines are one-component HS.
functions. What seems quite surprising is that this limit is closely approached at a value of ε of only 2, as can be seen from Figure 8. In Figure 8, the equation-of-state data for the fluid and ZnS crystal phases of the NA sphere model is superimposed upon the corresponding equation-of state pressures and coexistence line of pure hard spheres. The results show that the NAHS model properties are already beginning to scale like one-component additive hard spheres at a surprisingly low value of ε ) 2. Here, the two components A and B, are almost oblivious to the presence of each other, and the equation of state of two superimposed hard-sphere systems is recovered. Figure 8 shows the equations of state are almost predictable from the onecomponent HS fluid. At this stage, the reduced freezing pressure/ temperature of the two-component NAHS is the same as onecomponent HS. 4. Phase Coexistence and Scaling In Table 2, we summarize the coexistence data from these preliminary simulation studies so far undertaken. The data for pure HS is taken from Chang and Sandler.11 Also given in that table are the coexistence parameters for the bcc fluid coexistence of pure spheres.12 The latter is obtained from the Ree-Hoover single occupancy technique that stabilizes the bcc structure down
232
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
Table 2. Coexistence Properties of 1:1 Nonadditive Spheres ε
lattice
P*
T*
T#
F* (liq)
F* (c)
F# (liq)
F#(c)
% ∆Vf
0.00 0.00 0.0287 0.10 0.10 0.20 0.25 0.25 0.2599 0.30 0.50 0.6329 0.75 0.75 1.00 2.00
fcc bcc fcc f CsCl bcc CsCl CsCl CsCl NaCl CsCl f NaCl NaCl NaCl NaCl f ZnS NaCl ZnS ZnS ZnS
11.92 14.5 infinity 15.05 8.77 7.79 7.36 8.13 infinity 6.105 4.735 infinity 2.7568 2.8433 2.3961 0.5581
0.08289 0.06897 0 0.0665 0.1145 0.1283 0.1359 0.1229 0 0.1638 0.2112 0 0.3627 0.3517 0.4173 1.7918
0.08292 0.06897 0 0.0499 0.0857 0.0743 0.0696 0.0629 0 0.0745 0.0741 0 0.0677 0.0656 0.0522 0.0664
0.9458 0.955
1.0463 0.990
0.943 0.955
1.040 0.990
10.29 3.66
0.9966
1.0804
1.326
1.438
8.41
0.3782 0.360 0.365
0.4386 0.400 0.415
0.654 0.703 0.713
0.758 0.781 0.811
15.97 11.11 13.70
0.331 0.477
0.370 0.573
0.725 1.610
0.813 1.934
12.12 20.13
0.300 0.310 0.225 0.065
0.350 0.350 0.252 0.073
1.608 1.661 1.800 1.755
1.876 1.876 2.016 1.958
16.67 12.90 12.00 11.54
to low density; for unconstrained hard spheres, the bcc crystal is mechanically unstable with respect to slip along the 1:1:1 Miller diagonal slip plane. Using this knowledge of hard-sphere phase coexistence parameters for fcc and also bcc crystal structures, the various equation-of-state data discussed in section 3, and the zero-K maximum densities, we can now begin to construct something resembling a phase diagram. Figure 9 is a preliminary sketch of the NAHS model phase boundaries for all values of ε from -1 to +1; there seems to be nothing much more to discover between +1 and infinity. Although this analysis for ionic liquids is concerned only with positive nonadditivity, we note that for values of ε from -1 to 0, the stable crystal is two-phase; that is, separated fcc crystal lattices for the two components A and B. This is simply because the maximum packing density, or minimum pressure, is obtained if the system phase separates. At this stage, we cannot say to what extent, if any, phase separation takes place in the liquid region as the nonadditivity approaches -1. It is interesting to note, however, that there is also an experimental counterpart for the negative nonadditivity model in systems of binary metal alloys.13 Metal atoms always see a different metal’s atoms as being larger than they see themselves; that is just the opposite
of ionic liquid ion-ion pair potentials. Binary metallic mixtures readily supercool down to below a glass transition temperature to form homogeneous amorphous solids. The reason for this general phenomenon can be seen in the NAHS model. Since the thermodynamically stable phase is two separate coexisting fcc crystals of A and B, the crystal nucleation requires phase separation, assuming complete mixing in the liquid phase. The phase separation time scale becomes too long for the rate at which the systems are cooled, and the glass transition, rather than crystallization, occurs. Then this two-phase crystal system, at a very small positive value of the nonadditivity, becomes unstable with respect to the ionic salt CsCl body-centered cubic structure with 8-coordination. The data for the freezing transition in this region appears to follow a scaling law that T* ) T*(1 + ε)3 0
(9)
where T0 is the freezing transition for the bcc hard-sphere system. For negative nonadditivity, T0 is the fcc hard-sphere fluid phase transition, and the scaling law appears to be exact when separation of A and B into two identical phases occurs to counteract the effect of negative nonaddittivity. The observations about scaling from Figure 9 lead to the interesting conclusion that all ionic liquids should have the same reduced freezing transition temperature if the reduced temperature is expressed in units of the characteristic length, σAA. For ionic liquids at positive nonadditivity, when the reduced melting temperature is given by T# ) kBT/PσAA3
(10)
it appears to be almost constantly independent of ε, starting with the HS freezing point at ε ) 0 and approaching the same constant value at ε ) ∞. Whence, for all positive ε, T#(NAHS) ∼ T#(HS)
Figure 9. Phase diagram for NAHS model with reduced temperature in dimensions of σAB: the points are the estimates from MD computer simulations as discussed in the text. for -ve NA spheres, the solid line is the Tf(1 + ε)3 scaling law that passes through the hard-sphere fcc melting temperature up to ε ) 0; for +ve NA spheres, the scaling constant appears to be the HS fluid-bcc coexistence.
(11)
Ionic liquids, the molten alkali halides in particular, differ from all other liquidsssimple liquids; liquid metals; network liquids, such as water or silicasby having a large volume decrease on freezing, around 20%. The data in Table 2 and Figure 10 above suggest that the freezing parameters of ionic liquids, for example, which freeze to the same simple cubic (NaCl) ionic crystals, to a first approximation are independent of temperature or pressure. 5. Conclusions The forgoing analysis and thermodynamic properties from MD computations show that the general NAHS model can play
Ind. Eng. Chem. Res., Vol. 50, No. 1, 2011
233
structures for 1:1 salts without invoking long-range Coulombic forces or ionic Madelung constants. Finally, we note the industrial importance of various classes of organic ionic liquids, that there are many such ionic liquids whose structural and freezing properties are central to the functionality, and that these liquids can be designed according to their predicted properties. Design criteria can be obtained from correlated knowledge of existing ionic liquids14 or from ab initio quantum simulations.15 Predicting freezing points on the basis of statistical mechanics using simple models such as NAHS can compliment these “search and design” techniques at the outset. For predicting freezing transitions, perturbation treatments, starting with the scaling properties of NAHS models, may be just as important an alternative approach as that of the QSPR (quantitative structure-property relationship) method, which correlates the properties of ionic liquids that are presently known16 to predict freezing points of yet-to-be-synthesized ionic liquids. Figure 10. Equation-of-state and phase behavior of a NAHS model sodium chloride.
a role in describing and understanding the freezing properties of ionic liquids. It is often said that when computer simulation is used as a research tool, it is not what we put into simulation that informs us about the science of real complex systems, but what we leave out. This approach to understanding and describing structures of complex molecular systems can be termed “Hamiltonian surgery”. Using the full Hamiltonian of a pairwise Coulomb model, we already have a corresponding states explanation of the freezing points of alkali halides and the existence of lowtemperature organic ionic liquids with large ions.2,7 For any given value of nonaddittivity, the simple corresponding states scaling of all the alkali halides is easily predicted. With the NAHS model, we now have an explanation for why the cesium halides, for example, have a structure different from the other alkali halides. The NAHS model also offers an explanation as to why the lithium halides, for example, fall within a different conformal group of their own, compare from the halides of K and Na, for example. It has been known for many years that the Li halides follow a scaling law based only on the characteristic size of the anion (not r0). This type of scaling is found here with the NAHS model with the freezing points scaling as σAA3 (eq 5). The nondimensional nonadditivity parameter can be obtained from conventional spherically symmetric model Coulombic pair potentials of ionic salts. The NAHS model, thus defined, can be used to explain many of the general freezing properties of ionic liquids, molten alkali halides in particular. When the nonadditivity is low, for example, with cesium chloride, the crystal structure is CsCl-body-centered. As the nonaddittivity increases, the NaCl structure prevails; when the nonadditivity is high, as with ZnS, the tetrahedral sphelerite structure becomes the more stable, in accord with real experiment on ZnS type compounds with high nonadditivity. The NAHS model offers an alternative understanding of the general trend of ionic crystal
Literature Cited (1) Hansen, J.-P.; Mc Donald, I. R. Theory of Simple Liquids, 4th ed.; Academic Press: London, 2008. (2) Barbaro, R.; Jiang, J.-W.; Woodcock, L. V. Corresponding states theory for the freezing of ionic liquids. Ind. Eng. Chem. Res. doi: 10.1021/ ie101653n. (3) Clarke, J. H. R.; Smith, W.; Woodcock, L. V. Short-range effective pair potentials for ionic liquids. J. Chem. Phys. l986, 84, 2290. (4) Ballone, P.; Pasture, G.; Tosi, M. P. Structure and thermodynamic properties of molten alkali chlorides. J. Chem. Phys. 1984, 81, 3174. (5) Weeks, J. D.; Chandler, D.; Anderson, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237. (6) Ben-Amotz, D.; Stell, G. Reformulation of Weeks-ChandlerAnderson Perturbation Theory directly in terms of hard-sphere reference system. J. Phys. Chem. 2004, B108, 6877. (7) Woodcock, L. V. Corresponding states theory for the fusion of ionic crystals. Proc. R. Soc. (London) A 1976, 348, 187. (8) van Swol, F. Private communication: Sandia National Laboratory, 2002). (9) Ross, M.; Alder, B. Melting curve at high pressure. Phys. ReV. Lett. 1966, 16, 1077. (10) Hoover, W. G.; Ross, M. Statistical theories of melting. Contemp. Phys. 1971, 12, 339. (11) Chang, J.; Sandler, S. I. Determination of liquid-solid transition using histogram reweighting method and expanded ensemble simulations. J. Chem. Phys. 2003, 118, 8390. (12) Woodcock, L. V. Computation of the free energy for alternative crystal structures of hard spheres. J. Chem. Soc., Faraday Discuss. 1997, 106, 325. (13) Gerhard, K. Non-additive hard-sphere reference system for a perturbative liquid state theory of binary systems. J. Chem. Phys. 1990, 93, 5105. (14) Larsen, A. S.; Holbrey, J. D.; Tham, F. S.; Reed, C. A. Designing ionic liquids: imidazolium melts with inert carborane anions. J. Am. Chem. Soc. 2000, 122, 7264. (15) Turner, E. A.; Pye, C. C.; Singer, R. D. Use of ab initio calculations towards rational design of room temperature ionic liquids. J. Phys. Chem., A 2003, 107, 2271. (16) Eike, D. M.; Brennecke, J. F.; Maginn, E. J. Predicting melting points of quarternary ammonium ionic liquids. Green Chem. 2003, 5, 5323.
ReceiVed for reView July 27, 2010 ReVised manuscript receiVed September 21, 2010 Accepted September 29, 2010 IE101601V
