Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A statistical mechanical theory for the crystallization of alkali halides from aqueous solution Ursenbach, Charles 1992

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1992_spring_ursenbach_charles.pdf [ 6.63MB ]
JSON: 831-1.0061751.json
JSON-LD: 831-1.0061751-ld.json
RDF/XML (Pretty): 831-1.0061751-rdf.xml
RDF/JSON: 831-1.0061751-rdf.json
Turtle: 831-1.0061751-turtle.txt
N-Triples: 831-1.0061751-rdf-ntriples.txt
Original Record: 831-1.0061751-source.json
Full Text

Full Text

A STATISTICAL MECHANICAL THEORY FOR THE CRYSTALLIZATION OF ALKALI HALIDES FROM AQUEOUS SOLUTION By Charles Ursenbach B. Sc, Hon. (Chemical Physics) University of Calgary A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF D O C T O R OF PHILOSOPHY  in THE FACULTY OF GRADUATE STUDIES CHEMISTRY  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA  December 1991 © Charles Ursenbach, 1991  In presenting this thesis  in partial fulfilment  of  the requirements for an advanced  degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department  or  by  his  or  her  representatives.  It  is  understood  that  copying  or  publication of this thesis for financial gain shall not be allowed without my written permission.  (Signature)  Department of  Chemistry  The University of British Columbia Vancouver, Canada  Date  DE-6 (2788)  December 21, 1991  Abstract  T h e crystallization transitions and physical stability of model aqueous alkali halide solutions are investigated.  A second-order density functional theory is developed for  molecules of general symmetry, and is applied to the crystallization of anhydrous alkali halides, ice Ih, sodium chloride dihydrate, and lithium iodide trihydrate from the appropriate solutions.  Thermodynamic stability theory is clarified and then applied,  using Kirkwood-Buff methods, to ionic solutions and a simple binary liquid mixture. The liquid correlation functions are obtained for two liquid models. In the first, ions are modelled as charged hard spheres, and the solvent as a polarizable multipolar hard sphere. The second is similar but involves a modification to the hard repulsive core of the ion-ion potential, and this is found to have a strong influence on solution properties and crystallization transitions. The density functional theory presented is the first such theory for nonlinear molecules based on generalized spherical harmonic expansions. It is also the first molecular theory with the density given in direct lattice vectors instead of reciprocal lattice vectors. Because of the long-range correlations in ionic solutions, a technique involving Ewald sums is developed to aid the convergence of sums over direct lattice vectors. The resulting theory is shown to be superior for these models to a reciprocal lattice vector method. In its final form, the theory for molecules such as H20 requires a number of one-dimensional integrals to be performed numerically. When applied, the density functional theory yields minima for nearly all systems, showing the possibility of a phase transition. In many of these systems, liquid-solid coexistence is predicted as well. Model dependence is shown to be important in determining whether the equilibrium is under conditions similar to  ii  those of real systems, and whether the crystal parameters are similar to those of real solids. The theory predicts reasonable values for the oscillation widths about average orientations. In applying stability concepts, it is important to distinguish between the strictest criteria of stability, a set of conditions which are collectively necessary and sufficient for thermodynamic stability, and the weaker criteria, which are necessary but not sufficient conditions. We use this to clarify some misconceptions in the literature, and then give its nontrivial application to electrolytes. The results indicate that salts with only large ions, such as Csl, and those with a small ion, such as Na+ or K+, behave differently near the absolute stability limit. The former act hydrophobically, and appear to undergo demixing from the solvent, while the latter, which bind the solvent more tightly, do not show clear signs of any such demixing, but do appear to become mechanically unstable.  iii  Table of C o n t e n t s  Abstract  ii  List of Tables  vii  List of Figures  viii  Table of S y m b o l s  x  Acknowledgement  xiv  1  Introduction  1  2  P h a s e Transitions  9  2.1  The Density Functional Theory of Freezing  9  2.2  Thermodynamic Stability Theory  20  2.2.1  Strictest Stability Criteria  23  2.2.2  Instability and Integral Equations  28  3  4  G r a n d Canonical T h e o r y of I n h o m o g e n e o u s S y s t e m s  33  3.1  Grand Canonical Density Functional Theory  33  3.2  Application to Pair Structure Determination  39  3.3  Application to Freezing Theory  45  3.4  Application to Stability  48  Calculation of Liquid P r o p e r t i e s  50  iv  5  4.1  Integral Equation Calculations of Molecular Liquids  50  4.2  Models of Aqueous Alkali Halides  55  Calculation of P h a s e C o e x i s t e n c e  66  5.1  Derivation of the Molecular Density Functional Expressions  66  5.1.1  Ideal Contribution  68  5.1.2  Interactive Contribution: k = 0  70  5.1.3  Interactive Contribution: r-space  72  5.1.4  Treatment of Long-Range Potentials  75  5.1.5  Degenerate Terms of the Functional  79  5.1.6  Effect of Proton Disorder  81  5.1.7  Effect of Solvent Polarizability  83  5.1.8  A Parametrization of the Angular Distribution  84  5.1.9  Summary  90  5.2  5.3  5.4  5.1.10 Comparison with the fc-space Expression  91  Anhydrous Alkali Halides  94  5.2.1  Methods and Crystal Structures  94  5.2.2  Results for Charged Hard Sphere Models  97  5.2.3  Sensitivity to Model  103  5.2.4  Direct Space vs. Fourier Space Methods  105  Icelh  Ill  5.3.1  Crystal Structure  Ill  5.3.2  Results and Discussion  115  Alkali Halide Hydrates  125  5.4.1  Crystal Structure of NaCl-2H 2 0  125  5.4.2  Crystal Structure of LiI-3H 2 0  127  v  5.4.3 5.5  6  7  Results and Discussion  130  Practical Implementation of the Density Functional Program  136  5.5.1  Programming Tests  136  5.5.2  Series Convergence  139  5.5.3  Computing Resources  141  Calculation of S t a b i l i t y L i m i t s  144  6.1  Derivation of Stability Criteria for Electrolytes  144  6.1.1  Total Correlation Function Expressions  145  6.1.2  Direct Correlation Function Expressions  149  6.2  Hard Sphere and Csl Solutions  154  6.3  KC1 and NaCl Solutions  175  Conclusions  193  Bibliography  199  vi  List of Tables  4.1  Input parameters for model electrolyte solutions  56  4.2  Average numbers of contact and solvent-separated ion pairs  64  5.3  Parameters describing crystal/solution equilibria for different models of anhydrous alkali halides  100  5.4  Spatial and angular distribution widths for water molecules  123  5.5  Contributions to Proton Disorder  124  5.6  Experimental coordinates for NaCl-2H20  127  5.7  Nearest neighbour surface-to-surface distances in NaCl-2H20  134  5.8  Nearest neighbour pair configurations in ice Ih  137  5.9  Comparison of nearest-neighbour energies with Coulson and Eisenberg. .  138  5.10 Contributions to the electrostatic interaction of hydrates  138  5.11 Convergence of lattice vector expansion  140  5.12 Resource requirements for the molecular density functional program. . . .  142  6.13 Total correlation peaks, liquid configurations, and crystal configurations for aqueous Csl  174  vn  List of F i g u r e s  1.1  Phase transitions of the NaCl-water binary mixture  2  4.2  Schematic for RHNC calculation of simple liquids  51  4.3  Schematic for RHNC calculation of molecular liquids  53  4.4  The effect of short-range repulsions on mean ionic activity coefficients. . .  59  4.5  The pair distribution function g+-(r) for 2.0M NaCl with various repulsion parameters  62  5.6  Representative plots of ( / ? A f ) / V ) m m for CHS models  98  5.7  Variation of ( / ? A f 2 / V ) m m with concentration and A* for CHSR9 models.  5.8  The dependence of ( / ? A f i / V ) m i n upon the number of sets of lattice vectors employed for different methods of calculation  5.9  Evaluation of functional (ftAQ/V)  101  106  (no minimization) for selected solid  state parameters  108  5.10 Crystal Structure of ice Ih  112  5.11 Minimized functional values for proton-ordered ice Ih  116  5.12 Spherically averaged direct correlation functions of water near contact.  .  120  5.13 Crystal structure of LiI-3H 2 0  128  5.14 Functional values for alkali halide hydrates  131  5.15 Arrangement of quadrupoles in LiI-3H 2 0  139  6.16 Stability indicators of hard spheres in water  155  6.17 Stability indicators of Csl in water  157  vm  6.18 Weak stability indicators of hard spheres in water and Csl in water. . . .  160  6.19 Total correlation functions for the hard spheres in water  163  6.20 Ionic total correlation functions for Csl in water  165  6.21 Total correlation functions involving solvent for Csl in water  167  6.22 Composite total correlation functions for Csl in water.1  169  6.23 Composite total correlation functions for Csl in water.II  172  6.24 Stability indicators of KCl in water  176  6.25 Stability indicators of NaCl in water  178  g^f-J  near the instability  180  6.27 Ionic total correlation functions for KCl in water  183  6.28 Total correlation functions involving solvent for KCl in water  185  6.29 Composite total correlation functions for KCl in water.1  187  6.30 Composite total correlation functions for KCl in water.II  189  IX  Table of Important Symbols  Symbol  Description  A*  repulsion strength parameter  Cij  integral of c™%(r)  u(Hl2>0.1>ii.2)  liquid direct correlation function  c  ij\LliL2)  solid direct correlation function  c^Uru)  projection of liquid direct correlation function  di  diameter of species i  dij  (di + dj)/2  E  internal energy  e  E/N  F  Helmholtz free energy  m  intrinsic Helmholtz free energy  c  f  (Ch. 2,6, F/N)(Ch.  3, probability density)  r  nonequilibrium probability density  y/(2m + l)(2n + 1)  fipix.)  spatial distribution function  Gij  integral of h™%(r)  G  Gibbs Free Energy  9  G/N  9ip(&)  angular distribution function  ffij{Ll2iilliil2)  pair correlation function  H h  Enthalpy (Ch. 2, H/N){Ch.  X  3, Planck's constant)  liquid total correlation function  "•ij (Hi 2 5 0 . 1 , 0.2 )  hij(Ei,L2)  solid total correlation function  K^iiiru)  projection of liquid total correlation function  h,i  k  moment of particle of species i  la  constant, 8n2 for nonlinear particles  KE  kinetic energy  k  Boltzmann's Constant  k-n  reciprocal lattice vector  M  Madelung constant  mi  mass of particle of species i  mz  order of principal axis of symmetry  Ni  number of particles of species i  ne  number of independent extensive variables  np  1 or 2, depending on symmetry of particle  ns  number of species  P  pressure  Ujij  component of quadrupole moment tensor  Qi  charge of species i  %  generalized coordinates of particle i  R  general thermodynamic potential  r  R/N  U  spatial variables  XI  R  direct lattice vector  R'mni®  generalized spherical harmonic  Rr  vector connecting a pair of particles  s  Entropy  S  S/N  T  temperature  Tr d  classical trace  U  configurational energy  Uij  pair potential between particles of species i and j .  V  volume  V  V/N  Vi  partial molar volume of species i  "" ip  integral of gip(£Li) times -fiyM(B.i)  Wi(q)  Hi - (j>i(q)  X2  mole fraction of component 2  %i i yi  general intensive variables  Xi, Yi  general extensive variables  Xi  general extensive density  Z{R)  degeneracy of lattice vectors of magnitude R  Xll  OLn  d^abc  Fourier coefficient coefficient of Fourier/gen.sph.harm, expansion  P  1/kT  7±  activity coefficient  7ew  Ewald constant  8  width of angular oscillation of the principle axis  tip  Gaussian width of oscillation of i,p particle  Vij(r.i2,  Q>i,Q.2)  liquid series function  X  coupling parameter  p-i  chemical potential of species i  a  dipole vector  £  width of angular oscillation about the principle axis  s  grand canonical partition function  PiL  liquid density of species i  pi(g)  solid density of species i  Pis  average solid density of species i  PC  density of primitive cells in the crystal  p*  reduced density, pdzs  /#W)  two-particle density of species i, j  pi{l)  single-particle density operator for species i,j  lip  location of i,p particle in unit cell  ^(fli,02,£i2)  rotational invariant  $  external potential of system  ^i  external potential of particle i  n  grand canonical potential, — pV  a  angular variables  Xlll  Acknowledgement  I wish to thank my academic advisor, Dr. G. N. Patey, who has given me both guidance and independence throughout the course of this research. Financial assistance has been provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) and by the Chemistry Department of the University of British Columbia. I am also indebted to Dr. P. G. Kusalik for access to programs and results, and to Dr. D. Q. Wei for the use of a number of results. Most of all I would like to thank my wife, Laurel, who thought we were just coming out here for a couple of years. I am especially grateful to her for the many sacrifices she has made during the final months of thesis preparation.  xiv  Chapter 1  Introduction  Phase transitions provide a multitude of intriguing problems for the physical theorist, and are one of the remaining challenges of a general nature in classical statistical mechanics. Consider the phase transitions possible for a solution of NaCl in water at 1 atm pressure, indicated in Figure 1.1. We can explain in terms of thermodynamics the freezing and vaporization of water, and the freezing point depression and boiling point elevation that result as the salt is added. Thermodynamics also describes the saturation concentrations and eutectic points associated with the salt crystals. At higher temperatures than shown in the figure the anhydrous salt melts, still in coexistence with the salt water vapor, and if other pressures were included we would encounter critical phenomena as well, all accounted for within the laws of thermodynamics. The thermodynamic picture is within itself complete and self-consistent, but there is clearly much more to be learned about these coexisting phases. To explain by the methods of statistical mechanics the existence of phase transitions, we can formally begin at the partition function, from which all equilibrium behaviour should be predicted.  In the early part of the century, there was some debate as to  whether the partition function, an integral of continuous functions, could predict the inherent discontinuities of a phase transition [1]. It is clear now that phase transitions are macroscopic events, defined only in the thermodynamic limit. The limit of a function may be discontinuous, even if the function itself is not, a simple example being the delta function defined as the limit of a Gaussian distribution. The famous work of Yang and Lee  1  Chapter 1.  120  Introduction  r(°C)  100 806040-  20 0 -20-40  2  3 4 molality  Figure 1.1: NaCl in water at 1 atm pressure. The enclosed region gives the range of temperature and concentration where the liquid state is stable. The labels along each segment of the curve describe the coexisting phase. Data are taken from the International Critical Tables.  Chapter 1.  Introduction  3  [2] proved that it is possible for the partition function to have a discontinuous equation of state in the macroscopic limit, by showing that zeros of the partition function, which are imaginary for finite systems, may approach the real axis as the system grows. This was shown explicitly for the two-dimensional Ising model [3]. Having been assured that statistical mechanics can describe phase transitions, let us consider two simple cases, the condensation and freezing of spherical particles. The pair potential of a pure liquid of atomic particles is generally modelled as a steep repulsive potential, delimiting the particle diameter, and a somewhat less steep attractive potential, the two combining to form a potential well. When the particles are sufficiently dense they reside in the potential wells of their nearest neighbours, forming the liquid state in which the structure is largely determined by the repulsive cores. At sparse densities, very few reside in the wells at any given time, giving rise to the gaseous state.  At  low temperatures there is a range of densities in which the homogeneous phase is not stable. This is because on average the particles would be at a distance from their nearest neighbours which is close to the potential wells, and some of them must "slide" down into liquid state configurations - releasing room and energy for others to go into gas configurations far from the wells. At high temperatures, above the critical point, where the kinetic energy of the particles is large compared to the well depth, the distinction between these types of configuration disappears, and the liquid-gas transition vanishes. To obtain a liquid-solid transition we can use the same potential as above, but even a simple hard sphere "billiard ball" potential will suffice. As a liquid of hard spheres is compressed, the spheres become increasingly confined by their nearest neighbours. Forming an ordered crystal makes more efficient use of space, giving more room for each particle to move around in. When an attractive finite tail is added to this potential, then the two phases are of different energy, and this enters into the balance as well. At high temperatures, though, the potential well again becomes negligible, and all liquids  Chapter 1.  Introduction  4  approach the hard core model, as long as the temperature is not so high as to change the structure of the molecules themselves. The hard core potential, being infinite or zero, is not affected by temperature, so that fluid-solid transitions have no critical point. Any such freezing transition is thus related to packing efficiency, and it is of interest that a recent claim has been made [4] to have proven that the most efficient packing of hard spheres is in fact a lattice, as opposed to one of the infinite number of random packings possible. This is one of the outstanding questions of mathematics, and its solution would be a step toward solving the more difficult problem of proving that a lattice arrangement is the lowest energy configuration of particles bearing an attractive potential. This brings us to perhaps the most difficult problem of crystallization, namely, spontaneous symmetry breaking, or how the particles form the lowest free energy structure.  All theories of  freezing to date, as well as the present theory of crystallization, are designed to assess the possibility of a phase transition once a solid symmetry is selected. The ab initio prediction of the correct symmetry for any given interparticle potential is certainly beyond our present capacity. Having introduced the transitions of condensation and freezing, we consider what statistical mechanics requires in order to describe them. We have mentioned the Ising model above. This is equivalent to the lattice gas model of fluids. In this model particles are placed on lattice sites, and the repulsive core is modelled by allowing only one particle per site. When the attractive tail is allowed to extend only to nearest neighbours, the partition function can be solved exactly in two dimensions, and numerically to high accuracy in three dimensions. A somewhat different model is obtained if we represent the repulsive core as a hard sphere, and the attractive tail is approximated as infinitely weak and long-range. Ornstein [5] showed that the van der Waals equation of state can be obtained from the canonical partition function of this model, so we refer to it as the van der Waals model. Let us consider how each of these models, lattice gas and van der  Chapter 1.  5  Introduction  Waals, performs in the prediction of phase transitions [6]. The van der Waals model is accurate for describing dense fluids where particles have similar environments, in terms of nearest neighbours, etc. so that the potential they feel is well approximated by some constant background potential while the structure of the liquid is essentially determined by packing of the hard cores. At very low densities it can also be useful as the background potential is a good description of dilute gases. It does predict a liquid-gas transition, but at the critical point where individual particles experience widely differing environments, the van der Waals model predicts incorrectly both the value and dimension dependence of critical exponents. It is exactly here where the lattice gas model is most useful, as it treats accurately the attractions between particles, and this gives an excellent description of critical phenomena. At higher densities where packing is important the lattice gas is not at all descriptive of a real continuum system, and cannot give a liquid-solid transition. On the other hand, Longuet-Higgins and Widom [6] have used the van der Waals model to successfully calculate the freezing of noble gases, and although this is not a first principles approach applicable to general systems, it demonstrates the appropriateness of the van der Waals model at high densities. Thus statistical mechanics requires an accurate model of the repulsive core in systems where packing is important, and an accurate model of the attractive tail when the particles are in a wide range of environments. Returning to Figure 1.1, we consider what statistical mechanics would require to describe the crystallization transitions of an electrolyte solution.  Because the liquid  and solid phases are dense, the repulsive core must be treated accurately, and so any sort of lattice model is inappropriate.  However, simply because the system is dense  does not mean that all particles will have similar environments, as was the case in hard sphere freezing. Because there are three kinds of particles, a number of configurations would be possible around any given particle. In addition, the solvent has an anisotropic potential, thus again increasing the variety of environments each particle might find itself  Chapter 1.  Introduction  6  in. Because of these factors, we expect a van der Waals model to give poor descriptions of the liquid and its phase transitions. In our study, then, we will require a model which represents both the repulsive core and attractive tail properly. We could not expect to solve the partition function for such a model, but fortunately, this is not necessary. The information that will be of interest to us can be obtained from pair distribution functions of the liquid, and these can be obtained from integral equation theories. Liquid theories based on integral equations have been developed extensively during the past few decades, and the more accurate versions give good descriptions of non-critical fluids, including electrolyte solutions [7-9]. Also, during the past decade, a first principles density functional theory of freezing has been developed [10,11] which makes use of the results of liquid theory. This has been applied successfully to the freezing of hard spheres, and other systems, but has never been applied to crystallization from electrolyte solutions. This theory, as mentioned earlier, does not predict equilibrium structures, so that we must consider reasonable choices and compare their relative stability. The variety of possible crystal configurations is more extensive for this system than for simple hard spheres, so in general we shall only consider experimentally observed crystal structures. One feature not demonstrated in Figure 1.1 is the limit of stability of the fluid phase. Experimentally it is known that the liquid may exist for a period of time in regions beyond the phase transition lines. We should be able to use statistical mechanics to investigate such metastable states, and to do so, let us consider the following. Equilibrium is achieved in a system which is subject to certain constraints, such as walls defining its volume, or a bath defining its temperature. Equilibrium would no longer exist if the constraints were removed or changed. For instance, if we have a container with a 2:1 mixture of H2 and O2, this is not normally considered to be at thermodynamic equilibrium. However, by neglecting the possibility of chemical reaction, we may calculate the properties of the mixture as if it were an equilibrium mixture of diatomic particles, and these results will  Chapter 1.  Introduction  7  agree with experiment for a long period of time. The activation energy barrier provides a generalized "wall" which allows the mixture to be at equilibrium until a catalyst is added. Similarly, when a liquid is supercooled or a solution is supersaturated, although it is unstable with respect to a solid state, there is a free energy barrier which prevents formation of the solid until a critical nucleus is formed. Thus a metastable fluid can exist within a local free energy minimum for periods of time which are long compared to the time of measurement. The pair distribution functions we use are calculated from nonlinear integral equations which can have multiple solutions, not all necessarily of the same symmetry. One of these solutions corresponds to the equilibrium phase, and the others to metastable phases. In practice, when the integral equations are solved for an isotropic liquid, the isotropic symmetry must be imposed on the solution beforehand to make the calculation tractable. If this same symmetry is imposed past a liquid-solid coexistence line, we obtain a solution which does not describe an equilibrium state, but corresponds to the unique and reproducible metastable liquid. This procedure ceases when a state is reached where the homogeneous phase is locally unstable. Such a point defines the limit of stability. By this means we should be able to find the metastable regime of the liquid, and to examine the behaviour of the liquid near the limits of this regime. The object of this study is to apply the methods of statistical mechanics to determine the crystallization transitions of model aqueous electrolytes, and to examine the high concentration limits of stability. In Chapter 2 we review the density functional theory of freezing and thermodynamic stability theory. In Chapter 3 we formally show that all of the theoretical tools we shall need - liquid state theory, density functional theory of crystallization, and stability theory - can be derived from the grand canonical density functional theory of inhomogeneous systems. In Chapter 4 we review integral equation liquid theory, and introduce the models that will be employed. In Chapter 5 we make  Chapter 1.  Introduction  8  contributions to the density functional theory of crystallization, and apply this to the phase transitions of aqueous alkali halides. In Chapter 6 we derive the explicit expressions required for a statistical mechanical theory of stability for electrolytes, and use these to study solution behaviour near the stability limit. In Chapter 7 we summarize the principal results of Chapters 4, 5, and 6 along with the insights they provide into the phase transitions of electrolyte solutions.  Chapter 2 P h a s e Transitions  Let us consider a phase transition such as crystallization from an electrolyte solution. Salt is added to a pure solvent, and eventually a concentration is reached where the solution is in equilibrium with a crystal. Experimentally it is possible to proceed past this point to create a metastable supersaturated solution, but only up to a certain concentration limit. As this limit is approached, the system becomes unstable, and can no longer exist as a homogeneous phase. To locate these two points, the equilibrium coexistence and the stability limit, is the first object of our research. To locate coexistence we extend the density functional theory of freezing to crystallization transitions, in preparation for which we review the freezing theory in Section 2.1. To locate the stability limit by a statistical mechanical theory, we first require an understanding of the principles of thermodynamic stability, which we discuss in Section 2.2.  2.1  T h e D e n s i t y Functional T h e o r y of Freezing One of the most difficult problems in the field of statistical mechanics is to formulate  a theory of freezing from first principles [5], and during the past ten years or so most effort has been devoted to density functional methods. Since this area of study is relatively new, we give here a selective review of this work, with the object of assessing the viability of a theory of crystallization for model electrolyte solutions based on hard sphere plus embedded multipole potentials.  9  Chapter 2. Phase  Transitions  10  Density functional theory was developed to deal with problems of inhomogeneous fluids, such as the electron "gas" in the presence of nuclei, as explored in quantum mechanics [12,13], or a liquid in an external potential, as studied in classical statistical mechanics [14,15]. The theory of homogeneous liquids is in a relatively advanced state, and this can be used to advantage, for instance, by expanding correlations in the inhomogeneous phase about their isotropic counterparts. This fundamental approximation has its origins in the work of Kirkwood and Monroe [16,17], who formulated a theory of freezing in the canonical ensemble which replaced the angle-averaged radial distribution function of the solid by that of the liquid at the same temperature and density. This did not produce spectacular results, but continued to be recognized as a useful concept. Ramakrishnan and Yussouf [18] were able to obtain much better results by studying the problem in the grand canonical ensemble, and representing the direct correlation functions of the solid, c(r1,r2;p(r)),  by those of a liquid, c(r 12 ;/9[,), at the same temperature and chem-  ical potential. (r t is the position of particle i and ri2 = r 2 — r 1? and we have shown an implicit dependence on the single particle density, which is p(r) for the solid and pi, for the liquid.) Haymet and Oxtoby [19] reformulated the Ramakrishnan-Yussouf method as a density functional theory, increasing its flexibility (to study interfacial structure, for instance), and expressing the theory as it is commonly used today. In Chapter 3 we derive this theory in a general way, resulting in Eq.(3.78), but for the purpose of this discussion, we present the simplest form of the theory for a single spherical component. This requires the approximate functional  PMl = Jdr  | p ( r ) In ^  - [p(r) - pL]^-\  J dr_x J dr2c(r12] pL)[p(Ll)  - pL}[p(r2)  -pL], (2.1)  where /3'1 = kT, the Boltzmann constant times the temperature. This is referred to as a second-order density functional because the density is present to second order in the final  Chapter 2. Phase  Transitions  11  term. The two terms are generally referred to as "ideal" and "interactive" contributions. As will be shown in Chapter 3, AO is a functional of p(r) which becomes equal to the difference in grand canonical potentials for the liquid and the solid [—(pV)soii,j + (P^Oliquid] w h e n minimized with respect to p(r) for a given liquid structure, (p is pressure, V is volume.) The resulting p(r_) possesses the same (//, V, T) as the liquid, where p is the chemical potential. In a practical calculation, the liquid direct correlation function, c(ri2', PL)-, is obtained for some density, pi, and this, along with a specified symmetry for the solid, constitutes the input for the theory. The functional is then minimized with respect to p(r), or rather, its parameters. Because the minimum value of the functional is the grand canonical potential difference between the solid and liquid at equal (//, V, T ) , the liquid density is varied until the grand canonical potential difference so obtained is equal to zero, thus satisfying the requirements of equilibrium coexistence, which are msolid _ T-iliquid  solid  P p  liquid  =P solid _  —p  4  ,  liquid  The identification of the functional minimum with the grand canonical potential is only true for the global minimum, and so it is necessary to examine all reasonable solid structures in this manner in order to say with confidence what the true equilibrium solid is. Local minima which are not global minima are generally associated with metastable solid structures. As mentioned previously, this functional is only approximate. An exact expression can be obtained by making the substitution (see Section 3.3) |c(r12;pL) —  / J d\(l  - X)c(Ll,r2;px),  (2.2)  Chapter 2. Phase  Transitions  12  where px — pi, + X[p(r) — pi]. This new quantity is referred to as a coupling integral. (Note that the approximate form is regained by letting the direct correlation function at each point of the path of integration be replaced by that of the liquid.) To obtain improvements on Eq.(2.1) via Eq.(2.2), there are two general approaches one can take. In one approach [10], the quantity c^,^;  px) can be expanded as a functional Taylor  series in powers of A about pi. The coefficients are the second- and higher-order direct correlation functions of the liquid, and an exact result is given by the second-order expression, Eq.(2.1), plus the correction terms °° 1 f - S -U r J d<r-vn=3 -  t / * n C ( n V i - • • • >r«;PL)[P{T.\) - PL] • • • [p(rn) - pL]J  The third-order term, or portions thereof, have been obtained for some simple systems, but fourth and higher-order terms are unknown. The other approach [11] is to replace c(Ei>I.2) P(L)) by c(ri2j Poyt),  a  liquid direct correlation function at the same temperature,  but some other average density, popt- A number of theories based on this concept have been implemented, though none use the form of Eq.(2.1). Instead, they use a functional which gives the Helmholtz free energy, and then employ a double tangent construction to give equality of chemical potentials. These are all somewhat more complicated than the second-order approach given here. of choosing p0pt>  an  They differ from each other in their method  d are classed under the name of nonperturbative theories.  This  distinguishes them from perturbative theories based on the assumption that properties of the solid can be usefully expanded about those of the coexisting liquid. In applying any approximation, there are different methods which have been explored for evaluating the integrals. The first method [20] is to expand the periodic solid density in a Fourier series, so that the solid parameters consist of the average solid density and a subset of the infinite set of Fourier coefficients. This allows a reduction of the double integral in Eq.(2.1) to a sum over reciprocal lattice vectors. An approximation to this  Chapter 2. Phase  Transitions  13  method is to assume harmonic motion of any particle about its lattice site [21]. This constrains the Fourier coefficients to all depend on a single parameter, the width of the Gaussian distribution. It also suggests a third method of evaluation - expanding p(r) in direct lattice vectors and integrating in r-space [22]. Thus we have the following densities: 1. The general Fourier expansion. This is given by  p(r) = ps+ Yl "ne-"'21. The { a n } are independent expansion coefficients, and the {kn} are reciprocal lattice vectors of the specified crystal. ps is the average density of the crystal. 2. The Gaussian approximation in reciprocal space. This is the same as the above item, but an is given by a(|£n|) = ^ e x p ( - ^ e 2 / 4 ) , where e is the width of the Gaussian distribution. 3. The Gaussian approximation in real space. This is given by  Ft  where {FQ are the direct lattice vectors of the crystal. We now consider some of the general trends that have emerged from a number of studies during the past decade, and assess the possible usefulness of density functional theory for the present problem. Most effort in density functional freezing has been devoted to the study of single component systems of spherical particles, especially hard spheres. In order to freeze, a particle must possess a repulsive core, and we may classify as "sharp" repulsions those which decay faster than r - 6 , and as "soft" repulsions those  Chapter 2. Phase  Transitions  14  which decay as r - 6 or slower. This division is motivated by computer simulation results [23] which show that soft spheres freeze into a body-centred cubic (bcc) lattice, while sharp repulsions give rise to a face-centred cubic (fee) lattice. The latter is close-packed, so that sharper repulsions give more strongly first-order transitions, with a larger change in volume. It is naturally the goal of density functional theory to predict these trends, which it has not yet been able to do in a fully satisfactory way. The second-order theory (Eq.(2.1)) is quite capable of giving a transition to an fee lattice for sharp repulsions, such as the hard sphere and Lennard-Jones potentials, but does not give correct results for soft spheres [24]. Even in the successful fee transitions, the stability of the liquid is overestimated, giving liquid and solid freezing densities which are too large [24]. Nonperturbative theories are generally able to improve upon second-order theory in terms of quantitative agreement with experiment for systems of hard particles, but do no better with soft spheres [25]. Returning to the perturbative approach, but taking it to thirdorder, one obtains rather different results [26,27]. Because of the difficulty of obtaining triplet correlation functions, only a few studies have included this term in its entirety. These have shown [26] third-order theory to be successful at predicting the freezing of soft spheres into a bcc solid, but to also be worse than second-order theory at predictions for harder spheres. The implication seems to be that the perturbation expansion is not rapidly convergent [27], but that cancellations among higher-order terms result in some truncations being successful for certain classes of systems. Various discussions have been given of the importance of the triplet correlations in producing the non-close-packed bcc structure, but at any rate, these contributions to the functional must serve to destabilize the fee solid, even if it is the desired structure. These matters are apparently not fully understood, and the present understanding has not yet been generalized to a wide variety of other crystal structures. The similarity of results for second-order and nonperturbative theories reminds us that Eq.(2.1) can be considered as an approximation,  Chapter 2. Phase  Transitions  15  not only to the full perturbation expansion (by second-order truncation), but also to the coupling integral (by substituting c(rur^]  px) —> c(r\2',PL))- On the other hand, the  third-order theory can only reasonably be considered as a perturbation truncation. Thus, the second-order theory of Eq.(2.1) does play a unique role, and is the logical starting point for the study of any new system. Finally, concerning the expansion of p(r), the Gaussian approximation is found to have a nearly negligible effect on the thermodynamics of hard sphere freezing [21]. Similar results are obtained by both Fourier space and real space methods. When working with hard sphere systems, though, the discontinuity results in a long-range structure factor, so that considerable care must be taken to see that the Fourier expansion converges [21]. The expansion in r-space, on the other hand, converges with just a small set of lattice vectors [22]. Considering the success of the second-order theories in predictions for systems with sharp repulsions, it is not surprising that all studies of mixtures have used this class of potentials. These results are of interest to us since most of our calculations involve mixtures. Studies of binary hard sphere mixtures have examined the role of diameter ratio and shown that, as the diameters become increasingly disparate, the temperaturecomposition phase diagram changes from one showing simple behaviour to one showing a minimum and then to one showing eutectic behaviour.  Both second-order [28] and  nonperturbative [29] calculations show this progression, but the nonperturbative methods give better agreement with simulation. The second-order theory goes through this progression too quickly, as it also does for Lennard-Jones potentials [28] with parameters of real systems. Here it performs well for particles of similar diameter, but as the size difference increases it can give a eutectic while the actual system is at the previous stage of the progression. In such cases, the agreement is still good for small concentrations, but worsens considerably near the centre of the phase diagram. We conclude that secondorder theory is capable of giving reasonable results for such mixtures, but the accuracy is  Chapter 2. Phase  Transitions  16  eroded as difference in diameter increases, and as concentrations of the two components become similar. All of these mixtures described above have size ratios between .85 and 1. For very small size ratios, the ability to evaluate density functional results decreases, as the accuracy of the input liquid correlations becomes very poor. This brings us to one of the principal criticisms of most density functional calculations thus far. Nearly all have used a simple theory for the liquid structure, such as the mean spherical approximation (MSA), although a number of more accurate methods are available. This prevents us from knowing exactly how much error is due to the density functional theory itself. MSA solutions are appealing because they have a known analytic form for a number of systems, but it is desirable as well to develop density functional methods which can be applied to a variety of systems, many of which have solutions only in numerical form. This concern applies as well to two other classes of systems that have been studied, namely, nonspherical potentials, and charged hard sphere mixtures. The studies of these two systems [30-34] have the following points in common: 1. Relatively poor liquid state theories were used. 2. Few computer simulations of the freezing transitions are available. 3. Interparticle potentials were crude approximations to those of real systems. 4. Many have potentials with both hard and soft contributions. 5. All have used second-order density functionals. 6. Most give some reasonable transitions. 7. All have used only fc-space evaluations. It is precisely the results for these systems which would be of the most interest in guiding the application of density functional theory to model aqueous electrolytes, but there is  Chapter 2. Phase  Transitions  17  little that can be concluded from existing studies. The first three items make it difficult to know how much error can be ascribed to the density functional theory itself. The combining of both hard sphere and soft electrostatic potentials makes it unwise to use experience with hard and soft spheres as a means of guessing the reliability of the second-order method. These systems must then be considered in their own right, and it is therefore encouraging that any transitions at all are obtained. Some of the structures obtained are NaCl and CsCl salt structures for the charged hard spheres [30], and orientationally disordered close-packed solids [31,32] for nearly spherical ellipsoids and dumbbells. (We do not deal with the extensive literature on liquid crystal transitions.) The only orientationally ordered crystal that has been correctly predicted so far is the freezing of water to proton-ordered ice Ih [33], a non-closed-packed structure. This is a positive sign for this study, even though it was obtained by a very different method [34] than the one employed here. Ref. [33] will be discussed further in Section 5.3. Turning to the final item above, one reason for A;-space evaluations is that the direct correlation functions for systems with long-range potentials are themselves long-range in r-space, which results in a conditional convergence over the direct lattice vectors. In this work (see Section 5.1.4) we show the correct way to deal with this problem and present a simple method of performing r-space calculations which is much more accurate than fc-space methods, at least in this system. For completeness we add a few additional perspectives expressed in the literature. The most pessimistic view is expressed in papers by Lovett [35] and by Lovett and Stillinger [36]. The concern expressed in these papers is that only exact correlation functions will result in unique predictions by the density functional theory. This appears to be more of a formal concern, as all density functional calculations on real systems so far have, if they predict freezing at all, predicted clearly what the equilibrium solid is, and it is always a reasonable one with physical significance. Therefore, these views have not  Chapter 2. Phase  Transitions  18  prevented density functional theory from continuing to be widely used. A more obvious concern, expounded by proponents of the nonperturbative theories [37], is that the significant difference between second- and third-order theory predictions should make any perturbative approach untrustworthy, as the convergence is far too slow. As we have seen, the nonperturbative theories in general are only successful for the same systems as second-order theories (although they are generally more accurate than second-order for those systems), implying that if second-order theory is inappropriate, then nonperturbative theories may be so as well. In a recent study Rosenfeld [38] has attempted to explain the puzzling success of hard systems and failure of soft systems when treated by nonperturbative methods. Essentially it is argued that the large volume change on freezing makes the calculation insensitive to details, while the small volume change for soft systems amplifies the effect of small errors in the representation of solid correlation functions. Although pessimistic in tone, this study allows us to rationalize from a practical viewpoint the use of second-order or nonperturbative methods on hard sphere systems, and, presumably, on systems with similar features. A final point of view is that only a few "well chosen" reciprocal lattice vectors are necessary for the calculation of freezing properties. This appeals to the Verlet freezing rule [39] which says that the height of the first peak of the structure factor of a simple liquid coexisting with its solid is quite uniform for many systems, having a value of ~ 2.85. Thus, the first peak in the liquid structure factor is important in characterizing this transition. The lack of convergence of the lattice vector sum at this early truncation, it is argued, is only apparent, and is due to careless truncation at arbitrary numbers of vectors. A number of studies [40] have claimed accurate results using this approach, but it has not been widely accepted, as its application is not standardized. However, the idea that a few well chosen terms of the reciprocal lattice expansion should give good results (just as two or three terms in the perturbation expansion can give good results) is appealing, and it may yet prove useful  Chapter 2. Phase  Transitions  19  as a "rule of thumb" theory of freezing. Of similar spirit are studies which use a full lattice vector set at second-order, but a few carefully chosen third-order contributions [41]. W h a t does this all imply for the prospects of a density functional theory of crystallization for aqueous alkali halides? The first conclusion is that a second-order theory can be expected to give physically reasonable results, and can be realistically attempted. Because of the approximation of solid correlations by those of the liquid, it is important for the success of second-order theory that the environment of molecules have a certain degree of similarity in both phases. Even though alkali halides are highly soluble in water, prediction of crystallization will be testing the limits of this requirement. An advantage for us in these calculations is access to relatively accurate liquid correlations. This is important since in some cases diameter ratios differ considerably from unity, a condition which decreases the reliability of less accurate liquid theories. We also conclude that at this point it is not useful to use a more sophisticated density functional theory, especially since simulation results are not available to compare the results of various theories. It is against this background that we develop a second-order molecular density functional theory, as described in Chapter 5.  Chapter 2. Phase  2.2  20  Transitions  Thermodynamic Stability Theory Our approach to the stability of solutions is somewhat different from that generally  taken in the literature, so we give a careful discussion of thermodynamic stability to introduce it. The essential results of strictest criteria may be found in Refs [42,43] and the method of presentation is similar to that of Ref. [44]. We will derive certain conditions which must be met for a system to be locally stable. A system may be locally stable with respect to infinitesimal perturbations and still not be globally stable with respect to all finite perturbations. The realm of thermodynamics traditionally applies to global stability (equilibrium), but because we describe only local stability, our expressions are equally valid in the metastable regime, if it is assumed that thermodynamic functions remain well defined. To discuss thermodynamic stability it is useful to introduce thermodynamic potentials, such that the system is stable when at the minimum of a potential. Each potential is appropriate to a particular set of conditions. For instance, when entropy, S, particle numbers, { M } , and volume, V, are fixed, then the system is stable when energy, E, is minimized. Thus energy is a thermodynamic potential for this case. If instead, the entropy is not fixed, but temperature, T = (•gg)  , is held fixed, then by the second  law of thermodynamics we can show that E — TS is the potential to be minimized. This is a Legendre transformation which changes the independent variable from S to T, as shown by the differentials  dE = TdS - pdV + J2 VidNi, i  d(E - TS) = -SdT  - pdV + Y^ VidNi, i  where /u,- is the chemical potential of species i. Generalizing this method, let us consider  Chapter 2. Phase  21  Transitions  a thermodynamic potential of the form R = R(N;XUX2,  ...;y1,y2,...)  = E-  y1Y1 - y2Y2  .  We have chosen {x{, yi} as intensive variables, {Xi, Yi} as extensive variables, and  {xiXi},  {yiYi} as conjugate pairs, meaning that a;,- = (§x~j- We can define r = R/N where x,- = Xi/N, throughout  =  r(x1,x2,...;yl,y2,...),  the extensive density.  (For concreteness,  this section with a thermodynamic  the Helmholtz free energy of a simple binary  potential  f = F/N =  is mole fraction,  of particular  usefulness  contact to us,  mixture  F = F(N; N2,V,T)  where the conjugate pairs are {T,s}  we will make  = E-  TS,  f(x2,v,T),  (s = S/N),  {p, v} (v = V/N),  where the favoring of index 2 is  and {fi2 —  p-i,x2}(x2  arbitrary.))  Because a system is stable if its potential is at a minimum, a process in which {Xi}, {yi} are held fixed is spontaneous if SR < 0 for the process. (Similarly, cess in which {N2,V,T}  are held constant is spontaneous  a pro-  if SF < 0.) Let us partition  a system at equilibrium into two portions, containing NA and NB particles respectively, and adjust the system so that iV^, NB,X^, to 8NA + 6NB = 8Xf  Xf  change to new equilibrium values, subject  + 6Xf* = 0. The surroundings are the same before and after so  that all intensive variables are unchanged. Any such modification of the system must result in SR > 0 since { X , } , {y,} are fixed, and, since after removal of the partition, the reverse process is necessarily spontaneous. We can see what effect this restriction has on the derivatives of the potential, which appear in the Taylor expansion of SR, SR  =  (NA + SNA)(r  + SrA) + (NB + SNB)(r  + SrB)-(NA  + NB)r,  (2.3)  Chapter 2. Phase  22  Transitions  " e ' dr\„  A  1 J* J* / d2r  *•- - E(£>*fl + 5EE(^)W&*' + - ,  (2.5)  where n e is the number of independent extensive densities. In the derivatives all independent variables are implicitly held constant except the one with respect t o which differentiation is performed. Because {a;,-} are the same in both systems, the {Sx{} are zero and the Taylor expansion is only over the ne extensive densities for which constancy is expressed as 0 = (NA + 6NA)(xi  + 8xiA) + (NB + 6NB)(xi  + Sx{B) - (NA +  NB)xt  or 0 = (NA + 6NA)6xiA  + (NB + SNB)SxiB.  (2.6)  We have made use of SNA + SNB — 0. Using Eqs (2.4)-(2.6) we find that the zeroth- and first-order terms of Eq.(2.3) cancel and we are left with  2  NA +  SNA^yydxidxjJ  In discussing local stability, SxtB may be as small as we wish, so that higher-order terms are neglected and we obtain the quadratic expression  0<EE(^)^x^.  (2.7)  (For our example we have  Algebraic manipulation of Eq.(2.7) gives the following restrictions for the first few cases of n e-  Chapter 2. Phase  23  Transitions  Case ne — 1: ' d2r dx{'  >0.  Case rae = 2: Completing the square we may obtain either d2r ' <dxldx2> ( d2r\ ' d2r \ <9x a 2  , ( d2r n > 0 and 2  aX2  >0  (2.8)  >0.  (2.9)  ^xx2  or  <92r  a2r ,dx 2  \  2  , ( d2r n > 0 and dxS  dx\dx2 { d2r\ Ux 2 2 J  K  Case n e = 3: Proceeding similarly we obtain ' d2r \ <9xi  , / <92r > 0 and <9x 2 2  d2r ' dxxdx2/ ( d2r\  >0  3xx2  and d2r dx1dx2J (  d2 2  ' <9 r' dxj  d2r dxidx3 ' d2r  21  ,<9x 2 <9x 3 /  d2r \ \dxxdx3j 2 d r\ dx{' (  d2r ' d2r 2  ,dx2  >0,  dxxdx2 \fPr_ dxi'  or any of the other five permutations of the indices 1,2,3.  2.2.1  Strictest Stability Criteria  We have shown that certain expressions involving second-order derivatives of thermodynamic potentials are indicators of stability, and that their positivity is a criterion  Chapter 2. Phase  24  Transitions  of stability. We next show that the conditions for the case ne = 2 may be rewritten as  '0V  > 0 and *2  or  d2^/ dx22  >0  (2.10)  >o,  (2.11)  Xl  . /dV'N d2 r \ > 0 and dx 2 2 X l 5Xl2  2J 2  where if dr = Xxdxi + x2dx2 + • •  (2.12)  then we define the Legendre transforms r' = r — £1X1,  r" = r — x2x2, dr' = —X\dxi + x2dx2 + •  (2.13)  dr" = +xi dx\ — x2dx2 +  (2.14)  To prove Eqs (2.10),(2.11) we use standard identities and Eqs (2.12)-(2.14) to give '82r2„/ dxj  Xl  d dx2  'dr^ dx2/  _d_ dx2  ' dr ' dx 2  d  xi  Xl  Xl  + dxi _d_ dx\  Xl  dr  ' dxi  dx2  Xl  X2  'dr\  9x2j ' dx-i  Ml  Xl  x2  Xl  dx{  dXlJX2\dx2JXi  d2r  ' ay 2 dx2  Xl  d2r dx2dxi  dx2dxi I ( d2r dxi'  x2  2  8r 2  ' dr dxj  Xl  dx2dx\ ' 82r dxr2  (2.15) x2  We proceed similarly for r" and thus obtain agreement between Eqs (2.8),(2.9) and Eqs (2.10),(2.11). We note that because Eq.(2.10) and Eq.(2.11) are equivalent, they must  Chapter 2. Phase  25  Transitions  both be satisfied or violated simultaneously. (We can use the same example as in the previous section to demonstrate  that &  2  2  'd f'\  fd f\  \dv2)T \  / T,M2-Ml  !LV  \dvdx2).  (d2f\  \dv2) T x  ' T,x2  I  '  i. 1  WJT<V where f = f - X2(H2 ~ (J>i), df = —pdv — sdT + (/z2 — fii)dx2, df = — pdv — sdT — x2d(fi2 — Hi)-) By comparison with Eqs (2.8),(2.9), we see that  (Q^-T)  > 0 is a criterion of stability  not only in the case where the intensive variable x\ is held constant, but also in the case where the extensive density variable xi is held constant. Eq.(2.15) that for a stable system (g^-r)  < (9Xri)  since \dxrs)  We can also see from is also a stability  indicator and must be positive. Let us define a strictest stability indicator to be the second derivative of a thermodynamic potential with respect to an extensive density, where all the variables held constant are intensive variables. By extension to the case of ne = 3 and higher, we can show that no matter which variables are used to define a system, the strictest stability indicator will always be positive if the system is stable. Satisfaction of this criterion will in turn ensure the positivity of all similar derivatives with any extensive densities held constant, which we will refer to as weaker stability indicators. When a system becomes unstable, at least one of the strictest criteria of stability must be violated, while the weaker criteria may or may not be violated. Thus the simultaneous positivity of all strictest stability indicators is a necessary and sufficient condition for a system to be stable, while the positivity of weaker indicators is only a necessary condition. (For the simple binary mixture we obtain the following  relationships  Chapter 2. Phase  between stability  26  Transitions  indicators:  Thermal  '0V ds2  > < v,x2  \ds\ 'd2t'\ ds*  P>*2  «>M2-Ml  'd2h' > > ds2  >0 P.M2-M1  Mechanical  \dv\ 'd2e'\ > < dv2  dv2  T,x2 f/y  2  > > SyfJ.2-^1  &f_ dv2  >0 T,H2-t*l  Compositional  d.  Jbn  dxl  > <  dx\,  } > T,v  dxl  >0 T,p  where e = E/N, h = H/N,g = G/N are the internal energy per particle, the enthalpy per particle, and the Gibbs free energy per particle, and we have defined e' = e — x2{^2 — /^i) and h' = h — £2(^2  —  related to the familiar [dv?)  an,  d (al^)  /^l)- The quantities  (g^f) and fg^") for single components are  isochoric and isobaric heat capacities (Cv,Cp)  are  re e  ^d  to the adiabatic and isothermal  and the quantities  compressibilities  (XS,XT)-  To similarly express the strictest criteria of thermal and mechanical stability for a binary mixture, we would need to define an isobaric heat capacity with constant chemical  potential  difference, CPilt2-lil,  potential  difference,  Xr,«-/ii-)  and an isothermal  compressibility  with constant  chemical  Chapter 2. Phase  Transitions  27  Eq.(2.15) and its analogue for r" may readily be combined to give (d2r'\ \dx22) - 2^ ^ dr \  fd2r"\ \dx1y - ---^-. f d2r  \dx22JXi  (2.16)  V^iVxa  If the two numerators are strictest stability indicators, then this equation limits the possible relationships between different types of instability, eg. thermal, mechanical etc. (A pertinent  example of this for the binary mixture is the relationship between  mechanical  and compositional  stability of a binary mixture.  It is often stated (eg. Refs [44,4$])> that  a binary mixture  must become compositionally  unstable before it becomes  unstable.  The rationale behind this is that, because (9^2 )  indicator, then by specializing Eq.(2.15)  mechanically  is a mechanical  stability  as  (dV\ _(dy_\ \dxi)-\d ' 2 / T,pxi)^ \ ^"2JT,V  \dx2dv) (dy\ '  U  v-U)  dV2 J rj, '  it appears that before (g^£)  can reach zero, (g-f)  neglects the fact that the system  to very different conclusions,  7^ 0. All  to  (d2g\  dv2 I „  >2/  unstable with ( ^ 2 )  = 0. By specializing Eq.(2.16)  fd2f'\  we can see that consideration  must first become zero. All this  may be mechanically  that is required instead is that (g^f)  I ,X2  • M2-/*l  T f l  of the strictest  \d  V  i  '  1,P  (2.18)  \dxVT,v criterion  of mechanical stability brings us  which we briefly address now for general  systems.)  Each side of Eq.(2.16) can take on a value from 0 to 1. Therefore, if we allow the numerator on the first side to approach zero, the various mathematical behaviours of the system are:  Chapter 2. Phase  Transitions  28  1. Each side approaches zero. This means that on the second side either the numerator is vanishing, or the denominator is diverging. 2. Each side approaches a non-zero value. This means that the first denominator is vanishing at a rate comparable to its numerator, and that on the second side, if either quantity is vanishing or diverging, the other one is also. 3. Each side approaches 1. This is actually included in the above item, but with the stronger condition of equality between the numerator and denominator We are interested in the physical behaviour of the system, and not just the mathematical possibilities above. Simultaneous violation of both stabilities is consistent with the first possibility. If, for some reason, one of the stabilities is not violated, then the second option holds. The third choice occurs only when there is some simplification or symmetry in the system. A simple example is the ideal gas, as discussed at the end of Section 6.1. Another example is an equimolar binary mixture in which the self-interaction of both components is identical. A series of studies was performed [46] in which (g^y)  was  used as the indicator of mechanical stability. As we have shown, this is not the best choice, but a number of the results obtained were for this symmetric mixture we have mentioned, so that ( ^ 2 )  becomes the strictest indicator for those results. In Chapter  6 we shall give expressions for the quantities in Eq.(2.18) which will allow us to discuss such physical phenomena further.  2.2.2  Instability and Integral E q u a t i o n s  Integral equation methods are an important part of liquid theory [47]. They provide the input for our density functional work, and they are also used in determining the thermodynamic derivatives which define stability conditions for the solution.  Integral  equation techniques have been used previously to study mechanical stability of pure  Chapter 2. Phase  Transitions  29  liquids [48] and compositional stability of binary mixtures [49], often with particular interest in the critical point. Ionic solutions differ from these studies in that the phases investigated are not unstable with respect to other isotropic liquids, but to an ordered solid. This difference influences our route of calculation. Because integral equation theories are approximate, different routes of calculating thermodynamic quantities are inconsistent. When calculating both coexistence and stability properties for a given system, it is appropriate to use the same route for each. For instance, in a number of studies [50], demixing transitions of binary liquids were investigated using free energies obtained by the energy route. Stability conditions were then evaluated by numerical differentiation of these free energies. When density functional theory is used to describe an order-disorder transition, thermodynamics are calculated by a procedure analogous to the compressibility route. The compressibility route for liquids involves integrals over thermodynamic derivatives obtained by Kirkwood-Buff theory [51], the same derivatives as we have seen in a number of expressions in the previous section. It is therefore appropriate for us to use Kirkwood-Buff theory to calculate quantities related to mechanical and compositional stability. The thermodynamic derivatives of Kirkwood-Buff theory [51] are obtained from expressions for density fluctuations in the liquid. As one nears an instability these fluctuations can grow rapidly, yielding correlations which become infinite in range.  This  results in a divergence of the transformed correlation at k = 0. For a single component system such fluctuations give rise to divergence of the compressibility. In the case of the liquid-gas critical point, both this and divergence of the heat capacity have been well characterized by critical exponents. That both \T and Cp should diverge simultaneously  Chapter 2. Phase  Transitions  30  is at least consistent with the familiar thermodynamic expression CP  XT  Cv  Xs1  (2.19)  which we recognize as an example of Eq.(2.16). These effects are closely related to the solution of integral equations, all of which involve the Ornstein-Zernike (OZ) relation M r i2) = c(r X2 ) + pj  h(r13)c(r23)dr  where h(r) is the total correlation function. The OZ equation is a convolution, and may be Fourier transformed to yield  In performing integral equation calculations of XT (the other quantities in Eq.(2.19) are not accessible by a grand ensemble Kirkwood-Buff method), its divergence, and the breakdown of the numerical method, are ensured in the limit 1 — pc(k = 0) —> 0 by Eq.(2.20) and the relations P~XPXT  = 1 + ph(0) = [1 - pc(O)]" 1 .  (2.21)  The last equality is a restatement of Eq.(2.20) for k — 0. In a simple binary mixture we now have the possibility for compositional instability in which fluctuations are described by  MO), M O ) --+00, / h 2 ( 0 ) —• — 0 0 .  This expresses the notion that like components experience an overall average attraction to each other (a negative potential of mean force), while unlike components feel an  Chapter 2. Phase  Transitions  31  average repulsion (positive potential of mean force). These quantities are given by the Ornstein-Zernike relations [47] 1+Pi~hn(0)  = [l-P2C22(0)]/Dc,  (2.22)  I+P2MO)  = [l-piCum/Dc,  (2.23)  = c12(0)/Dc,  (2.24)  = [1-Pi2n(0)][l-P2^2(0)]-/Wi2(0)2.  (2.25)  h12(0) Dc  To have Dc —* 0 is the obvious mechanism to produce long-range correlations and also interrupt the numerical solution of the integral equations. As will be discussed in Section 6.1, we can use Kirkwood-Buff theory and the relations above to derive the results  fd2f'Y *VT«-  Sc l+xlX2VS/Dc  W  (PgY  £c  =  &r?L_ 2/T,p  (2.26) (2.27)  Sc'  Sc  = 1 - p(x2cu(0)  + x22c22(0) - 2x 1 x 2 c 1 2 (0),  Vc  = P2(c22(0)-c12(0))-Pl(c11(0)-c12(0)).  (2.28) (2.29)  where the asterisk indicates that the quantity has been divided by its ideal counterpart. Suppose now that as we approach an instability we have Dc —> 0, since we have shown that this is one possible mechanism, and assume Sc, Vc remain non-zero. From Eqs (2.26) and (2.27) we can see in this case that both compositional and mechanical instability occur. However, we may similarly derive the weaker indicators of stability,  T,x2  2  'd f\*  \»  2,  -  D,c,  -^(l + x1x2V2/Dc).  (2.31)  Under the aforementioned conditions Eq.(2.30) would not indicate mechanical instability and Eq.(2.31) would not indicate compositional instability, which is another illustration  Chapter 2. Phase  Transitions  32  of the importance of using the strictest criteria of stability. In chapter 6 we will develop such relations explicitly for electrolyte solutions, and apply the results of integral equation calculations to the determination of stability limits. We address two issues before concluding this section. The stability theory derived above is appropriate to systems at equilibrium. Points on a spinodal are not generally at equilibrium, so that we are necessarily making an assumption in applying the theory to our system. We are assuming that thermodynamic functions exist for metastable, homogeneous liquids. The integral equations impose an isotropic symmetry on the system, and it is reasonable that this would, beyond coexistence, correspond to a supersaturated ionic solution, as argued in Chapter 1. The second point is that for liquid-gas transitions of pure liquids it has been shown [52] that approximate integral equations do not always give a true spinodal, but in certain cases and regions only "diverge" to a finite value. To be more precise, there are a number of lines we may consider. There are the lines defined by x " 1 = 0 as calculated by the various routes, and there is a line beyond which solutions to the integral equations cannot be obtained. As shown above, the line defined by X? 1  =  0 f ° r the compressibility route is also a line at which numerical solution of the  integral equations is interrupted, but in some approximations this line can be preceded by another line where the compressibility is large but finite, and past which solutions do not continue. This requires very exacting procedures to discern, and our calculations have not been performed with an eye to distinguishing such cases. However, in the liquid-gas transition studies referred to [52], it is clear that divergent behaviour is being approached, even if it is not actually reached. We therefore make an assumption that the same is true of these binary mixtures, and that the results calculated will be faithful, though not exact, indicators of the behaviour of these systems.  Chapter 3  G r a n d Canonical T h e o r y of I n h o m o g e n e o u s S y s t e m s  In this chapter we lay the formal groundwork for the applications of Chapters 4 through 6. The essential results of this chapter are not original, and can be obtained from various sources. Our purpose is simply to show the theoretical unity of the methods employed in this work by deriving all formal expressions from grand canonical density functional theory. This is not to say that some elements could not be derived by other means.  The essential results are Eqs (3.67) and (3.68), whose application is briefly  reviewed in Section 4.1, Eq.(3.78), which is extensively developed in Section 5.1, and Eq.(3.84), which is a simplified version of the results to be obtained in Section 6.1.  3.1  Grand Canonical D e n s i t y Functional T h e o r y Density functional theory was initiated by the work of Hohenberg and Kohn [12],  and extended to the grand canonical ensemble by Mermin [13]. Evans [14] has given a review of the general theory and some applications to statistical mechanical systems. His notation is followed here in part, although we generalize to mixtures of non-spherical particles. We shall consider a mixture of n species, which possesses Ns particles of the s species, and N = Yls Ns- The generalized coordinates and momenta of the i are given by q. = (r,-,0»)  an  d p = (p  particle  •,Pn)i where 0 4 denotes the Euler angles and  p .p^ . the linear and angular momentum of particle i. We will use the vector notation to represent a set of the scalars for each member of the species, eg. N_ = {N\, N2,.. •, Nn}. 33  Chapter 3. Grand Canonical Theory of Inhomogeneous  Systems  34  The kinetic energy is given by KE, the configurational energy by U, and the energy of an external potential by $ . The central quantity in the grand canonical ensemble is the partition function [47], E = Tr c l exp[-P(HK  - £ • jV)],  (3.32)  where the classical trace is, for particles with d degrees of freedom (3 for spherical, 5 for linear, 6 for non-linear),  Trd =  l ? \.. Nn\h^ JV--^ Nl N  Jr<fe-^  and HN  U + KE + $  (3.33) n  u  dv---'iN)+Ei:  Ns  „2 Pr,i  +, x,y,z V>  2 P0,i  +  ? it <•<*>.  (3.34)  where /i is Planck's constant, and ms and 7fc>s are the mass and moments of inertia of a particle of species s. The partition function is also the normalization constant for the grand canonical distribution function, / = H-1exp[-^(//iv-/i-iV)], which is used to obtain macroscopic properties as averages of a microscopic operator 0,  (6) = TvdfO. Two examples of this are the one- and two-particle distributions N,  ps{g) = (53 % - £ , ) ) = (ps(o)), t=i  Chapter 3. Grand Canonical Theory of Inhomogeneous  Systems  35  which are related to the macroscopic observables of density and scattering structure factor. Another average which can be trivially performed is  {HN-H-N + p-^nf) = =  (HK-£•  N_-HK  +£•  (3-35) g_-kThiZ)  -kTlnZ  = n(v,i», which is just the grand canonical potential, identical to — pV. Let us define a functional  W  = Trcl [F{HK - E • N + 0-* In /»)],  which is identical to Eq.(3.35), except that /  n  is any distribution, not necessarily equi-  librium, which has the same normalization as / , namely, Tr c ] / n = Tr c i / = 1. We can show that f)[/ n ] is minimized by / . ^[/n]-tt[/]  = Tvd{(r-f)(HK-H,.N)  = Trd [(/» - / ) ( % - /£ • jV + r = -r  1  r1(r\nr-flnf)}  + 1  In /)] +fi-lrFrd(/n In / n - /» In /)  In 2(1 - 1) + ^ T r c l ( / n l n ^ ) - 1 + 1  = ^ T r c l { / ( f l n £ - ^ + l)}. One can readily show that (a; In a; — x + 1) is non-negative for all non-negative values of x, and that it has a global minimum of zero for x = 1, i.e., /  n  = /.  Using these results, it is our object to show that each single-particle density implies a unique external potential. Since /  n  is a functional of $ , then this would mean that  quantities Tr c j fO are functionals of the single-particle density. Suppose we have two distinct single-particle potentials, <bA / 3>B which give rise to fA, fB and to p^ig), £_B(q)Then, letting (• • • ) A , (• • -)B refer to averages in systems A and B respectively, we have  Chapter 3.  Grand Canonical Theory of Inhomogeneous  nA = (Hl-i£-N + = Sl" +  Systems  36  r1^fA)A  ($A-$BB-)  Similarly we can show that ClB <QA + ( $ s -  $A)A.  Adding these results gives 0<TrcI{(/A-/B)(<DB-<^)}.  (3.36)  Now we note that  Tr c l /EE^(£,) 5=1 t= l  n  =  Tv s=l  Na  .  dfJ2J2jdis(a.-ii)Mi) t'=l  = Jdqp(q)-l(q),  (3.37)  which, combined with Eq.(3.36) yields  0 < J dq [^(q) - £{$} • [<f (£) - <£%]. Obviously we cannot have fA = fB or pA(q) = pf(g)-  Therefore each one-particle density  corresponds to a unique $ , and we can consider / and f2[/] as functionals of p_(q). Using the same procedure as in Eq.(3.37) we can write ()[/] as « [ / ] = / dqpiq) • (l(q) - //) + Tr c l f(KE  + U + / T 1 In / ) .  Chapter 3.  Grand Canonical Theory of Inhomogeneous  Systems  37  More generally we can show that 0 [ / n ] is a functional of p^(q) where £n(<7) is not necessarily the equilibrium density associated with 4>. First we define p^(q) = Tr c ]  fnps(q)  and we will then show that we can write  tt[/n] = / dqEn(q) • {±{q) - E) + ^ n ( £ ) ] ,  (3.38)  where = Tr c l [fn(KE  f\£(s)]  + U + /T1 ln/n)].  The first term is explicitly a functional of £^(q)-  The second term, -T7^11 (<?)], has no  explicit dependence on $ , but is determined by / n , which is the equilibrium distribution / ' for some other potential $ ' . This in turn is a functional of some other equilibrium density p_'(q) which we identify with the non-equilibrium £^{q)- Then, defining ws(q) — Us - <f>s(q), % n (<7)]  =  -jdqf?{q)-w(q)  > Tvd[f(HK-li-N  + F[f?{q)]  +  (3.39)  l3-1lnf))  = «l^(a)], so that 0[^ n (£)] has a minimum at fP(q) = p{q)- We can now give an interpretation to the quantity F[p]iq)], which is generally called the intrinsic Helmholtz free energy. The reason for this can be seen by writing  %(<?)] +/%£(£)•/£ = -pV + G = F  (3.40)  = r\eti)] + /dg4!ia)-£(a),  (3.41)  where we see it contains all free energy contributions which do not come directly from the external potential.  Chapter 3.  Grand Canonical Theory of Inhomogeneous  Systems  38  The essential result of this section is that there exists a functional of the nonequilibrium single-particle densities which is minimized by the equilibrium single-particle density, and whose minimum is equal to the grand canonical potential. This result is of considerable importance in the remainder of this chapter.  Chapter 3. Grand Canonical Theory of Inhomogeneous  3.2  Systems  39  A p p l i c a t i o n t o Pair Structure D e t e r m i n a t i o n Using the technique of functional differentiation with the results of Section 3.1 we  can obtain expressions for the distribution functions [14,53,54]. Differentiating Eq.(3.39) with ws(q) fixed we obtain T-^TT = r  7  - ws(q).  (3-42)  Since ri[^ n ] is a minimum for p_(q), the equilibrium density corresponding to w(q)., we can write  sn [Pn]  Wis)  =o=  8:F[ ]  As discussed in the previous section, T\p^\  £  - ws(q).  (3.43)  does not depend upon w(q), but is fully  determined by p^_(q), so that its functional derivative is not affected by holding w(q) constant. Therefore we have p_  '8 * *  (3.44)  Ps{q)  Using the result from Section 3.1 that w(q) is a functional of p_(q), we can now use Eqs (3.43) and (3.44) to write  «T0 Sps(q)tpt(q')  _««-.(*)  (3.45)  &pt{£)'  fl[pn] is also a functional of w.(q), and we can again use Eq.(3.39), this time to obtain the expression  ^ ( i ) = ? / *. ftMg" ( W ( i J - ""(i'J - '• <2).  (3-46)  For pf(q) — pt(q) the first term vanishes and we can also write  m\d  Sps(q)  6ws(q)6wt(<£)  8wt{^)"  Writing ps(q) as  Mi)  T r d { [ E f % ~ £,)] exp[-P(KE  Trcl exp[-P{KE + U-U  + U - E ? E * u>«(j.)]}  £f' ^(fl,-)]  (3.47)  Chapter 3. Grand Canonical Theory of Inhomogeneous  40  Systems  the functional derivative of Eq.(3.47) may be obtained as  rl  80 (a)  £§j  N  =  Nt  °  (E%-i)]E%-a,-)])-/».(a)^)  = <£ *(i-iW-s,.)) + M E % - i W - 4 ) > - * ( ^ ) (s,»V(t,j)  i 6  = P${g,SL) + sAg-i)Ps{q)-Ps{g)pt{q).  (3.48)  It is convenient to consider 82J7[pJ/8ps(q)8pt(cr') from Eq.(3.45) in an alternate form. If we write T = JF id + J**, where [14] Tid  = Tr c l / i d (A'£; + / ? - 1 l n / i d )  = p-1T,/da.P'(s)(M^p.(a)) - !)>  ( 3 - 49 )  and [75]  2  ,2  1 \ 2  ^-UmJ V^U^J U W U ^ ' ' l  z,s,  (350)  '  then  ^T[p]  8st8{c£-q)  8ps(q)8pt(£) C  ^  ~  pa{<£)  -cst{q,q),  = -8ps(q)8pt(qjy  (3.51)  ^  Using the definition of functional derivatives, we can formally write  *«VtfO= E / jj^*P*U)<V-  ( 3 - 54 )  Substitution of Eq.(3.54) into Eq.(3.53) and inspection of the result reveals that  M(g-g) = E y ^ ( g W ) M g l ) ^ .  (3-55)  Chapter 3. Grand Canonical Theory of Inhomogeneous  Systems  41  which combined with Eqs (3.45), (3.48) and (3.51) allows us to write / > it ) (£,£ / ) - Ps(q)pt(J) =  ~  ps(q)pt(q^)cst(q,g)  Ps(q)Y: Jr^(£,£")[^)(£",£)  - PS{<()PM')W,  (3-56)  or hst(q,j)  - Cst(q,q!) = E / ^ ( £ " ) c t r ( £ , £ " ) ^ ( £ " , £ ) % " ,  (3-57)  where ,(»)/  /x _ hst(q, q ) = h  Pst(q,q)  f  (3.58)  - 1.  Ps(q)pt{q') Eq.(3.57) is the inhomogeneous OZ equation. We can obtain a second relation for hst(q,<f)  and cst(q,<f)  by the generalization  of a method due to Percus [47] for isotropic liquids. This requires restricting UN = f/(o , • • • , £ „ ) to a pair potential, so that the presence of an additional particle may be treated as a one-particle external potential. Consider the grand canonical partition function for a system of particles with chemical potentials {/J.s}, interacting as before by c  ^v(£ 1 ) • • • J£AJ) within the external potential $ , but now also under the influence of an  additional particle of species s at q^,  Uo(%;qv---,qN)  = UN+l - UN = Y, E M £ t - -£o)> t  t=i  where ust{q — q[) is the pair potential between particles of species s and t. For an external potential UQ + $ , we have the grand partition function Es  =  Trdexp[-P(KE OO  =  E  OO  ••• E  + UN + [U0 + 1  .  71  *]-iiN)]  Nt  TTT—TT7 / I I I I *;&•) exp[-/9(lfr +  where the superscript s on S 5 denotes it is in the field of the (N + 1)  U0)]df,  particle, a particle  of species s. For convenience we have performed the integral over momenta and have  Chapter 3. Grand Canonical Theory of Inhomogeneous  42  Systems  denoted  <(&) = exp[/?(/^ - <Mi,))]Mt. Continuing, we obtain  1  ~  oo  oo  y ... y .  _^i Z  s\%o)  Af„=0-'Vl-  ^ ATi=0  n  :(&) Jv  Nt  /iIIKWexph^iv+i]^ n- "  t=li=l  S 1 * _ _ ~ _ N.z*a{%) NA...N \ z*(a ) ^ 1~> '*' 2 ^ " " 2 ^ n  Nt-6st  / n n'^-)^-/?^]^- 1 .  X  (3.59)  We compare Eq.(3.59) to the quantity /?s(£) of the system without the additional particle 1  oo  oo  1  Ly ... y -  Ps(g)  N, X -i  ~  /  n  E%-i) U=l  Nt  AT  / I I I I <(£••) exph/JC^^, • • -,qN)]dq;  J  t=\ t'=i  oo  oo  ~ l~f " ' 1-/ " ' 1-J -" Ni=0  Ns=l n  M \... M \  ATn=OiVl-  iv  «-  Nt~Sst  / I I n"<(2 i ) ex P[-W£i.-d J v- 1 .£)]¥" 1 '  X  (3-6°)  Comparing Eq.(3.59) and Eq.(3.60) we see explicitly that S s is dependent upon CL, and Similarly we may compare prs(q), the single particle density  that H*(^) = E,ps(q0)/z*(q).  of species s in the field of an additional particle of species r, -i  Pl(9)  =  oo  oo  AT ~*(n\  n  r  ^ E - E l v T ^ / n ^  N1=0  Nn=0 °- °  £*(<7 )  -ivl- ' ' '  JVn  °°  - ^  Nt—Sst  <=1  i=l  1 / n n'^dOexpM^+j^2-" N,' • • • N l °-°  sz;(q) /VN^*^/7^  iyi ~or\%) (a ) N^x=o '" N ^ 3 =1 '" N ^P ATS=\ JVnn=o -  i  i  °°  •=E-  M 2 o ) •=- AT1=0  <(£,-) e x p [ - W o + l / J V ) ] c ^ - 1  n  iV  /• "  «-  «=i  oo  E  ATS=1+S3t  E - NE=0 N,=l+S, t  n  Nt-Sst  i=i  JV.(JV,  - «„)*;(£)*:(£„) iVi. • • • JVn.  (3.61)  Chapter 3. Grand Canonical Theory of Inhomogeneous  to pyj(£,£'),  43  Systems  the two-particle density of species r and s with the additional particle  absent,  1  i  rf?(2.sf) = = ? E<—> - ETV-, I . . .  N I n  x  /  oo  1  oo  -  = E - £ ••• £ ••£ "JV1=0 X  Nt  £%-i.)%'-^.) J] II <(£,) exp[-/?^( £l , • • • ,qN)} Ndc£  AT s =l+5 st n  JV t =l+5, t  N.(Nr -  SrJzMz:^)  AT„=0  Nt—Srt—Sst  Jtl' fl ""z;(qi)ex1?[-pUN(qv---,qN_2,q,j)]dq:N-2  (3.62)  Comparing Eqs (3.61) and (3.62) we have PpX&Qn) pr(q)  (3.63)  Defining zrs*{q\ £Q) = ^(tf) e x p f - ^ u ^ ^ - ^ ) ] , we expand l n [ ^ ( £ ; %)/zrs*(q-, qj] as a functional of pKqjq^) about urs(q — £Q) = 0 to obtain In Z  T{SJ%) v - [(  z  *s(l) 6 34/  ?/  Srttjg-g')  r {q;, PPrs.(S}2o)  ~s v x ' i o ' '  6 In  tprt(q';q„W + «ri(£-90)=0  z:*(q-%)  \prtti;9o)-p*ti)W + ~<ZM) Urt(g-^)=0  Note that ^*(£,'£ 0 ) = exp[/?(/xs — <^s(£) — urs(q,£0)]/A:S  so that if we define u ; ^ ; ^ ) as  to s (£) in the presence of an additional particle of species r at a we can write In zrs*(q\ £Q) = flwrs(qj,q)  — l n A s . Then the second term in the parentheses in Eq.(3.64) is equal to  fitsfi(q — <f)/pl(q]qj0) — c L(£>2''£o)' ^(£,io) Priq^jpsiq)  ~ -°  an<  t  ^ * n e expression becomes  •/  ^(a'.ao) MfiJ  -  /><(</)  <v +  (3.65)  Chapter 3. Grand Canonical Theory of Inhomogeneous  In (hniz,^)  + l)  =  -Purs(q,%) -Purs(q,qX))  44  + J2 I Pi(^)cAl^)hrt(q^,qX))d^ J  t  =  Systems  + hrs(q,q0)-crs(q,q0)  If we neglect higher-order terms, then, for a given ps(q),  + •••  + ---.  (3.66)  Eqs (3.57) and (3.66) allow us  to calculate the pair structure of the system. This allows us to completely determine the thermodynamics of a system which interacts only by pair potentials. For a homogeneous system (ps(q) =constant) Eq.(3.57) reduces to the Ornstein-Zernike equation and a truncated Eq.(3.66) to the HNC (hypernetted-chain) closure of liquid statistical mechanics. When higher-order terms of Eq.(3.66) are approximated by those of a reference system, then it is known as the RHNC (reference hypernetted-chain) closure. An exact, non-perturbative density functional theory of freezing would require us to solve these equations for inhomogeneous systems as well as for homogeneous liquids. In practice this is very difficult, and, as discussed in Chapter 2, we shall make approximations that allow us to use only the pair correlations of the coexisting liquid state. As a result we will only require the isotropic versions of Eqs (3.57) and (3.66), which are  Kt(Li2,Q.l,£h)  - ^ ( Z j ^ f ^ O j )  = TTlYLpr  J  C  ir{L23M2^)hsr{Ll3,Q.\-,Q3)dr3dQ3,  (3.67) Mhst(Li2,0.1,0.2)  + 1] = -^u3t(r12,£Li,Qa)  + hst(r12,Q,1,Q,2)  ~ cst(r12,QiiQs)  H  • (3.68)  The solution of these equations will be reviewed in Section 4.1.  Chapter 3. Grand Canonical Theory of Inhomogeneous Systems  3.3  45  Application to Freezing Theory To develop the density functional theory of crystallization we require an expres-  sion for the grand canonical potential difference between two different states of the system, denoted A and B, in terms of the single-particle densities, pf(q) and pf(q)We have previously defined /362Jrex[pj/Sps(q)Spt(q^) pSJ^^pj/Sps^q)  = -c^\q).  = — cst(q,(f) and so we also define  It is also convenient to define p*(q) = pf(q) + X[pf(q) -  pf(q)]- From the definition of functional derivatives we have W*V]  8  -^^hxs{l)dq_.  = £/  (3.69)  Using basic notions of calculus we can write  =  dX[pf(q)-pf(q)),  (3.70)  and, noting that the functional ^jFex[/9A] is also a function of A,  Substituting Eqs (3.70) and (3.71) into Eq.(3.69) and integrating the result with respect to A we obtain  /w*V] - /w*V] = -E[1^f s  Similarly we can show that c ? W ) - c?>(9;/) = £ t  %4% zxM(a) - ps(i)}-  (3-72)  J  JO  [XdX  ldq'cst{q_,q'-p*)[pf(q>)  Jo  J  - p?(q>)].  (3.73)  ~  Substituting Eq.(3.73) into Eq.(3.72), and using Eqs (3.43), (3.44), and (3.49) to show that c a)  i (£;£) = HAsps(q))  ~ Pw,(g),  Chapter 3.  Grand Canonical Theory of Inhomogeneous  Systems  46  we can obtain  /?A^X[£J = E/^W(2)-MA^(£))}[^(£)-^(£)]  ~T,f0dx jf dx' j di\di = ^[p8]  where AJ^^p}  - J^x{p_A}.  (3.74)  c  st(i, ?'; £X'M (i) - p^i)M(i) - tf(£%  According to Eq.(3.39), we can obtain 0ACI by  adding Eq.(3.74) to PA^p)  = £  fdqpf(q){ln[Aspf(q)}  - 1} - £  J  s  s  / % ^ ( i ) { l n [ A s ^ ( £ ) ] - 1}  (3.75)  J  and t o d  - E \ s  l [Pf(3)fiwf  (2) - PiidPwHsL)],  (3-76)  J  which gives  s  J  + Ejdi  [pfisD^(^jfy) -ipfd)~Psd)}}  - E £ dX £  dX  ' J dl\H  (3-77)  c  st(q,i\ pX')[pf(q) - pA{q)][p?{<£) ' P^)]-  This is the fundamental result we desire. To consider how we might employ it, let us assume that wA(q) = wf (q), and that we have knowledge of an equilibrium state A, ie. we know pf(q),  and that we can also determine cst(q, (/; px) for any given state B. Then  if we create a new functional by replacing pf (q) by p^,B(q)  (we know we can do this from  the result of Section 3.1) then this new functional will have a minimum and b e equal to the grand potential difference for p^'B(q)  = pf(q)-  In general this difference will be  non-zero, but when state A is chosen such that minimization results in (3A£l[pJ = 0, then we have satisfied TA = TB,  pA = PB,  wA(q) = wB(q),  Chapter 3.  Grand Canonical Theory of Inhomogeneous  Systems  47  which are the requirements for the inhomogeneous states A and B to be in equilibrium. For vanishing external potential, the last condition reduces to the more familiar p,f = fif. The functional which we shall use is simpler in three respects than Eq.(3.77). One simplification is that our reference state is an isotropic liquid, so that pf(q)  = pf,  a  constant. Another simplification is the absence of any external fields, and setting the chemical potentials equal in the two phases. This results in wf(q)  = wf(q).  The final  simplification results from approximation of the coupling integral, which, as discussed in Section 2.1, may be done in two ways which give identical results. In the perturbative approach, c(q, q'\ £A ) is expanded about A = 0, and truncated after the first term. In the nonperturbative approach, c(q,(/;p_x ) is assumed to be equal to c(q,(/;p_A)  throughout  the range of integration, so that f0 dX f0 d\' may be performed trivially. We therefore have a functional, which, when minimized, is equal to  -[pf(r,n)-Pt}\  This functional is parametrized into a calculable form in Section 5.1.  (3.78)  Chapter 3.  3.4  Grand Canonical Theory of Inhomogeneous  Systems  48  Application to Stability From Eq.(3.48) we have -1 M ( j ) _  (2)  r 1 J ^ J y = p\V{q,q) + 6stS(q - £ > s ( £ ) - pa(q)pt(q!).  (3.79)  This gives the response of the density, £^q) to changes in the chemical potential or external potential. If the integral  (3-8°)  M(i) = £ / | ^ M ( 2 ' ) < V becomes divergent for any value of £, even when the {6wt(q')}  are as small as we wish at  every point <?', then this indicates that the system is unstable. For a single component, the inverse of this functional derivative, as shown in Section 3.2, is defined by the relation c(  ,x  S{q  q ) =  f M r ) M g " ) j /, dq  -- -  J 6W(qJ')8p(qJ) -'  f)in(ti\  fi( n — n'\  ,«,an  (3 81)  "  where  Spilt)  Mi)  -c(q,q').  (3.82)  Approaching an instability implies that the integral  ( 3 - 83 )  Mi) = / j^jWW  becomes as small as we wish at every point q even though p(q^) may be divergent at any point q'. This means that  % - j) Mi)  c  (?,i') ^ o  for all such £'. There is not much use we can presently make of this formal result for the solid state, but we may easily recover results for the homogeneous liquid in the <f>(q) —> 0 limit. With w(q) —> p,p(r)  6  —> p Eq.(3.83) becomes  "%-i')  "=I --c(k P  -<\q-q'\)  = 0) Sp,  Spdq (3.84)  Chapter 3.  Grand Canonical Theory of Inhomogeneous  Systems  49  so that 1 — pc(k = 0) —• 0 as an instability is approached. This allows us to make contact with the discussion surrounding Eq.(2.21). Similar relations are obtained in Section 6.1 for binary mixtures, both simple and ionic.  Chapter 4  Calculation of Liquid P r o p e r t i e s  Liquid correlation functions are the input for our density functional and stability work in Chapters 5 and 6. Because the methods for obtaining the liquid solutions are well established, we give only a brief review in Section 4.1, which by itself is insufficient to fully describe the procedure. These details may be obtained elsewhere [7-9,58,60-63]. The results of Section 4.2 have been previously published [55]. The initial results of Ref. [55] were obtained by Dr. D. Q. Wei and Dr. G. N. Patey. My role was extending the results to additional systems, and studying their effect in the context of crystallization (Chapter 5).  4.1  Integral E q u a t i o n Calculations of M o l e c u l a r Liquids A unique homogeneous liquid mixture can be defined by sets of pair potentials and  particle densities, and the thermodynamics of such a system can be obtained from the grand canonical radial distribution functions, gij(r12,Q.i,Q.2)-  These can be obtained in  the RHNC approximation by solving simultaneously the Ornstein-Zernike (OZ) equation Vij(r.i2iQi,Qa)  = g-^^Pk  J [riik(Li3,Qi,Q3) +  Cik(r13,ni,Q3)]cjk(r23,Q.2,Q3)dr3d£l3, (4.85)  and the RHNC closure,  ln[»fo-(El2, A l . i k ) + C , - j ( r 1 2 , O l , G 2 ) + 1] = % ( r 1 2 , O l , 0 2 ) ~ ^ M y ( H l 2 . 0 l , f l 2 ) +  50  B  ,  (4-86)  Chapter 4. Calculation of Liquid  Properties  51  Figure 4.2: Schematic for RHNC calculation of simple liquids  RHNC  T)ij(r)  -  OZ  FT-1  fjij(k)  where B represents the higher-order terms of a reference system.  In this work, the  reference system is a hard sphere mixture, and the higher-order terms are obtained from the computer simulation radial distributions of Verlet and Weis [56]. Eqs (4.85) and (4.86) can be identified with Eqs (3.67) and (3.68) using the relation  where the distribution function gij(r12,0.i,fi2) rfij(r12,0.1 >H2)  an  d Cij(rl2,0.1,0.2)  ma  1S  equal to hij(r12,0i,02)  + 1- Although  y be- considered as defined by the above equations  for the purposes of the calculation, but they can also be systematically derived through cluster expansion analysis [47,57], which sheds more light on their properties. The series function, Tjij(r12,Oi,0a)j  nas  the particularly useful property of smoothness at contact,  which allows it to be accurately back-transformed.  For a simple mixture of spherical  particles, the solution of Eqs (4.85) and (4.86) can be accomplished by the iterative scheme of Figure 4.2, beginning with some initial guess for Cij(r). Such straightforward application is not possible for non-spherical particles and long-range interactions, and we give here only a brief review of the procedures required. It is necessary to expand the  Chapter 4. Calculation of Liquid  Properties  52  correlation functions in terms of rotational invariants [58], e.g.,  <7.-i(r12,&,£2) = E/ mB 'Cft( r ")*r l (fii»fla.£")  (4-87)  mnl  where fmnl  is can be defined arbitrarily, but it is most convenient [58] to define it as  S< 2m + l)(2n + 1). The rotational invariant is given by ^'(a,fl2,r12)=  W ti'i/\'K"  where (  m i  "  "f ' )R^M)KM2)Ko(ti2) v  (4.88)  J  't J is the 3-j symbol [75], Rlmn(Q)  is a generalized spherical harmonic  as defined by Messiah [59], and r is the unit vector pointing in the direction of r. The general theory for such expansions, and their use in solving the OZ equation is given in a series of papers by Blum and Torruella [58(a)] and by Blum [58(b)]. The RHNC closure can also be expanded, but not in an obvious manner, due to the logarithm of the expansion. The correct method for doing this was given by Fries and Patey [60], along with a summary of Blum's OZ method. Using these techniques, one can find the projections c™J!^(ri2) by an iterative method based on the cycle given in Figure 4.3. The Fourier transform of the full correlation function [58] involves the Hankel transform of each of its projections. The zeroth- and first-order Hankel transforms are numerically feasible, and the purpose of the "hat" transform is to reduce all higher-order Hankel transforms to these simpler cases. In addition, when long-range electrostatic forces are involved, it is necessary to remove these in some manner to avoid divergence of t h e k = 0 value of the transform. We denote the short-range function with superscript s. The OZ equation for non-spherical potentials takes the form of a matrix equation, potentially of high order. The x-transform is a means of grouping projections so that sets of smaller matrices are obtained. After reversing the various transformations to obtain ^^"L(^i 2 ), the RHNC closure  Chapter 4.  Calculation of Liquid  Properties  53  Figure 4.3: Schematic for RHNC calculation of molecular liquids  , , ,  hat transform  short ,  range  .  ,  Hankel , transform  ^ni(s)(  fdr  mn,{s)  chi transform  czrw imnx(s)  drc^;ij\r)  RHNC  OZ  »x(»)/  d_ mnl ( \  JV£E W (*)  inverse chi transform  3r  fjmn .(r)  inverse hat transform  f)mn'\s>  -  long range  ** '  J  r)  mnKS ~mnl(s), <n > k)  -  inverse Hankel transform  Chapter 4. Calculation of Liquid  Properties  54  may be applied [60], where the derivative with respect to r provides the means of dealing with the logarithm of the expansion. For this method to be useful, the expansion in rotational invariants must possess sufficiently rapid convergence to allow truncation at a reasonable number of projections. It has been shown [61] that imposing n, m < 4 gives good convergence and is numerically feasible for this water model.  Chapter 4. Calculation of Liquid  4.2  Properties  55  M o d e l s of A q u e o u s Alkali H a l i d e s Details of extensive RHNC calculations for aqueous alkali halides are given in the  work of Kusalik and Patey [8,62]. In this section the basic model is reviewed, and a modification of the model is introduced. This is followed by a discussion of the resulting changes in solution properties, and their relevance to the study of crystallization. Perhaps the simplest molecular model of an electrolyte solution is that of hard spheres with embedded multipoles. The charged hard sphere (CHS) model of aqueous alkali halides used in this study consists of representing ions as hard spheres with embedded point charges, and the water molecule as a hard sphere with a point dipole and quadrupole.  The ions are considered non-polarizable, but polarizability of the water  molecules is taken into account indirectly through a modification of the dipole using the self-consistent mean field (SCMF) theory [63]. This involves an iterative process in which the initial-guess dipole is gradually changed until a dipole is obtained, which, along with all other multipoles in the system, generates a particular mean field. The mean field desired is one which interacts with the experimental polarizability of an isolated water molecule to result in the same dipole as the one giving rise to the field, to some level of convergence. The dipole obtained by the SCMF theory varies little with salt concentration for this model [8], so all finite concentration solutions use the same dipole as that obtained by SCMF calculations for the pure solvent. These are all at 25°C. This same dipole, which is ~ 50% larger than that for an isolated molecule (see Table 4.1), is also used for the pure solvent at other temperatures. This is an additional approximation, but as will be evident in Chapter 5, it will not affect our essential conclusions.  The  quadrupole moment we use for input is also modified, but with a different motivation. Initially, its off-diagonal elements and the sum of its diagonal elements are all zero, but the diagonal elements are individually all non-zero [65]. It is noted, however, that \QZZ\ is  Chapter 4. Calculation of Liquid  Properties  56  somewhat smaller than |<5xx|, |<5yj/|5 and so we set Qzz = 0 and QT = QXx = —Qyy This is referred to a tetrahedral quadrupole, and corresponds to a molecule of higher symmetry. It results in slightly simpler calculations, and has been shown to be an excellent approximation [61]. Summarizing, then, the CHS model consists of hard sphere and electrostatic potentials, where the parameters determining these are given in Table 4.1. The hard sphere  Property d  B.20  Value 2.80A  d  U  1.904A  d  Na  2.352A  dK  3.024A 3.584A 3.248A  4 fiz (gas phase) tiz (SCMF)  4.032A 1.855 x 10- 18 esu cm, Ref.  [64]  2.768 x 10"18esu cm 2.63 x lG- 26 esu cm2, Ref. [65]  ^ixx  **izz  -2.50 x l(T 26 esu cm2, Ref.  [65]  -0.13 x 10- 26 esu cm2, Ref.  [65]  2.57 x 10~26esu cm2  QT  Table 4.1: Input parameters for model electrolyte solutions, particle of species i.  di is the diameter of a  diameters for the ions are taken from x-ray data [66]. Particle densities for the integral equation calculations were taken from experimental values at 1 atm. Extrapolated values for the density were used for concentrations beyond saturation.  Chapter 4.  Calculation of Liquid Properties  57  For this model it was found that some of the calculated thermodynamic properties do not agree well with experimental data for aqueous solutions. An important example of this problem is the mean ionic activity coefficient, j±, obtained for aqueous alkali halide solutions. For these solutions it was found that the models considered in Ref.  [8] do  not deviate sufficiently from the limiting law at high concentrations and are generally in poor agreement with the experimental values. One possible way for improving agreement is suggested by previous calculations using McMillan-Mayer (MM) level theory [67,68]. The MM approach to electrolyte solutions does not treat the solvent explicitly, but incorporates average solvent effects into ionic potentials in a systematic way. Friedman and coworkers [68-70] and later Pettitt and Rossky [71] showed that the activity coefficients predicted by such theories were very sensitive to small changes in the effective potentials, and that very good agreement could be obtained by small, arbitrary modification of the potentials. In view of these earlier observations, it would seem likely that some relatively minor adjustment of the simple CHS ion-ion potentials used in previous calculations might give thermodynamic properties in much better agreement with experiment. If this is the case, then such adjustments could well have an effect on density functional predictions for transitions from these liquids. The calculations of Ramanathan and Friedman [69] and Pettitt and Rossky [71] suggest that our ion-ion potential requires some repulsive interaction in addition to the hard sphere term. Therefore, we consider models in which the solvent-solvent and ion-solvent interactions are retained exactly as in Ref. [8], but the interionic potential is modified to oo, u  ij(r)  = \  r < dij , ,  where Aij is a positive constant. Also, d^ — (di + dj)/2, where the hard sphere diameters, di, of the various ions are again those used in Ref.  [8]. Thus, the only change we make  to the models considered previously is the superposition of the Aij/rn  repulsion onto the  Chapter 4. Calculation of Liquid  Properties  58  charged hard sphere interactions of the ions. Following Ramanathan and Friedman [69] in all calculations presented here we have used n = 9 and Aij — Ad-  ,  where A depends only upon the salt and has the same value for all three ion-ion interactions. This model we refer to as the CHSR9 model. The motivation for this particular functional form comes from solid state theory which defines a unique value of A for each crystal structure of a salt [69,73]. However, since in the CHSR9 model the hard sphere interactions are retained, and since an investigation of sensitivity is our principal interest, we simply view A as an adjustable parameter. We did carry out a number of calculations with values of n ranging from 7 to 11 and found that our ability to "fit" experimental data was not sensitive to the actual power law employed. Also, it is worth remarking that although the additional repulsive term is included between all ionic species, its influence upon the solution properties is exerted mainly through its presence in the cation-anion interaction. Ions of like charge in the CHS model already have a negligible correlation in the near-contact region where the repulsive term is significant. In calculations it is convenient to express energies in units of kT and distances in units of the solvent diameter, ds. The strengths of the repulsive interactions are then determined by specifying the reduced parameter A* = pAlda =  pAiil{dZfidl,).  The ln7± results obtained for model NaCl, KC1 and Csl solutions with different values of A* are shown in Figure 4.4.  In all cases it is clear that the activity coefficients are  very sensitive to the value of A* and that curves which closely resemble the experimental plots can be obtained. For NaCl the value A* = 6.5 given by the solid state formula employed by Ramanathan and Friedman [69] produces results which lie very close to the  Chapter 4. Calculation of Liquid  Properties  59  Figure 4.4: The effect of short-range repulsions on the mean ionic activity coefficients of NaCl, KCl and Csl. The solid curves are experimental results and the dashed curves are for the present models. The straight lines denote the limiting laws. The theoretical curves are labelled with the value of A* used in the calculation. The concentration, c, is expressed in moles/litre.  Chapter 4. Calculation of Liquid  Properties  60  Chapter 4. Calculation of Liquid  Properties  61  experimental values up to concentrations of 3.5 moles/litre (M). However, for KC1 and Csl the values of A* (i.e., 1 and 0.7, respectively) which give reasonably good fits to the experimental data are much smaller than the parameters estimated from the solid state assuming power law repulsions, so the fact that the values agree for NaCl is likely just coincidence. For NaCl (A* = 6.5) and Csl (A* = 0.7) the added r - 9 repulsion accounts for approximately 5% of the total cation-anion interaction at contact, and for about 1% in the case of KC1. Thus relatively small perturbations can bring about such significant changes in the behaviour of In j±. Two additional comments concerning the results given in Figure 4.4 are in order. We note that for NaCl with A* — 6.5 numerical solutions to the RHNC equations were not stable up to the experimental saturation concentration (i.e., 5.41M). Therefore, on the NaCl plot we also include results for A* = 2 which is approximately the largest value of A* for which numerical results could be obtained up to the experimental saturation concentration. Also it should be noted that for the present model the dielectric constant of the pure solvent (93.5) is somewhat larger than the experimental value (78.5) for water at 25°C. This difference is evident from the theoretical and experimental limiting laws shown in Figure 4.4. If a model solvent with a dielectric constant closer to that of water were used, then slightly larger values of A* would be needed in order to fit the experimental data. However, this would certainly not influence our basic observations concerning the sensitivity of ln7± to the interionic potential. As would be expected, the added repulsion also influences ionic correlations, the most significant structural change being a decrease in the contact value of the cation-anion pair distribution function <7+_(r). This effect is most dramatic for NaCl and results for 2M NaCl solutions with different values of A* are shown in Figure 4.5.  We see that the  contact value decreases very sharply from 100.5 for A* = 0 (i.e., charged hard spheres) to practically zero when A* = 6.5. Again, this is particularly striking in view of the fact that  Chapter 4.  Calculation of Liquid  Properties  62  Figure 4.5: The pair distribution function g+-(r) for 2.0M NaCl with various repulsion parameters. The solid, dotted and dashed curves are for A* = 0, 2, and 6.5, respectively.  Chapter 4. Calculation of Liquid  Properties  63  CO U3 CO  o c<\ II II II * ^  # ^  * ^  -co 05 CN • 00 00  o or-4  CO  r-l  C^  •  II I  + i  + -c\?  iO  I  CO  C\2  Chapter 4. Calculation of Liquid  64  Properties  the contact interaction energy changes by only ~ 5% over this range of A*. The decrease in contact value is accompanied by an increase in the height of the second peak which results from structures where a cation and an anion are bridged by a water molecule (Note that the cusp-like behaviour occurs at a value of r corresponding to collinearity or, in other words, at the separation where all three species can no longer remain in contact). From the curves we would expect a decrease in the number of contact ion pairs and an increase in the number of solvent-separated ion pairs as A* is increased. A rough measure of these two numbers can be obtained by integrating /9_47rr2<7+_(r) (/>_ = p+ is the number density of negative ions) from contact to the first minimum, and from the first to the second minimum giving the results summarized in Table 4.2. The numbers Table 4.2: Average numbers of contact and solvent-separated ion pairs calculated as described in the text. In all cases the salt concentration is 2.0M. The solvent diameter d, = 2.8A and di/d, = 0.84, 1.08, 1.28, 1.16, and 1.44 for Na+, K+, C s + , C I - , and I", respectively.  salt  A*  contact  solvent-separated  NaCl  0  .735  1.20  NaCl  2  .247  1.43  NaCl  6.5  .019  1.59  KC1  0  .975  1.06  KC1  1  .618  1.23  Csl  0  1.38  1.08  Csl  .7  1.15  1.35  given in the table indicate a clear shift from contact to solvent-separated ion pairs and, in the case of NaCl, essentially the elimination of contact pairs altogether. Thus for the present model the strong deviations from the limiting law which occur as A* is increased  Chapter 4. Calculation of Liquid  Properties  65  are associated with increasing numbers of solvent-separated as opposed to contact ion pairs. For NaCl solutions it is necessary to essentially eliminate contact pairs in order to obtain good agreement with the experimental activity coefficients. Of course here we have looked at only the CHSR9 model, and there may well be other physical mechanisms which would bring about agreement with experiment.  Nevertheless, our calculations  do suggest that solvent-separated ion pairs may be much more important than contact pairs in determining the thermodynamic properties of at least some aqueous alkali halide solutions, such as those with small cations. To conclude this discussion, let us summarize the two main points in the context of crystallization. Firstly, ln7± is a key quantity in determining the concentration at which crystallization can occur.  Therefore one would expect the results of a theory  of crystallization to be sensitive to details of the model which strongly influence this thermodynamic quantity. Secondly, because our density functional theory approximates correlations in the solid by those in the liquid, it would be expected to perform better (for an approximate theory) if the local structure of liquid and solid have some degree of similarity. Therefore, when two models give drastically different local solution structures, as is the case with NaCl above, it is expected that this would strongly affect the choice of structures which density functional theory is capable of predicting as the end product of crystallization. These points will be explored more fully in the following chapter.  Chapter 5 Calculation of P h a s e C o e x i s t e n c e  As discussed in Chapters 2 and 3, we apply the density functional theory by parametrizing the single particle density, and then minimizing the functional with respect to these parameters for a given liquid input. Building on work such as Refs [22,31,32], we will derive the detailed expression necessary for r-space calculations of a molecular system, and apply it to various transitions involving model aqueous alkali halides. This is the first time that r-space expressions have been given for molecular systems. It is also the first theory for molecules of general symmetry which uses correlation functions expressed in a generalized spherical harmonic expansion. Some of the results of Section 5.1.4 and most of the results of Section 5.2 have been published previously [74].  5.1  D e r i v a t i o n of t h e M o l e c u l a r D e n s i t y Functional E x p r e s s i o n s From Eq.(3.78), the second-order functional to be minimized is 0ASI  =  E/^r|joLl(r,a)lnf^^^-[  / 3 l  (r,0)-  / 3 ! X  //Q]}  (5.89)  - 2 Y, J ^ i / <*& / dZ* / d£L2Cij(L12,£2X, 0 2 ) x [ # ( r i , f i i ) - piLIIu][pj{r2lQ.2)  ~ PJL/IQ]  ,  where ns is the number of species, /?;(r, Q) is the solid single particle density of species i, pn is the average density of species i in the liquid, ^ ( r ^ j O i i I ^ )  66  ls  the direct correlation  Chapter 5. Calculation of Phase  Coexistence  67  function between particles of species i and j , and 0. = Euler angles,  IQ  = 8TT2 for molecules of general shape,  fi = spherical polar angles, IQ = An  for linear molecules,  0. = 1,  for spherical particles.  IQ = 1  (5.90)  Crucial to the implementation of this theory is a suitable parametrization of the single particle density.  This could be done exactly in terms of a Fourier expansion  and generalized spherical harmonics, but such a procedure would require a prohibitive number of minimization variables. We use instead a density of a form which is physically reasonable, and whose implementation is feasible. Following previous workers [31,32], we impose the condition that spatial and angular variables be independent of each other by choosing a density of the form TTLi  Pi(r,a)  = J2Mr.i)9iP(ai),  (5.91)  p=\  where p is the index over the m, particles of species i in the primitive unit cell, and fip(Li),9ip(£Li)  are  both normalized to 1. The most common approximation for the  spatial density is a Gaussian distribution [22],  fM  = liTT-TE ex P[-fc - LB + I*D 2 /QP 2 ],  (5-92)  where {R} is the set of primitive lattice vectors, and r,-p gives the average location of the p  particle of species i within the unit cell. The Gaussian form of fip has been  demonstrated to be an accurate approximation for spherical particles in crystals [21]. Its efficacy in molecular crystals has not been similarly demonstrated, but it is the obvious first level of approximation. The parametrization of the angular distribution, gip(Q.i), will be dealt with in Section 5.1.8, where we will again assume harmonic motion and impose a roughly Gaussian form.  Chapter 5. Calculation of Phase  Coexistence  68  The theory we develop will be specific for c t j(r 1 2 ,ili,0.2) expressed in the rotational invariant expansion  co-O^k.fla) = EClM*)*™'^^^).  (5.93)  mnl fit/  In Eq.(4.88) $^"'(0.1,0.2) £12) w a s defined in terms of the generalized spherical harmonic, R'mniQ.), which for t h e case of n = 0 reduces to the complex conjugate of a regular spherical harmonic, Rlm0(Q)  = J-^^Y^iO^).  It possesses t h e orthogonality property  [75]  /'<fiiflL(fl)/Cn»(Q) = 2nT*«'«m«'^-  (5-94)  It is convenient to give the second-order functional as /?A0 = (/?A0)i + (/3Aft) 2 + (/?A0) 3 ,  (5.95)  where (/?A0) X = E / * / « l | M i : , f l ) l n ( ^ ^ )  (/9Aft) 2  =  -[pi(r,n)-piL/Iu}\,  I «! - - £ J dLl j dn, J dr2 J c?0 2 c tj (r 12 , & , fi2) x[-/»t(Ei,Qi)/»jL//n - Pj{Za>Qa)piLlIa + piLpjL/Iu\  (0Afl)3 = - ^ E J dr^ Jd^i Jdr^ JdS^Cij^^MPii^MpA^M5.1.1  (5. 96)  (5.97)  (5-98)  Ideal C o n t r i b u t i o n  We consider first the term which arises from various ideal contributions in the derivation of the density functional theory of crystallization in Section 3.3. T h e key assumption in what follows is that the Gaussian distributions of fip(Li)  are sufficiently narrow that  integrals pertaining to each particle may be calculated independently. This, along with  Chapter 5. Calculation of Phase  Coexistence  69  the previous requirement that the angular and spatial integrals be separable, results in the considerably simplified expression,  (/?Afi)x  PiL  =  V£\piL-pia-pia\n^ t=i  fdrfdfL +£ »=1  J2  fip(^)9ip(Ql) In  p=l  (5.99) 9=1  n  =  VY\\piLns  Pis - Pis In -j-  rni  +E E  / drfMlnfiM  + V / <&</,•,,(&) In #,(&  where the average solid density is /?,-5 = j*y fdrfdQ  ,(5.100)  />,-(r,f2). The elimination of the  sum over q is as follows. Let us suppose that the {e, p } are sufficiently small that peaks are non-overlapping. Then let us divide the integration over r into regions each containing a single density peak. If we consider a region containing a peak of the function fip{r_i), then within this region we have ln/o^rj,£2 X ) ~ ^[fiP(r1)gip(Q.1)].  This eliminates the sum  over q. The integral J dr / i p ( r 1 ) l n / , p ( r 1 ) is identical to that in the spherical case (e.g. Ref. [22]) and yields — Vpc[3/2 + 3 \n(y/ireip)] for each particle in the unit cell to give TTli  =  PiL  X) { P^ ~ Pc £ p / 2 + m T~ + h =i  \~h ~ h  3 ln  (v^etp) - Bip  (5.101)  P  Bip = J dQLgip{Qa)\ngip(Q.TL),  (5.102)  where pc = Pis/i^i is the density of unit cells. Once the form of gip(Q.i) is chosen, the ideal contribution is determined. If there are Zs symmetry related points in the primitive unit cell, so that it possesses only mi/Zs  crystallographically distinct particles of species i, then Eq.(5.101) may be  rewritten by substituting mi P=I  mi/Zs  p=i  Chapter 5.  5.1.2  Calculation of Phase  Coexistence  70  Interactive Contribution: k = 0  Using the orthogonality property, Eq.(5.94), we can write Eq.(5.97) as (/?Afl) 2  =  ~  £  PiLPiL^Zm  + ^ E / * i / ^ i / dr2 J dU2  Cij(Ei2,Qi>Qa)[Pi(Li,Qi)pjL  + Pj(ra,Qa)PiL]-  (5.103)  We will now show that only the average solid density, /9,s, contributes to this expression. As mentioned before, we can, if we wish, formally express our parametrized density as a double expansion in a Fourier series and in generalized spherical harmonics, Pi(L,a) -EEalabc^T-[RlM)r,  (5-104)  abc  fc  =  where [•••]* denotes the complex conjugate, and we note using Eq.(5.94) that c*oooo Pis/Sn2.  We substitute this into Eq.(5.103) and perform the angular integrations to  obtain V  ns  87T2  ^ E  nvv'  J  ;  V  ;)4,  m  o  ^  ns  m  ^^/^^^2c^(r  1 2  )e'^^^0(r  1 2  )  )  ~r  where we use the convention n = —n. We have also used the properties \m—n\ < I < m-j-n and p,' + v' + A' = 0 of the 3-j symbol (  m >  "  \ J in order to eliminate / and A'. Now  consider the quantity  Jdrx  Jdr_2cll%{r,2)eik--^R^{u2)\  =  J dLle^^  =  V8kQ J dr12rl2ctl{rl2)  =  iKtikpSnoS^oV I  =  /^12c^(r12)e^-^^0(r12)]  V6koSnoK'oc00<ij(0)  J  dr12R^0(r12)}  dr12rl2c™0i:j(r12)  Chapter 5. Calculation of Phase  Coexistence  71  where we have used translational invariance and orthogonality of both plane waves and generalized spherical harmonics. We may treat j dr^ / ^ £ 2 c ^ 0 c | ^ ( r i 2 ) e 1 - ' L l i ^ 0 ( r 1 2 ) the same way, which reduces Eq.(5.105) to (/9Aft) 2  =  - \ v £  piLPiL^niO)  8TT 2 2  n  °  ij=l  and dividing by V gives ' BACl\ —T7-  1 ns = - « E [(^  -  ?")(ftx - ft.) - Pispjs]^Zj(°)-  (5.106)  Chapter 5. Calculation of Phase  5.1.3  72  Coexistence  I n t e r a c t i v e C o n t r i b u t i o n : r-space  The final term can be expanded in the form  (/?Aft)3 = - ^ £ j dLlJ'dShj mnl  dr2JdQ.2  fi'l/'X' ^  '  m;  (5.107)  x J2 h fo )9iv (Ql) £ / « (E2 )#g (Qa) p=l i  g=l n3  mi  m  j  -iEEEE/w E i,j=l p=l g = l mnl  M'i/'A'  m  n  f  /i'  *'  A'  yu-rnn' ip  V-Tirnv'v jq  2 /,  2  p  ^e-(ra-\E'+Z}q])2/c]q2 -w 3 / 2 ,t . 3 Z 31 R'  X  (5.108)  ^'0(^12)5  where  w?/» = jdm™,M)9>vm-  (5.109)  We can use translational invariance in a number of ways to give a more convenient form for this term. For pairs of particles assigned to lattice vectors R,Rf, the contribution to the functional varies roughly exponentially as — \R — Rf\. For all significant terms, then, \R — R[\ is small compared to the dimensions of the crystal, and we may set Rf = 0 and multiply by Vpc, the number of primitive cells. We are also free to make the substitutions £1 —* Hi +Zjq, La —* £2 + 1 J V For brevity we define R_T = — R + Zjq — Zip, so that we have  (^)  = -fEEEEE^E R  i,j=l  1  x [Tttiptjqy  p=l ( J = l mnl 1X1/  m  .1  n  ..1  ' w I  Svynv'ttJIfni/'i/ ip 39  (5.110)  ti'iy'X'  j dLl J dr2 c™?,.(r12) e x p [ - ( £ l + Rrf  /eip2 -  v\/ejq2]Rlxl0(ru).  Chapter 5. Calculation of Phase  73  Coexistence  Similarly we choose r a as the origin for the integral over r 2 , and set r 2 = £12 + 1.1 to give ns  V V )  m,  rnj  -f£££££/w£  //?A0\ 3  Z  R i,j=l p = l 9 = 1 mnl  1 yTTciptjqJ roo  m  n  '  1 W.m'J''iH^'™J/'1/  jxVA'  f°° JO  x / drir2exp[-r2/e,-p2 - r2/eig2] ./o  x / dr_\ / dr_1 2 exp  v  2  (5.111)  exp  El-^r  7Li-r12  RUtu).  ^jg  We substitute into the angular integrand using an analogue of the Rayleigh expansion [75], e-R-r  £(-)'(27+ !);,(/&•) £ /=o  where im{x)  RUmiR'Ut)]*,  (5.112)  A=-/  is the spherical modified Bessel function, related to the modified Bessel  function and to the Bessel function by [76,77]  ll{x) =  (5.113)  V^'+l0,0'  1,/+A*) = H )l +' i+ 2H-; ^ i (**)  (5.114)  The angular integrals are then straightforward and result in the expression  (g*a\ \  V  m  _e±Y.t£££/""" £ (; ;, •,)*?'*#"  / 3  ^  .R « J = l p = l ? = 1 mnl  x7  ^ /  {•Ktip€jqf  J  V  '  ^ 1 2 ^ c^,.(r12) exp[-r22/ei?2 -  RT2/eip2}  (5.115)  JO e  t 00  '  /i'l/'A'  c?ri r x ex]  0  tp  T f-jq 2 C j p C. tjq 2  \  2 r i  (4-)% f ^  ) ii ( ^ :  ) ^o(4).  Combining Eq. 6.633#4 from Gradshteyn and Rhyzik [76] with Eqs (5.113) and (5.114) gives the identity  I  00  dx x2e~ax  2  ii{^x)ii{^x)  1  / 7T  = —w — exp -3 4 Va  /52 + 7 2  4a  «/  2a/'  (5.116)  Chapter 5.  Calculation of Phase  74  Coexistence  which when applied yields  \  V  / 3  R  »,J = 1 P=l 9=1 mnl  47T  I tip  x exp  M'l/'A'  ^  '  l-OO /-00  J  HT  - 2 e-  + Cjg 2 0  2RTr12  f\2  2 , e "„ 4-  | tj  f  . 2 i  \  (5.117)  -  . 2}  f  (ri2€i 2+R Cjq2)2  ity |<eexp ,. 2„p 2/cr . 2 , t 2\ It is convenient to introduce the quantity 1 c c c C »p  jg l »p + J9 ^  } and its reciprocal,  and to combine them with the two exponentials in Eq.(5.117), to obtain ns  \  V  ) 3  £  mi "V  «J=1 P=l 9=1 mnl  47T  fjVA'  ' ^ 2 ^ o  X  2i? T r 12 X exp  f  \  M  "  /  (ri2  exp  e L  #r) 2  2  ip  4- e \  2 /  J9  2RTrl2  . 2 i f . 2 T^ i *-7 *-jq /  H  *-ZP ^p  \  W(r*2)  (5.118)  2 i ,. 2 e. \Mp i \?g  The quantity e~xii(x) is convenient to work with, and has a finite series expansion [76,77]  e  (+ "M = ^ * £ .» l JL» {HM-fr-*} t 2x ^Jifc!(/-ib)!(2x) i/(  (5.119)  A  x r, x < i . (2/ + 1)!!' (The odd factorial is defined as (2/ 4- 1)!! = 1 • 3 • 5 • • • • (2/ + 1).) Some rearrangement allows us to write Eq.(5.118) in a form convenient for numerical computation, '/?A0\  ns  Pc  mi  mj  47T  V /•oo  x / ./o C(ijpqRT;  ( r ^ - i ? . ) 2 " C{ijpqR ;r ), T 12 . 2 _L 2  c?ri2r22 exp  r 1 2 ) = ^ Ci(ijpqRT; /  r12) exp  e c  T  tJ9  2#rr12 . 2 i ,. 2 Ctp T^ c J ? ^ f  (5.120)  f  tp  ' «/  2i? T r 12 e  2  4- e  \ 2  (5.121)  Chapter 5.  Calculation of Phase  Ci(ijpqRT-,r12)  =  75  Coexistence  j:c™>J(r12)  (5.122)  mn 1LV  The indices in Eq.(5.120) define a pair of particles in the crystal for which the integration over ri2 is to be performed. The function C(ijpqRT;  r\2) is a linear combination, specific  for the given pair of particles, of projections of Cj_,(r12, H x , Q_2). The integral in Eq.(5.120) is written as a normalized Gaussian times C(ijpqRr;  r\2), and as e;p, tjq —•» 0, the Gaussian  becomes a three-dimensional delta function selecting the value C(ijpqRT; of terms in Eq.(5.120) for which R = 0, i = j,p  RT)- The subset  = q depend on RT = 0. In this case  we see from Eq.(5.121) that only the / = 0 terms contribute (ii(0) = Sio), implying A' = 0,ra = n , / / = —v' as well. In Eq.(5.122), the quantity in square brackets reduces to just the rotational invariant, $™ n '(f^ p ,Q.j g , i£ T ), if <7;p($2i) is an angular delta function. Although the quantities c™™^(r12), W£"' M , W^v,  R[y0(RT)  may in general be complex,  the functional will, of course, always be real.  5.1.4  T r e a t m e n t of L o n g - R a n g e P o t e n t i a l s  For many systems of interest, the interparticle pair potentials are long-range, decaying as 1/r 3 or slower. Since asymptotically the direct correlation function varies as —(3uij(r12, On, Q.2), it will have similar long-range behavior, and the sum over lattice vectors in Eq.(5.120) will be conditionally convergent [78]. In this section we present a method of dealing with this aspect of r-space calculations. We first consider a simple case where all interactions are spherical, and the only long-range interaction is the Coulomb potential. Then the Eq.(5.118) is given as  Chapter 5. Calculation of Phase  (PMT\  =  76  Coexistence  - ^  _p£  4?r „  (r 12 - RTf  roo  Jo  c?r 12 ri 2 exp ( \  ce  2i? T r 12 e  e  ns  mi m>  «'p +  jg  . 2 i  tp  i  £c .  jg  „oou (  \  2  2#Tr12  \.  . 2 ^ip  e  /  (5.123)  i e. 2 ' ' I *-jq  47T K 2 E E E E 2 T^i P tl g t 1 [^(^ + ^ 2 )] 3 / 2 X  J  z-oo  ' drl2r\2 o  exp  (r 12 - # T ) 2 . 2 i Cjp T e  . 2 cJ?  £  000 r ''OO , i j ( 1 2 )  c  ip  ~  c  jg  4i? T ri2  (5.124)  The approximation made in the last step is accurate for RT > 0 and typical values of tip, Cjq. For 72T = 0 we have &o(0) = 1, but we will not need this result here. We define a potential, f»j(r), such that  0 0Vij(r)  =<  r < d. v  where d^ = (di + dj)/2, and introduce a short-range direct correlation function, c 0 0 t j ( r ) , defined by c  oo%(r) = cw°i'jir) - Pvij(r) •  The required integrals involving Coo°j(r) can be evaluated in two steps, using  V v ) The contribution from c00jj(r)  3  - {V )  + 3  { V )3  can be calculated as shown above, and that from /?u,j(r)  is done as follows. The contribution of f3vn(r) to terms where RT = 0 is negligible, since (3vn(r) and e x p ( — J T T ) do not contribute significantly in the same regions for typical values of e,-.  Chapter 5.  Calculation of Phase  77  Coexistence  Therefore the following will refer to terms where RT / 0, as indicated by a prime on the sums below. For RT ^ 0 we have  '0ASIV  n,  PC  = -*r\ fiE«,i=l E P=1 E 9=1/ E 2  ^ip  n,  \ A ( e t p 2 + e jg 2 )  itT  + e•J. 9 2  n-%  m, m> \ '  = - ^ \' E EEE £ «,j=l P=l 9=1/ 2  x^ip2  r(-)*«  (r 12 - # T ) 2 "  Jri2 exp /  i  mi rrij  + £jg2 I  dll-RT  2  e  2  \J*{Up + jq )  RT  exp(-x2)dx  •% E E E E ' 1it ^ K t'J=l p=l 9 = 1 /  T  PC  (5.125)  4  R  i,j=lp=lq=lj  XlT  The first set of sums is the reduced electrostatic energy per unit volume when the particles are fixed at lattice sites. For a simple crystal this can be conveniently expressed using the nearest neighbour distance, RNN,  an  d the Madelung constant, M, as —TPCM/RNN  [79].  M has been tabulated for a number of simple crystals of high symmetry. The second term gives the correction due to vibrational motion of the ions about their lattice sites, and falls off rapidly as a function of RT. To calculate the long-range contributions to the functional for more general interactions and crystals, we do the lattice summation in three steps. First the sum is calculated straightforwardly as in Eqs(5.120)-(5.122) for some truncated set of lattice vectors. Second, for each pair of particles included in the truncated sum, the interaction energy of the two particles held fixed to their average positions and orientations, 0Uij(RT, f^p, 0j ? )> is added to the first sum. Third, the energy of a crystal of particles fixed to their average positions and orientations is calculated by an Ewald sum method [78], and this  Chapter 5. Calculation of Phase  Coexistence  78  quantity is subtracted from the whole functional. This method again allows us to use a relatively small set of lattice vectors in the summation of an effectively short-range correlation function, and allows us to use the Ewald sum method over the long-range part of t h e correlation. T h e Ewald sum replaces the Madelung constant, since M is simply a number which may be calculated by methods such as Ewald summation for crystals in which the electrostatic energy varies linearly with the density. We use the Madelung constant method for the anhydrous alkali halides study described in Section 5.2. T h e other calculations, Sections 5.3 and 5.4, use the more general Ewald method. We apply the Ewald method to sum all multipole interactions in the crystal. Although some of these (dipole-quadrupole, quadrupole-quadrupole) are unconditionally convergent, we find it t h e most convenient approach. The explicit expressions for our Ewald sum calculations are [78,80] rpelec  JV  3  4  v71" ,-=i  z?self  2-v 3  N  V 7 1 " ,-=i  . pk-space  =  1  V  , rir-space  N  A<y5  ,-=i  45  / r i oc\ N  v ^ ,-=i — —  — (5.127)  £k-sPace  =  2TT ^ -L e ~ fc2/47 iw  {  (5.128)  • N  N  i  cos  Y^qi s i n ( i • n) + YliBi'k) " N  +  J29«  cos  N  (£'&) ~ £(/£i'&)  .1=1  ^r-space  =  n ) - o 53(4.-Qi-k) s i n ( i • r,)  (k'  sin  J JV  ( £ • & ) - o 12(k-Qi-k)  »'=1  /<£ JA  E  ^ ^ R , ij)  ->2  JV  cos(k • Zi)  "* t = l  + £%££(& i,j) + K ^ ^ R , i,j) (5.129)  V E. id ) + £ d i P P d i p e G £ ' *J) + C  ^SSS(&i,i) = «?i/o(r«),  S  p  (  ^  *'•?') + £ q u a d , C q u a d ( ^ > * , » ,  (5-130)  Chapter 5.  Calculation of Phase  Coexistence  79  ^S^dfrC&i.J)  =  Mrij)[Qi(N-!Uj)-qj(jii-!Uj)l  (5.131)  *EE£d(&».»  =  f2(ri])[qi(rirQl.rij)  (5.132)  =  -Mrij)^-H~  =  -2{f2(rij)(Hi-Ql-Lij-H-Q±-rij)  Eti^aipi&hJ) ^SSp(&i,i)  + qj(ruj.Qj:.Lij)}/3, f2(rij)(i^-rUj)(H-Lij),  +/3(ry)[(£i-&i)(&i -Ql-Uj)  (5.133) (5.134)  + {fH-JUj){Ui - £ i - £ , i ) ] } / 3 ,  ^3K-d(&i.i) = WfritiQi-^ + iHrdurQrQrruj  (5-135)  +/4(r,j)(r 4 j • g | • JUj){Ui • Q±• H,j)]/9, where 7 e w is the Ewald constant, qi is the charge of the i  particle, A c is the volume of  the unit cell, the prime on the sums in Eq.(5.129) indicates that i ^ j for R = 0, and Mr)  = - — /„_i(r), r dr  /o(r) = 5.1.5  erfc(7e W |r + i£|)  —Fi—•  D e g e n e r a t e T e r m s of t h e Functional  There are two ways that terms of the functional are degenerate for general crystals (specific crystals may have higher degeneracy). These are through particle exchange, and through elimination of complex quantities. A fundamental symmetry of the direct correlation is under exchange of particles Qj (Zj 2 , flj, Oj) = Cji(r2i,Q.2,£Li),  (5.136)  which can be shown to result in the condition c  nml  /  \  /  \m+n  ^J«(r2i) = ( - )  rnnl /  \  V,ti(ri2)-  Consider the effect on Eq.(5.108) of exchanging the pairs (i,j) (211,212) (0.1,0.2) {KiBl)-  ( r i o7\  (5.137) (p, q) (m, n) (/i, v) (//, v')  We see immediately that at least three factors remain invariant,  Chapter 5. Calculation of Phase  namely, fmnl,  W^  Coexistence  80  ^W^"'", and fip(r.\)fiq{l.2)-  Collecting the remaining factors we see  that  (; ;< ;,)ci(^i)4'o(-ii2) = (-r+n+'(; ; ;, ^-r+^^x-y^w^), (5.138) so that each individual term remains invariant. Therefore we can effect the summations in such a way as to eliminate one term out of every pair representing a particle exchange, and multiply t h e surviving term by 2. The one term which may not be treated in this way is the one given by i = j,p = q,m = n,/j, = v,\i' = v',R = R! (or rather, R = 0, since we have set Rf = 0). Within each term of the summation, the quantities W / i™ M ,c™^(r 1 2 ), and Rl\i0(Li2) may b e complex, as may also be their product. Consider the two terms obtained through simultaneous sign change of /z, u, / / , v\ whose sum is  +(; ; -^w^wf^RL,^)^^)  x(-)A'[M'o(ii2)]*(-)m+n+^+,/ei(r12)]* =  2  ( ; ;  )MWr/"Wf^Rlyo{r12)c^tJ{rl2)\.  i  The symmetry property W™*1  M  = (—Y  +>i  [W^'i  M  ]* comes directly from a similar property  of the Wigner matrix element, R™lfi(£l) [75]. T h e symmetry properties of c™J^(ri 2 ) are given in Ref. [58]. These results give us a second method of eliminating one of a pair of particles. T h e only terms to which this does not apply are those for which H = v = / / = v' — 0. For this case, WJ£00 and Rl00(r12)  are real, and from  C l ( r i 2 ) = (-)m+n+'[c0X(ri2)]*,  (5.139)  Chapter 5. Calculation of Phase  we see that c^^{r\2)  Coexistence  81  is real for m + n + / even, and pure imaginary for m + ra + / odd.  Therefore m + n + / can never be odd for fi,u = 0. (This last statement is a general result for rotational invariant expansions.) Combining these two methods above, a subset of all the terms can be reduced in number by a factor of 4. In practice we take advantage of most, but not all, possible reductions. For the case when R = 0,i = j,p  = q,m  = n the combination of both  kinds of reduction requires considerable care to sort out properly. For this subset, then, we employ only reduction by sign exchange. The extra number of terms that must be calculated because of this is very small in comparison to the total functional.  5.1.6  Effect of P r o t o n D i s o r d e r  It is evident that a phase may contain certain types of order without being fully ordered. Liquid crystals contain orientational order, but lack the full spatial ordering of a solid. Plastic crystals are ordered solids whose freely rotating, non-spherical particles have no orientational order. Such disordered crystals are simply described by the present theory with gip(Q.i) — 1/87T2. Among crystals with orientational order, we must distinguish between those with regular order (the positions and orientations are the same in each cell) and those with random order (the positions are the same in each cell, but the orientations are randomly distributed among only a few fixed directions). A pertinent example of a randomly oriented crystal is proton-disordered ice, which will be discussed further in Section 5.3. The oxygen atoms of ice possess a well-understood spatial ordering which is the same in all cells. The protons are also understood to lie in directions so as to form hydrogen bonds with nearest neighbours, so that the water molecule is not freely rotating. However, each molecule has four nearest neighbours, and the two towards which the protons are oriented differs from cell to cell. This leads to a residual entropy (STe  ) of ice at low temperatures, which has been measured [86].  Chapter 5. Calculation of Phase  Coexistence  82  This random ordering cannot be straightforwardly described by the present theory. The general theory of course has no such restrictions, but here we rely upon a parametrization of a regularly ordered crystal density. The proton-ordered crystal was also the choice for the study of the freezing of water in Ref.  [33]. There are no doubt a number of ways  one might try to account for the extra entropy in these cases. We give one such method which is applied in our calculations. Suppose we know the correct second-order functional for the randomly ordered crystal which we write as f3(nTan-nL)  =  JdrJdCl  j /  a n  (r, fl) In (pmU  (£  ' ~) ]  (5.140)  -(/ a n (£,0)-^//n)} ~ 2 / dl^ J d$2a / dr2 J d 0 2 c ( r 1 2 , f^, 0 2 ) x I ^ C & . Q i ) - / W W a n ( r 2 , H 2 ) - PL/In] , where the superscript "ran" refers to the randomly ordered crystal, and for convenience we have assumed a single component. We can rearrange the ideal term as Jdrjdalpian(r,Q.)\n =  fdr  f d£L{pTAn(r,n)ln(ApTAn(r,Q.))  - J drJdQpTan(r,n) =  ( ^ j } ^ )  f3F™n'id™l-psV\n(ApL/IQ)  In (\pL/Ia)  ~ ( P ^ ^ Q ) - pL/IQ)\  (5.141)  - 1} + VpL  + VpL.  (5.142) (5.143)  where we have used Eqs (3.40), (3.41), and (3.49). Consider the effect on each term of replacing /9 r a n (r, fi) with pOT (r, 0.), where the superscript "ord" refers to the fully ordered crystal. First we would have pran,ideal _ pord,ideal _ mcresid  / r IAA)  Chapter 5.  Calculation of Phase  Coexistence  83  in which ^ o r d > l d e a ' j s evaluated at the average density of the disordered solid. Next, there is obviously no effect on V PL, and —psV hx^KpL / IQ) is also unaffected since minimization cannot distinguish between ordered and disordered solids with the same average density in this term. Finally, the interactive term would be expected to change in going from P™n(r,Q)  to pOTd(r,Q),  but as a simple approximation, we will assume that it does not  change, an assumption which can be tested. The result, given this approximation, is that the functional for the randomly ordered system is obtained by simply subtracting TSTe from the functional for the ordered crystal. This negative contribution should result in freezing at higher temperatures, as the increased entropy of the ice would stabilize it relative to the liquid.  5.1.7  Effect of Solvent Polarizability  We should briefly consider the role of molecular polarization in the density functional theory presented here. As discussed in Chapter 4, polarization is accounted for by the self-consistent mean field (SCMF) method [63] in which each solvent particle is assigned an effective dipole. In a real system this effective dipole is the result of some work done on the particle by the local field to change its internal coordinates, and is different for every particle. In the mean field approximation, every particle is assumed to experience the average field, giving rise to the the same dipole in every particle, the effective dipole of the system. In calculating the internal energy of a polarizable liquid within the SCMF method, one calculates the internal energy of the effective system in the usual way [47], and then adds a correction for the energy required to polarize each molecule. It is important to consider how this approximation affects the thermodynamic quantities related to equilibrium, namely, the equality of temperature, pressure, and chemical potential of the two systems. The temperature is trivially the same in the polarizable and  Chapter 5.  Calculation of Phase  Coexistence  84  effective liquids, and in the solid. The pressures of the polarizable and effective liquids are also equal, since the internal coordinates make no direct contribution to the pressure, and when requirements of density functional theory are satisfied, these are also equal to the pressure of the solid. Now consider the chemical potential, fi, the work required to introduce an isolated particle into the system at constant temperature and pressure. This can be divided into contributions from two separate steps. One is the work required to polarize an isolated molecule with permanent dipole ^ into a molecule with an effective dipole, m. The other is the work of introducing a particle with effective dipole m into the effective liquid. We denote these as the work of polarization, /xpo1, and the work of insertion, fi>ns. We thus have  li = r  l  + »in%  and a similar expression for the solid. Now if the solid were treated separately by an inhomogeneous SCMF method, it would have a different effective dipole, but we are assuming the particles in the solid to be the same as those in the liquid. (This is probably a good approximation for ice, and less so for hydrated crystals.) Therefore fig0 = n^° and since the density functional theory imposes fi1™ = [i1™, we necessarily have ^s — I^L-, and no correction needs to be made at this level of approximation.  5.1.8  A P a r a m e t r i z a t i o n of t h e A n g u l a r D i s t r i b u t i o n  All previous molecular density functional methods based on harmonic expansions [31,32] have dealt with linear molecules. Most of these have used a parametrization of the angular distribution of the form exp J2 anPn(cos 0) where 6 is the deviation of the symmetry axis of the molecule from its average orientation. This is usually truncated at one term (e.g. Ref.  [32]), with the lowest allowed term  Chapter 5. Calculation of Phase  Coexistence  85  corresponding to either n = 1 or n = 2, depending on whether or not the particle has a plane of symmetry perpendicular to the symmetry axis. Another study [31] employed the form exp [ax cosn"(9 - 0,) + a2 cos(</> -  fc)],  where 9 and <f> are now the space-fixed angles describing the orientation of the axis, and #,• and 4>i give the average orientation. Again, np equals 1 or 2, depending on the symmetry. Our parametrization is closer in spirit to the first form, but is generalized to account for rotations about the unique axis. Any rotation of a three-dimensional body may be expressed as a combination of rotations about the three principal moments of inertia. In the body-fixed frame, let these three moments coincide with the x, y, z axes, and let the 2-axis also be the principal axis of symmetry. Then a rotation of the molecule may be described as a rotation about the axis of symmetry, as well as a rotation of the axis of symmetry about the x and y axes. Let us consider these latter rotations first. Gaussian oscillations of the unit vector of the principal axis can be described (for small oscillations) by a distribution of the form exp[-Syr2x  - 8xr2y]  =  exp[—Sy sin 2 $ cos 2 <f> — 8X sin 2 0 sin 2 <f\  =  e x p [ / ( ^ ) ( ^ c o s 2 ^ + ^sin 2 <^)],  (5.145)  where rx = sin 9 cos <f>, ry = sin 8 sin <j> are components of the unit vector expressed in _i  spherical polar coordinates, and where 8i  2  is the Gaussian width of oscillation about the  i-axis. For small 9 we can write f(9) = - sin 2 9 = cos 2 6 - 1 « -92 « 2(cos 9 - 1).  (5.146)  The first two expressions of f(9) have a periodicity of 7r, so that they give distributions  Chapter 5. Calculation of Phase  Coexistence  86  for molecules with a plane of symmetry perpendicular to the principal axis. The last expression for f(0)  has a period of 2TT, and describes the distribution of a molecule  without a plane of symmetry perpendicular to the principal axis. For the case when 6 = 8X = 8y the quantity in Eq.(5.145) simplifies to exp[8f(0)}.  From the above we can  construct a useful parametrization of the oscillation about the z-axis. If ij> is the Euler angle describing rotation about the principal axis, then we require a parametrization which is approximately equal to exp[—aip2], and which has a period of 27r/m z for a particle with an m^-fold axis of symmetry.  One possible choice which satisfies these  requirements is exp[8z(cos mzip — 1)]. We choose for consideration the angular distribution g(<f>,0,i/>) = Azexp[t(cos(mzip)-l)]Axytnpexp[(cosnp  0-l)(8y  cos2 <f>+8xsin2 </))], (5.147)  where A, and AXy,np 3-r^ normalization constants, and np may be 1 or 2. The angles 4>,0,(f) in Eq.(5.147) are angles which give the orientation of a particle relative to its average orientation. We denote them with the superscript "rel". Some relations of these to space-fixed and average angles are ^,rel  cos(0reI)  =  ^ space _ ^  =  s i n 0 s p a c e sin 0ip cos((^ s P a c e — <^ip) + C os 0 s P a c e cos 0, p .  (5<U8)  The angles (4>iPi0ip,ipip) give the average orientation in the space-fixed frame of the p particle of species i. We need to be able to calculate the quantities f?,p and W/t™M M from Eqs (5.102) and (5.109). For practical reasons it is desirable that these integrals involving 9ip(£Li) be reasonably tractable. Therefore, we restrict ourselves at this point to the case 8 = 8X = 8y. Physically this implies assuming that the restoring force to motion of the principal axis away from its average orientation is the same in all directions. The effect of this assumption will be discussed in Section 5.3.  Chapter 5. Calculation of Phase  87  Coexistence  M  Because the integrals in B{p and W^  are over all angles, they are invariant with  respect to rotation, and it will be often convenient in the work below to rotate to the body-fixed frame. This will sometimes involve rotation of the spherical harmonics, given by [75]  RlmMspace)  = E <,-(&p)/?- 0 (a reI )•  (5-149)  The normalization integrals (calculated in the body-fixed frame) are r2w  J 'o = (A%,tl)-1  (i^j)"1  #exp[6P(cosm2^-l)]  27re-**/o(6p), t2ir  =  fit  d<f> J  dO sin  =  4?re-6">i0(8ip),  =  d<f> Jo  Jo =  (5.150)  2TT. hf-e-**'  0exv[8ip(cos6-I)} (5.151)  d9sm0exp[Sip(cos29-l)} cerf (y/s~).  (5.152)  y vip  where cerf (x) = —i erf (ix), the complex error function. The ideal contribution integrals (also calculated in the body-fixed frame) are Bip  —  J + yo  dijj {Af exp[£ p (cos mzi/j - 1 ) ] | In ^jf o  |J4*P  exp[£ tp (cos rnzip - 1 ) ] |  ^ s i n 6 1 { A ^ i B e x p [ M c o s ^ - l ) ] } In { A * i n e x p [ M c o s « / - ! ) ] } 7  i(6 P ) _ A + ^ fh(M  l n ( A ; M * f l ) + 6 P ( T 7 ? 7 ~ 1 ) + ** tp ( r T T T " 1 ) , » P = 1 Jo(tip) J \io(f>ip)  (5-153)  For the integrals W™M M we consider separately the cases /i = 0 and // ^ 0. For /U = 0 we have  Chapter 5.  Calculation of Phase  88  Coexistence  w£"'° = J dnR^0(nsp&ce  )glp(asp&ce)  = / <m£ /$, (a„)/#(arel )9iP(arel) =  EK>MP)A1PJ0  xA%inpjo =  #exp[^(co8m,0-l)] fye-**  P J R™, 0 (a p )27r4 y>n  I £  d0sm0r™(0)exp[6ip(coSn0-l)}  dd sin 0 P m (cos 0) exp[^ p (cos n 0 - 1)]  — nm (ci \]™SzM. „ _ i •> " p io(8ip) e  —  ^n'oKshp)' x /  -*p  / J e - 5 ' P cerf ( ^ ~ )  c?xPm(a;)exp[^,p(a;2 — 1)], nv = 2 .  This last integral can be performed analytically using /  dxxneSx2  =  0,  n odd  (5.155)  (5.156)  Chapter 5.  Calculation of Phase  Coexistence  89  For fi ^ 0, np = 1 we have  w£"'" = |<m^(O s p a c e W(fl s p a c e ) /•27T  _  J£V 2  I  rfJ,  g-*V"/' e ^p(cosm z (V'-V'ip)-l)  Jo  x A * ?1 r  G?0 sin 0 r™lfl(0) exp[Sip cos 0 cos 0 tp - 1]  «/ 0 /-27T  J'  . ,  c?^ e - 1 ^ ^ exp[£,-p sin 0 sin 0,p cos(<^> — <j>ip)]  o  x27re~t'1''t'ipAXyA  /  d0sin0r™ M (0)exp[£cos0cos0 i p -  l]/^(6 8 p sin0sin0 i p )  e-^>/0(6P) x  2e  ,  1  ,_ , r d0 s i n 0 r » ; (^) c M«»-(«-«.>)-i) e -*.>™'«in^P/ , ( ^ sin0 sin0 i p ).  This last procedure is not applicable for np = 2. Thus, for molecules with a plane of symmetry perpendicular to the principal axis, the present method is restricted to linear models. We note that for molecules with an m 2 -fold axis of symmetry, /J, is always a multiple of mz [58]. Therefore the indices of the modified Bessel functions in Eq.(5.157) are integral.  Chapter 5. Calculation of Phase  5.1.9  Coexistence  90  Summary  We now collect together the results of the previous sections. Using only the case np = 1 from Section 5.1.8, which is appropriate for molecules with no plane of symmetry perpendicular to the axis of symmetry, we have  p&n  = =  v  f  JU,  mi/z  ,_  JL\PiL-Z iL spc Zs c £  '  [5/2 + l n ( ^ L 7 r 3 / V ^ , p / o ( 6 P ) e - 5 ' ^ o ( ^ > ) )  ^  l^i p -  p=l  1  ^(s&H-MSS- )**""]} ~\  5Z [(PiL - Pis)(pjL ~ Pis) ~ Pi,Pj,]$%!ij(0) i)j=l ns  PC  mi  mj  ££££  PUii  47T elec/ + (^'^1"^) [T(eip2 + e.g2)]3/2  R i,j=l p = l g = l roo  Jo  dri2r\2  ex  P  (r»-i2r)2" e-  2  C(ijpqRT;r12)  2  4- e  +PpcUelec, C(ijpqRT;  r12) = £  (5.158)  Ct(ijpqRT;  r 1 2 ) exp I  2 i-jp  ^  12  J  r >-jq  CKu>?4;n2) = E c ^ ( ^ ) _  =  jp  L  2 «p  2 1  L  (5.159)  J«  (5.160) M'I/'A  mn 11V  yymn'n  T 12  t/ I  e  'P^/mz\Qip)  1  e~^/o(6p)  2c-*.>to(^)  r-ini>iv-in'<t>iv  e  (5.161)  cos(e - ^ ) - 1 ) e - 6 " ' s i n e s i n e " ' / M , ( ^ p sin 0 sin 9ip), x / " c/0 s i n ^ / 2 ( ^ ) e ^ (  ^m\"ip)  ^-m/i'O  w£"u = j# 0 (a P ) Jo(^ip) 5resid  -lg ^  e  resi(juai  entropy per particle p of species i, ufj ec is the electrostatic pair  potential between particles of species i and ,;', and U vibrationless crystal.  (5.162)  c  is the electrostatic energy of the  Chapter 5. Calculation of Phase  5.1.10  Coexistence  91  C o m p a r i s o n w i t h t h e fc-space E x p r e s s i o n  We briefly outline the derivation of an expression which allows (/3An) 3 = - 2 E  / dr, J dfa / dr2 J dQ.2Cij(r12, Hi, Qa)Pi{Li, £Li)pj(r.2, 0a)>  (5.163)  to be evaluated in terms of reciprocal lattice vectors. In Eq.(5.104) we gave an expansion of the single particle density in Fourier and generalized spherical harmonic terms. One can straightforwardly show that parametrizing the density as in Eqs (5.91) and (5.92) leads to denning the coefficients of Eq.(5.104) as * U = ^TPc  E  W  t ^ ^ V ^ .  (5.164)  Substituting Eqs (5.104) and (5.164) into Eq.(5.163), expanding c;j(r 1 2 ,Q.i,0. 2 )  m  r  °t a ~  tional invariants, and integrating over £Li,Q.2, one readily obtains (fUM)3  =  -\pl  E E  E  x r ' E f :  e-<fcH+^V/V(^P+^)  / d * / dr2 £  B  , A Hr v ^T'' e "'- i - i e "- 2 ' r 2 ^vo(£i2).  We choose 2j as the origin for the integral over r 2 , and substitute £.2  =  (0An)3 = -i /J |EEE e " ( * ?e?p+ * ie3,)/4c, ' (4l " I, ' ,,+4a ^' ) [d^e-^+M*  cft(n2)  (5.165)  Li + £12 to give  (5.166)  fits  Orthogonality in the r_x integral requires k = k2 — ~ki,  an  d  we can  express e - ' - - 1 2 in  angular and radial contributions using the Rayleigh expansion [75], 00  /  c*-*a = X;«"'i/(*r»)(2/ + 1) E /=0  A=-/  ^Ao(i)[^o(ni 2 )]*.  (5.167)  Chapter 5. Calculation of Phase  Coexistence  92  Orthogonality in the f12 integral then results in  0?Aft)3 = - ^ ^ S ^ 6 " * 2 ' ^ ^ ' ^ ' " ' 1 " ^ ' ij  P<1  k  J  mnl  :  y  - — — —  _fc2(,  u'j/'A' *  .M  '  - y/'cr E E E e~"VNp' ~,q"'e"ViJ' """' ij  X  (5J68)  PI  + 2)/4  Xip)  k  Ed(*o/ w E (: ; ; )^7v^r'^U(£)  (5.169)  uf v1 A  mnl  where the square brackets in Eq.(5.168) enclose the definition of the I  Hankel transform  °f c)Ti/"ij(ri2)- The quantity in square brackets in Eq.(5.169) is equal to the Fourier transform, c(k, Qjp,Q.jq), in the limit that <7,p(Q.i) becomes a delta function. It is convenient to extract the k = 0 contribution, and to do this we return to Eq.(5.168) where the identity ji(0) = Sio gives us the restrictions /, A' = 0, m = n and / / = —v'. Then we have (* = 0 contribution) = - \ & E E E "  ij  C*(0) E h ^ ' W ^ W ^ ,  pq mill;  (5.170)  )j.i  where we have used standard identities for the 3-j symbol. We can obtain a further subset of Eq.(5.170) if we note that m = 0 implies fi = v = / / = 0, or (* = m = 0 contribution) = - y / £ E E O ij  pq  0  )  = ~J  E PisP^ZW-  (5-171)  ij  It is convenient to subtract Eq.(5.171) from (/3ACt)3 and add it to (0ACl)2, resulting in newly defined quantities  C¥)'2 = ~\ £ ( * * - P»)(PiL - ^ ) C , ( 0 ) , ij  (5-172)  Chapter 5.  Calculation of Phase  Coexistence  93  (^)l = -kEEEC(»)E(-r'^' , '' , < t]  pq  (5.173)  m^O  - - / & £ £ £ e-*2(4>+^>/4 e'lte^-r,,) »i  P9 i / o  X mnl mnl  „l,jl\> M'l/'A  V.  M  /  If we compare Eq.(5.173) to a previously obtained expression, Eq. 3.21 of Ref. [31], we see that the latter depends on the quantity ( T J + I J ) , instead of on ( r ^ —Tj g ). The expression of Ref. [31] appears to have a sign error (the difference in subscripts is unimportant), which, if incorporated into the computational procedure, would give erroneous results for crystals with non-zero {r,-}. We do not know if the results of Ref. [31] were affected in this way.  Chapter 5.  5.2  Calculation of Phase  Coexistence  94  A n h y d r o u s Alkali H a l i d e s  5.2.1  M e t h o d s and Crystal S t r u c t u r e s  The properties of these electrolyte solutions necessary for the density functional calculations were obtained by solving the RHNC equations as discussed in Chapter 4. All calculations were carried out at 25°C. As shown in Section 4.2, the CHS model gives activity coefficients which are in very poor agreement with experimental results, but these can be greatly improved by adopting the CHSR9 potential. We shall see below that crystallization is similarly sensitive to details of the ion-ion interactions. Because all particles in the crystal are spherical, there is considerable simplification in the form of the density functional, and we can also replace the Ewald sum treatment of long-range interactions with the Madelung constant method discussed in the first part of Section 5.1.4. Letting + , —, and s refer to the cation, anion, and solvent, we define pi = p+L = p-L , ps = psL , Cu  =  c + + (0) + 2c + _(0) + c__(0) ,  cn  = c;+(o) + 2c;_(o) + ci_(o),  Cis = 3?+.(o) + c™ s (o), c  ss  =  Coo°ss(0) .  For all solutions considered calculations are carried out for two crystal lattices. These are the NaCl and CsCl structures which are obvious choices for the model alkali halides considered below. Both are cubic lattices, with a basis of one ion at (0,0,0) (in unit cell coordinates) and the counterion at ( | , | , | ) , so that there are only three possible minimization parameters in this application, namely, />c,e+, and e_. To avoid confusion  Chapter 5. Calculation of Phase  Coexistence  95  with chemical substance names, we will refer to the NaCl structure as an interpenetrating face-centered cubic (ifcc) lattice, and to t h e CsCl structure as an interpenetrating simple cubic (isc) lattice. We have investigated three different methods of evaluating ftA£l/V, and for convenience we give the explicit expressions for each of these. In Method I all calculations are carried out in r-space using the expression ^  =  2pI + ps-pc[5 +pi(pc  + 2\npI  - pi/2)Cn  + 3ln(7re+e.)]  + ps(pc - pi)Cis  -  \p\Css  _rpcM/RNN+1££gEtimM^i«sfiL\ n„  + ~ yoo  —J=L  where Z(Rij)  I  (5,74)  /2  <,.(*)exp(--j)ta«ft,  is the degeneracy of direct lattice vectors of magnitude i? t J . Method II  also treats the Coulombic contribution in terms of the Madelung constant, but the short range contribution is calculated in fc-space by simplifying Eq.(5.173) to 2§2  =  2PI + ps - pc[5 + 2lnPI ~\pcCii  +  + Pi(pc ~ Pi/2)Cn  S\n(we+e-)] + ps(pc  -TpcM/RNN + —-}222 -—g  - pi)Cis  erfc  -  \p2sCss  (5-175)  ( A ')  -\plY,Z{kW++{k)e-^2l^ff__{k)e-^2'' k  +2P(k)cs+_(k)e-h2^+t2-)/4}  ,  where Z{k) is the degeneracy of reciprocal lattice vectors of magnitude k, and P(k) is defined as follows. When Eq.(5.173) is simplified for spherical particles, Rly0(k) no longer  Chapter 5. Calculation of Phase  Coexistence  96  depends on k. The only remaining dependence on the direction of k is in e*-'^ _ x >\ (We have dropped the p and q indices since these crystals have only one particle of each species per unit cell.) These terms may be grouped in pairs of the form (k, —k) and reduced to the real quantity cos(k- (r,- — Tj)). If the unit cell edge is of length L, then, for i ^ j [79], 27T  k = —(h,k,l),  h,k,l integers,  Li  £«'  —  T-j — ^ I j ' 2' 2)  '  cos(k • (n - r,-)) = ( - l ) h + k + 1 .  (5.176)  Since parity is the same for h + k + 1 as for h 2 + k 2 + l 2 , the sign of cos(fc • (T,- — Zj)) is also dependent only on the magnitude of k and from this comes the definition P(k)  =  (_l)h+k+l  is  Method III the k -space approach is followed throughout and /3ASI/V  In  evaluated by ^  =  2PI + Ps - / 0 c [ 5 + 21n / 9 / + 31n( 7 re + e_)] -\(pc  - PifCn  + PS(PC - pi)CIS  ~\P2C E Z(k)[c++(k)e-^>  - \psCss  +  (5.177)  c__{k)e-^2l>  k  +2P(k)c+_{k)e-k2(t\+t2-)'% A discussion of the three methods is given in Section 5.2.4. For the reasons given below only Method I is adequate for the present problem and, unless otherwise indicated, all results reported here were obtained employing Method I including 30 sets of direct lattice vectors. The procedure used in order to evaluate the integrals in Eq.(5.174) depended upon the ratio of {e,} to the grid width upon which c00Jj(r)  was defined.  For large  {e,} numerical integration methods were found to be sufficiently accurate, whereas for the opposite case polynomial fits (to a few points of c00Jj(r)) integration were found to give more accurate results.  followed by analytical  Chapter 5. Calculation of Phase  5.2.2  Coexistence  97  R e s u l t s for C h a r g e d H a r d Sphere M o d e l s  Some results for CHS models are shown in Figure 5.6 and numerical values for coexistence concentrations and crystal properties are given in Table 5.3.  For NaCl no  physical minima were obtained when constrained to the isc lattice. For the ifcc structure physical minima were achievable and the equilibrium condition j3AQ,/V  = 0 is clearly  met (cf. Figure 5.6). Thus NaCl modeled as charged hard spheres crystallizes into the experimental structure. However, it can be seen from Table 5.3 that the coexistence concentration is nearly twice as large as the experimental value. Also, the crystal parameters e + and e_ are about two orders of magnitude smaller than those obtained in computer simulations [82] of the Tosi-Fumi rigid-ion model [83]. We note that underestimation of these fluctuation parameters is typical of density functional calculations [84(a)], though not as severely as appears to be the case here. For KC1 and Csl both isc and ifcc lattices give physical minima, but for the isc lattice the minimum is always deeper than the ifcc result. This means that the isc lattice is the thermodynamically stable state for the ({fii}, V, T) of the solution. For KC1 and Csl it is difficult to ascertain whether Ar2j s c actually goes to zero for the CHS model. The highest concentrations for which A O j s c was determined (see Figure 5.6) are approaching the stability limit of the solution and the solution properties are changing very rapidly in this region. This is related to the divergences introduced in Section 2.2 and which are discussed more fully in Chapter 6. Because of the divergence at the stability limit, it is not unreasonable to assume that f)j s c does vanish before the stability limit is reached but we have not been able to prove this numerically. However, as discussed below (see Figure 5.7) for CHSR9 models the isc structure is clearly shown to coexist with solution. The fact that we can vary continuously from CHSR9 to CHS models by adjusting the value of the parameter A in the Hamiltonian supports our conjecture that coexistence  Chapter 5.  Calculation of Phase  Coexistence  98  Figure 5.6: Representative plots of ( ^ A O / V ) m j n = 0(p\[QUid — Psolid) ^ o r CHS models. As indicated on the figure, the circles and squares indicate results obtained for isc (i.e., CsCl) and ifcc (i.e., NaCl) lattices, respectively.  Chapter 5. Calculation of Phase Coexistence  99  I  NaCl (ifcc)  o  - CO >—3  c3 e  O  i O  1 O •*•  1  1  ID  o  o  CO  CO  1  1  o 8 01 <•-  O  •  - m  u  © o cv  o o CV  o o  o o CO  1  1  1  o o  CO 1  Chapter 5.  Calculation of Phase Coexistence  structure  A*  molarity  100  PC*.  M;1  e.d'1 — "s  KC1  isc*  0  ~9.01  ~.461  ~.0002  ~.0002  KC1  isc*  .74  6.81  .4611  .00038  .00037  KC1  isc*  1  5.15  .4609  .00049  .00047  KC1  ifcc  0  8.97  .3550  .00038  .00037  KC1  ifcc  .74  4.15  .3550  .00038  .00037  KC1  ifcc  1  2.54  .3550  .00028  .00028  KC1 (exp)  ifcc  -  4.17  .3518  -  -  NaCl  ifcc*  0  9.17  .4992  .00021  .00020  NaCl  ifcc*  .5  5.37  .4992  .00021  .00019  NaCl (exp)  ifcc  -  5.42  .4898  .068°  .066°  Csl  isc*  0  ~8.85  ~.2573  ~.0006  ~.0006  Csl  isc*  .7  8.29  .2570  .00086  .00082  Csl  isc*  1.3  2.82  .2569  .00099  .00093  Csl (exp)  isc  -  2.78  .2295  -  -  Csl  ifcc  0  6.17  .1980  .00077  .00075  Csl  ifcc  .7  3.11  .1978  .00089  .00085  Csl  ifcc  1.3  1.01  .1978  .00096  .00091  Table 5.3: Parameters describing crystal/solution equilibria for different models of anhydrous alkali halides. All calculations are at 25°C and the solution densities were taken to be the experimental (exp) values at 1 atm. The value of A* = 0A/ds = 0Aij/(d*jds) -9 determines the strength of the r repulsion. The thermodynamically stable structures are indicated with an asterisk. a From the molecular dynamics simulation of Ref. [82].  Chapter 5.  Calculation of Phase  Coexistence  Figure 5.7: Variation of ( / # A f i / V ) m | n = /?(pij q uid — Psolid) for CHSR9 models. The symbols are as in Figure 5.6.  101  W1  ^  concentration and A*  Chapter 5. Calculation of Phase  Coexistence  102  CO  oin ^ • ^  $  ID  O  H-3  • * ~6  *  Fi  u J? Z c3 c  O  1 O  o  1 T3  o  »  1  1  1  o 8 to <—  ^  n  -a?  on  N  g  o o  o o 1  o oCVJ  o oCO  1  1  Chapter 5. Calculation of Phase  Coexistence  103  occurs very close to the stability limit for the CHS models. For ifcc lattices it can be clearly shown that for all three salts Afijfcc = 0 before the stability limit is reached. However, as emphasized above, for KCl and Csl (modeled as CHS) the ifcc lattice is unstable with respect to the isc case and hence for these models ifcc crystals could only exist as metastable states. Of course, one could adopt a less rigorous point of view and incorporate the correct experimental lattice into the model itself.  The relevant question then becomes whether or not density functional theory  is capable of predicting a transition subject to this constraint. From this perspective, KCl (constrained to the experimental ifcc lattice) does coexist with solution but at a concentration which is much higher than the experimental value (see Table 5.3). To summarize, for CHS models crystallization transitions are likely predicted for KCl and Csl and definitely predicted for NaCl. By combining these results with absolute stability considerations, as discussed in Chapter 6, we also conclude that the concentration range of supersaturation is very narrow for CHS models. For NaCl and Csl, the density functional calculations predict the experimental structures but this is not true for KCl which occurs naturally as an ifcc lattice. For all three salts, however, the structures obtained are in agreement with the radius ratio rules which is not surprising for simple CHS models at 25°C. We remark that the radius ratio rules are frequently violated for alkali halides, the energy difference between isc and ifcc lattices being small enough to be reversed by small contributions. For example, the influence of ion polarizability is often suggested as a possible factor and this is completely neglected in our models.  5.2.3  Sensitivity to Model  We now consider the effect of model by discussing density functional results for the CHSR9 potential introduced in Section 4.2. Results for different values of A* are included in Table 5.3. It can be seen that although the r~ 9 repulsion clearly has a very  Chapter 5. Calculation of Phase  Coexistence  104  large influence upon the concentrations at which crystal/solution equilibria occur, it does not change the predicted crystal structures. That is, the thermodynamically stable structures remain ifcc for NaCl and isc for KC1 and Csl. It is apparent from the table that the crystal/solution coexistence concentrations can be brought as close as desired to the experimental value by adjusting the repulsion parameter A*. While this is not of great importance in itself, it does indicate that crystallization transitions may well be very sensitive to the short range part of the interionic potential. Adjustment of A* does not have as dramatic an effect on crystal densities and Gaussian widths. Although the {e,} are increased slightly towards simulation results, they remain two to three orders of magnitude smaller. We shall discuss the model dependence of Gaussian widths again in Section 5.3. We note that for KC1 (ignoring for the moment that the experimental ifcc structure is obtained here only as a metastable state) and Csl the values of A* needed to obtain the experimental coexistence concentrations are fairly close to the values (i.e., 1 and 0.7, respectively) which essentially reproduce the experimental activity coefficients of the solutions (see Figure 4.2). However, this is not true for NaCl where the value A* « 0.5 which gives the correct crystallization concentration is much smaller than the value (i.e., 6.5) necessary in order to fit the activity coefficients.  Indeed, with A* = 6.5 we did  not find NaCl crystals in equilibrium with solution at any concentration considered (i.e., above 0.025M). This is undoubtedly related to the fact that for model NaCl solutions with A* = 6.5 the dominant structures are solvent-separated ion pairs and there are essentially no contact ion pairs at all (see Table 4.2). This means that for this system the fluid state correlation functions (particularly for oppositely charged pairs) must be extremely poor approximations to those of the anhydrous crystal. Thus it is likely that for this model we are seeing a complete breakdown of the density functional theory at least as implemented by this approximate method. This is an expected result [84(b)] for  Chapter 5.  Calculation of Phase  Coexistence  105  systems where the local structure in the fluid is substantially different from that of the solid, and the present example serves to underline the limits of the density functional approach. Finally, we remark that another significant difference between the CHS and CHSR9 models lies in the predicted concentration range for metastable solutions. For the CHS models the difference between the coexistence concentration determined here, and the absolute stability limit determined in Chapter 6, tends to be very small, but the CHSR9 models can be metastable over much wider ranges.  5.2.4  D i r e c t S p a c e v s . Fourier S p a c e M e t h o d s  Let us return briefly to the question of r-space vs. fc-space calculations. From Eqs (5.92) and (5.164) it is clear that the r-space sums will converge with fewer terms as the tip decrease, and that the &-space sums will require fewer terms as the e,p increase. It has been generally concluded that (at least for freezing of hard spheres) it is more efficient and accurate to perform the calculations in r-space [84(a)]. This is also true for the present crystallization problem as illustrated in Figures 5.8 and 5.9 using some results for a 4M Csl solution with A* = 0.8.  In Figure 5.8 the minimum flAQ/V values for isc  crystal structures obtained using Methods I, II and III (defined in Section 5.2.1 and in the figure caption) are shown as functions of the number of sets of lattice vectors included in the calculation. It is obvious that the rate of convergence varies greatly for the three methods. The completely r-space approach (Method I) converges rapidly requiring only about ten sets of direct lattice vectors in order to yield a converged result. Methods II and III which require &-space sums do not converge even with 1300 sets of reciprocal lattice vectors [the upper limit imposed by our numerical evaluation of the c,j(fc)] and clearly are not reliable approaches for the present calculations. Further insight into the problem can be obtained by evaluating (with no minimization) /3AO/V as a function of  Chapter 5.  Calculation of Phase  Coexistence  106  Figure 5.8: The dependence of ( / ? A f i / V ) m m upon the number of sets of lattice vectors employed for different methods of calculation. The results shown are for 4M Csl (isc structure) modeled with the CHSR9 potential (A* = 0.8). In Method I the functional is evaluated using the Eq.(5.174), Method II uses the Eq.(5.175), and Method III uses Eq.(5.177). The abbreviations d.l.v. and r.l.v. refer to direct and reciprocal lattice vectors.  8 •8 c-f.  £ o c  Method I  Method II  -5-W  1 C\) —  a,  4-  o  0  o 5.5-  o° ° ° oo ° o  -130-  §  Method III  o  3.5-  0  8.  o  -140-  oo  3-  0  o  0  3  o  0  o  6.5 0  10  20  sets of d.Lv.'s  30  1 Du -\ ()  o ° ° 1  1  500  1000  sets of r.l.v.'s  1500  c.o  0  i  i  500  1000  1500  sets of r.l.v.'s  o  Chapter 5. Calculation of Phase  Coexistence  108  Figure 5.9: Evaluation of functional ((3Att/V) (no minimization) for selected solid state parameters. As indicated on the figure, the curves represent the three different methods of calculation defined in Figure 5.8. The results are for the same Csl system as in Figure 5.8.  Chapter 5. Calculation of Phase  i;  C\2  £; lO cv•  / /  - -i=:  K ,\  w  Tlo Q.  < c < ox  I  Qo:  ©  -?  1  i  o  o  O  lO 1  1  V  c3  1  ID  1  ^  u;  •«.  1 1 ,  **>  octf  s° o  /  / / /  c<  - w ^_^ 1 o  /  / >  bJD  ^o  =i=i  II «  ^ UJ  /  •  T3  /  o  1 !  TT  1  /  m  CO  ^^^  i i  o  iD " C\2  -7'  ; ) .'  II II CO  109  Coexistence  /  o  i—H _ ^H  1  l  o1  o  1  i: / •  -7 jj  II CO  w  1 < 0  i O  i ID  CO  C  1  Chapter 5. Calculation of Phase  Coexistence  110  e = e + = e_ for selected solid state parameters. The results obtained using the same Csl solution input are shown in Figure 5.9. We see from the figure that all three methods agree at large e (i.e., edj1 > 0.1) but part company and give very different results at the smaller values. Since the r-space method becomes more accurate as e decreases, it is clearly the only reliable method for the ranges of e, (see Table 5.3) relevant for the present models and theory.  Chapter 5. Calculation of Phase  5.3  Coexistence  111  Ice Ih  5.3.1  Crystal S t r u c t u r e  Ice occurs in a number of different forms [85]. The principal feature which all forms share in common is the hydrogen bonding of each molecule to four others in an approximately tetrahedral arrangement. They differ in the degree of distortion from this tetrahedrality. Some structures are proton-ordered, in which proton positions are identical from cell to cell, while most are proton-disordered. An important trend among the ices is the distortion of the water molecule with increasing pressure, which accompanies distortion of the hydrogen bonding. For this reason we restrict our choice of possible structures to ice Ih, which is the common hexagonal crystal found at atmospheric pressure. A preliminary discussion of ice Ih was given in Section 5.1.6, and we discuss its structure [86] in more detail here. Ice Ih has a hexagonal crystal structure described by the space group P63/mmc [87]. A convenient choice of primitive cell is shown in Figure 5.10.  Figure 5.10(a) shows  the geometry of the cell, and Figure 5.10(b) shows the placement of the tetramolecular basis within the cell. The molecules are shown as being tetrahedral in shape. One of the molecules lies at a vertex, with the remaining vertices occupied by replications of this basis molecule. Four vertices, labeled A, B, C, and D form the corners of a rectangle which passes through each molecule in the basis. This rectangle is shown in Figures 5.10(c) and 5.10(d), with two different sets of H 2 0 orientations. Protons are represented by + and lone pairs by —, but only those in the plane of the page are shown. These orientations, and others, are found randomly in cells of naturally occurring ice, which is orientationally disordered. To construct proton-ordered ice, though, there are far fewer arrangements possible than in proton-disordered ice, since each pair of nearest neighbours must have exactly one proton between them. The diagrams in Figures 5.10(c) and 5.10(d) are two  Chapter 5. Calculation of Phase  Coexistence  112  Figure 5.10: Crystal structure of ice Ih. (a) The primitive hexagonal cell, (b) Locations of four tetrahedral basis molecules. The two molecules inside the cell are on the line joining the centres of two equilateral triangles shown. (c),(d) Two different sets of orientations which can be used to generate proton-ordered ice. Protons are represented by + , and lone pairs by —. (e) Illustration of the sixfold symmetry collectively possessed by the unit cells.  Chapter 5. Calculation of Phase Coexistence  113  (a) ^  ^  \  ^  c=7.367A  .. •'I 60 "  r >^  b=a  a=4.523A  (d)  (e)  Chapter 5. Calculation of Phase  Coexistence  114  possible arrangements which may be replicated to give proton-ordered ice. Finally, in Figure 5.10(e) we indicate in a view along the c-axis how replications of the primitive cell result in a hexagonal crystal, even though the primitive cell itself does not possess such symmetry. A good approximation to the experimental structure of ice is to assume that c/a = Wo/3 and d/a = l / v 2 4 , which leads to all angles and orientations being exactly tetrahedral. These values for c/a and d/a differ from experiment by less than 0.2%. Any time we impose some restriction such as this on the parametrization of the solid, we constrain the minimization to give results which are an upper bound to minima obtained by less restricted parametrizations. We would expect the difference in this case to be small. As indicated above, we have out of necessity chosen a proton-ordered structure, as have other workers [33]. This is in contradiction to the experimental structure, and we wish to approximate the correction for this using the method of Section 5.1.6. To do this, we require the residual entropy of ice Ih. This can be obtained from existing measurements [86], but a well known calculation by Pauling [88] gives a convenient value which is within the limits of experimental uncertainty. In this calculation, each hydrogen atom is assumed to lie on a line joining two oxygen atoms, being bonded to one and coordinated to the other. This permits the H atom two possible locations, depending on which oxygen it is closest to. Therefore in an A^-molecule ice crystal there are possible configurations. These are not all permissible, however, as some would correspond to species such as O H - , H 4 0 + 2 , etc. For each 0 atom, there are sixteen configurations of the H atoms surrounding it, and six of these configurations correspond to the chemical compound H 2 0 . We can use this to reduce the number of configurations to give 5""esid = k In W = k \n[22N(6/W)N] (3(Tsiesid)  = S/(kV)  = Nk ln(3/2),  = />H2o ln(3/2).  (5.178) (5.179)  Chapter 5.  5.3.2  Calculation of Phase  Coexistence  115  Results and Discussion  We employ two different models to study the freezing of water.  The first is the  same model for water discussed in Chapter 4, which is a polarizable multipolar hard sphere (MHS). The second model differs by having an additional r~ 9 spherical repulsion, identical to the repulsion added to ions in the CHSR9 model. This second model of water we refer to as MHSR9. We will study the temperature and density dependence of the minimized functional for the MHS model, and then study model dependence at a single temperature and density using the MHSR9 model. All of these minimizations will be for a proton-ordered crystal with proton positions given by Figure 5.10(c). We also employ the approximation discussed in the previous section where all the distances and angles are tetrahedral. From these minimizations we will be able to explore reasons for the model dependence of both coexistence and solid structure. We briefly consider the role of residual entropy, and compare our density functional results for freezing of water with those of Ref. [33]. Water freezes at 1 atmosphere at 0°C and a density of 1 g/cm 3 , or a reduced density of p* = .7341 in particle diameter units. Agreement with experiment would be obtained for a given model if freezing were predicted at this state point. In Figure 5.11(a) we plot the value of the minimized functional against the temperature of a liquid at p* = .7341.  The  temperature dependence does not actually extend to 0°C, as we are unable to calculate liquid correlations below about 12°C. However, coexistence is predicted, in this case at a temperature about three times larger than the experimental absolute temperature. We discuss below the reason for this extreme discrepancy.  It should be emphasized that  the liquids corresponding to the different temperatures do not all represent real water at 1 atm. The density at all points is the experimental density of water at 0°C and 1 atm. Different results would no doubt be obtained using liquid correlations with the  Chapter 5.  Calculation of Phase  Coexistence  116  Figure 5.11: Minimized functional values for proton-ordered ice Ih. (a) Variation with temperature for the experimental density at 0°C and 1 atm. (b) Variation with density at 20°C. (c) Variation with parametrized r - 9 repulsion on the hard core. The temperature is 25°C and the density is the experimental density at 25° C and 1 atm.  117  Chapter 5. Calculation of Phase Coexistence  p* = .7341 lOO-i  ^__^ e _  0-  (a)  —-€)  -100-200-OUU T  0  (b)  (c)  |  200  |  T(°C)  400  600  Chapter 5.  Calculation of Phase  Coexistence  118  experimental density and effective dipole appropriate to each temperature. However, this would not change the essential nature of the results, and would not affect the direction of the temperature dependence of the functional. Because we are unable to obtain liquid correlations at low temperatures, the density dependence in Figure 5.11(b) is given, not at 0°C, but at 20°C. The minimized functional is not equal to zero between either of the limits in this plot, nor can we calculate any liquid correlations outside of these limits. In connection with Figure 5.11(a), this implies that, at least for the present density functional theory, the MHS model is everywhere metastable with respect to solidification at this temperature. In order to give some explanation of the poor agreement with experiment obtained for the MHS model, we plot in Figure 5.11(c) the value of the minimized functional for the MHSR9 model for a number of different values of the repulsion parameter, A*. The temperature is 25°C and the density is the experimental density at 25°C and 1 atm. We see in these results a strong sensitivity to model. We can explain this effect from at least two perspectives. First, if we recall that the quantity displayed in Figure 5.11 is /?(/)ijqujd ~~ Psoiid), it is apparent that for the MHS model density functional theory predicts solids with much higher pressures than the liquids. This is not surprising, since the virial pressure of the hard core liquid at 25° C is negative at the experimental density, while it is entirely possible that the solid selects a more reasonable density corresponding to a positive, or at least less negative, pressure. As repulsion is added, the pressure of the liquid will naturally increase as well, eventually rectifying the difficulty. As one would expect, pressure also increases with temperature in Figure 5.11(a). This suggests that the MHS liquid at p* = .7341 would actually freeze at a very high temperature, and somewhat different models, such as MHSR9 with a large repulsion, are required to obtain freezing under conditions closer to those of real water. The second perspective from which we can explain the trends in Figures 5.11(a) and 5.11(c) is from investigation  Chapter 5. Calculation  of Phase  Coexistence  119  of the spherically averaged direct correlation functions. As we can see in Figure 5.12, the change effected by temperature or short-range repulsion is dramatic in the region of contact.  The Gaussian peak in this region generally makes an important contribution  to the functional, and lowering the contact value results in making the functional less negative. By dissection of the functional, we can show that this quantity and similar ones for non-spherical projections are causing the large changes. This is then a distinct weakness in using MHS models of real systems in crystallization studies. Because of the cusp-shaped contact, the direct correlation in this region will be artificially large, and thus one would generally expect the functional to be unduly negative. Any model which corrects this anomaly, even such an ad hoc technique as used here, is almost certain to improve these values. We are not asserting the superiority of these repulsive systems as models of the liquid, since we have not performed a thorough investigation of their ability to predict liquid properties, but the lesson from this is clearly that freezing and crystallization transitions are very sensitive to model, and that models involving hard cores appear to be a poor choice for modelling transitions of real systems. These effects are related as well to the extremely small Gaussian widths predicted for CHS models in Section 5.2, and for MHS models as shown in Table 5.4. Since a Gaussian peak selects values of the correlation function in such a way as to minimize the functional, when the largest value at contact occurs at a cusp, rather than a rounded peak, small Gaussians and near-contact of nearest neighbours are very desirable features of the optimal density. The average densities are not shown, but we note that they are correlated to Gaussian widths in the sense that the distance between the surfaces of nearest neighbours in their average positions is always of the same order as the Gaussian width. It is sensible to conclude that a model whose correlation functions have rounded peaks at contact would be predicted to have wider Gaussian peaks, and a less close-packed  density.  Chapter 5.  Calculation of Phase  Coexistence  120  Figure 5.12: Spherically averaged direct correlation functions of water near contact. (a)Three different values of the repulsion parameter (cf. Figure 5.11(c)). (b)Three different temperatures (cf. Figure 5.11(a)).  Chapter 5. Calculation of Phase  (b)  p*  121  Coexistence  = .7341 25°C  200°C  500°C  1e 13  10-  5\\ \  0-  1  \\  \\ \\ \  ^^•"^^I'^C^rs-—_ 1 "  1.1  i  1.2  v rlA  1.3  -  i  1.4  1. 5  Chapter 5. Calculation of Phase  Coexistence  122  We consider briefly another feature of the predicted density. The density of ice at 0°C is 91.7% of the coexisting liquid density, which corresponds to a nearest-neighbour 0  o  distance of ~ 2.76A. Because we have used 2.8A as the diameter of H2O, even at close packing the density is only 88% of the liquid density at 0°C, so that complete agreement with experiment is not possible for physical densities. The fact that the density decreases upon freezing is one of the special features of H2O, but neither this density functional theory nor the molecular theory in Ref. [33] can be truly credited with predicting it, since it is a direct result of the input crystal symmetry. The final aspect of the parametrized density that we can compare to experiment are the angular oscillations. For sufficiently small oscillations we have 62 « sin 2 0 = 2(1 - c o s f l ) , which, from Table 5.4, we can see is true for all values considered. We can therefore obtain the mean square deviations from integrals similar to those in Eq.(5.153). The values of the root-mean-square (r.m.s.) deviation we obtain are in the ranges indicated in Table 5.4. The rock and wag modes are both oscillations of 0, one about the z-axis and the other about the y-axis, and their width is related to S. In general we would not expect the rock and wag to be the same (6X / Sy). However, because of the tetrahedral quadrupole we employ, and the tetrahedral coordination in ice Ih, the two modes should be nearly identical in our model. The twist is an oscillation of V> about the symmetry axis, and its width is related to £. The principal result of interest to us at present is to compare the magnitude of values obtained by the theory to typical experimental values, and it is clear from the table that for crystals near coexistence, the r.m.s. deviations predicted are of the appropriate order of magnitude, though they may be off by an order of magnitude away from coexistence. This is in contrast to the Gaussian widths, which are still an order of magnitude off at coexistence, and up to four orders off elsewhere. This is consistent with  Chapter 5.  Calculation of Phase  Coexistence  123  500° C p* = .65 p* = .95  A=0  A = 16  900cm -1  300cm"1  .000044  .0072  .000046  .000065  .000055  .015  -  -  Rock  .027  .17  .028  .034  .030  .10  .11  .25  Wag  .027  .17  .028  .034  .030  .10  .18  .39  Twist  .020  .30  .020  .023  .021  .065  .13  .28  Mode  12°C  Transl.  Table 5.4: Root-mean-sqare (r.m.s.) spatial and angular distribution widths for water molecules. The first six columns refer to theoretical results in ice Ih, and have reference to Figure 5.11. The last two columns are typical values at 25°C for librational modes of water in various crystals, when the modes are located at the wavelengths shown (adapted from Table I of Ref. [91]). For the translational mode the values are in units of the particle diameter. For the rock, wag, and twist librational modes all values are in units of TT. For purposes of comparison, the librational modes of water in ice are in the region around 650 c m " 1 (Ref. [92]). our explanation of the sharpness of the positional Gaussians, since the cusp at contact would not be expected to have the same effect on the angular distributions as on the spatial distributions. We conclude that the predicted angular density is not as sensitive to differences between MHS and MHSR9 models as is the spatial density. We briefly consider now the role of proton disorder in the crystal, using the MHS model. In Table 5.5 we list the minimized value of the functional at 25°C for the two proton orderings in Figures 5.10(c) and 5.10(d), both with and without the residual entropy term suggested from Section 5.1.6.  The value of the residual entropy term,  pc ln(3/2), is ~ .264, consistent with the results in Table 5.5. It is clear that the effect of changing the structure is about an order of magnitude greater than the effect of the added residual entropy term. In Section 5.1.6 we assumed that the difference in contribution from different structures was negligible, and since this assumption is invalid we conclude that this method of accounting for proton disorder is inappropriate, at least for ice.  Chapters.  Calculation of Phase Coexistence  Structure Fig. 5.10(c) Fig. 5.10(c) Fig. 5.10(d) Fig. 5.10(d)  Residual Entropy No Yes No Yes  124  Functional -185.93 -186.19 -183.36 -183.62  Table 5.5: Contributions to proton disorder. We note that, as required, the residual entropy contribution did lower the value of the functional which, in conjunction with Figure 5.11(a), implies a raising of the temperature of fusion. We conclude this section with a comparison to the theory for the freezing of water found in Ref. [33]. The similarities between our methods are that a second-order theory is employed, and that the crystal is proton-ordered, though the specific nature of their ordering is not indicated. Some differences are that in Ref.  [33] a density functional  theory based on site-site correlations is used, and evaluation is performed in &-space. Furthermore, Ref. [33] is not a Hamiltonian calculation from first principles, but rather, the correlation functions are obtained from inversion of neutron diffraction measurements of real water at 4°C and 25°C. The results obtained in Ref. [33] are more realistic than those obtained by our method for the MHS model. Freezing is predicted at —5°C, and the Gaussian widths are ~ .05 for oxygen atoms and ~ .25 for protons. While the many differences between the two methods make comparison difficult, from our study of the effect of different potentials we assume that the use of experimental rather than model correlation functions could by itself account for the difference. This emphasizes the need for better models when applying density functional theory to the study of real systems.  Chapter 5.  5.4  Calculation of Phase  Coexistence  125  Alkali H a l i d e H y d r a t e s The alkali halides which form crystal hydrates invariably contain a small ion (Li + ,  N a + , or F~) and a larger counterion. The smaller ion interacts strongly with the water molecules, providing the drop in energy necessary to compensate the decreased entropy of the oriented solvent. Thus, for a compound which possesses various hydrate structures, those with a lesser proportion of water are stable at higher temperatures, where entropy has a greater weight, as illustrated in Figure 1.1.  5.4.1  Crystal S t r u c t u r e of N a C l 2 H 2 0  The only stable hydrate of NaCl at 1 a t m is NaCl • 2H2O which possesses a monoclinic structure [93]. The primitive monoclinic cell has three unequal edge lengths, a, b, and c. The unique axis, c, is at right angles to the remaining two, which themselves intersect at an arbitrary angle, 7. The point group [87] of the crystal is 2/m, meaning that, in addition to the identity, the unit cell also possesses a two-fold axis (along c), a mirror plane (perpendicular to c), and the combination of these, an inversion. The specific space group [87] is P2x/6, meaning that the cell is primitive, the two-fold rotation is a screw axis which advances half a unit along c with each rotation, and the mirror is a glide plane which advances half a unit along b with each reflection. The combination of the screw axis and the glide plane gives an inversion about some point other than the origin. We represent these symmetry elements as identity: screw axis: glide plane: inversion:  (x,y,z), (—x, —y, z + j ) , (x,y+^,-z), (~x, j — y, - — z).  Chapter 5. Calculation of Phase  Coexistence  126  The center of inversion is about the point (0, j , j ) . It is conventional with this space group to move the center of inversion to the origin, so that we shall actually be using the expressions identity: screw axis: glide plane: inversion: If a particle is located at (x,y,z),  (x,y,z), (—x, \ — y, z + | ) , {x,y-\-\,\-  z),  (—x,—y,—z).  then three other symmetry related particles will be  located at the coordinates above. The orientations of symmetry related particles depend only upon the point group elements, and are given as  identity:  (0,(f>,tp),  rotation:  (0, <f> + ir, ij)),  reflection:  (ir — 0, <f>, —tp),  inversion:  (TT — 0, <j> + T , — •*/>),  where 0,<j>,i/> are the Euler angles [75]. We are concerned here only with non-chiral molecules. Otherwise the mirror and inversion give the enantiomer. In Table 5.6 we give the experimental coordinates for four of the particles in the unit cell of N a C l - 2 H 2 0 , namely, Na+, CI", H 2 0 ( 1 ) , and H 2 0 ( 2 ) . The coordinates of the other twelve are generated by the symmetry elements at each stage of the minimization, as the minimization parameters of symmetry related points are necessarily identical. From consideration of nearest-neighbour distances one can see that each water molecule is surrounded by two cations and two anions in a roughly tetrahedral configuration, and that each ion is surrounded by four water molecules and two counterions, in a roughly  Chapter 5. Calculation of Phase Coexistence  127  <t>  e  a  b  c  Na  .02374  .45503  .17020  —  —  —  CI  .29210  .21354  .12190  —  —  —  H 2 0(1)  .776  .169  .324  1.127T  .324TT  -.848TT  H20(2)  .238  .241  .492  .294TT  .4977T  -.648TT  a  Table 5.6: Experimental coordinates for NaCl-2H20 . Values for a, b, and c are in unit cell coordinates. octahedral configuration. Unlike ice Ih, the proton positions are fully ordered [94], so that no residual entropy would be predicted for these hydrates.  5.4.2  Crystal S t r u c t u r e of L i l 3 H 2 0 LiI-3H20 has a hexagonal unit cell [95], similar to that of ice Ih in Figure 5.10(a), o  o  but with the experimental cell dimensions a = 7.45A, c = 5.45A. In Figure 5.13(a) we give a view along the c axis showing the placement of the two chemical units within the cell. Vertical positions are indicated in units of c. In Figure 5.13(b) we show a different cross-section, in which the vertical axis is c and the horizontal axis is the long diagonal of Figure 5.13(a). From the distances indicated here, the only position parameter which is variable is the distance r, where the water molecule shown is located at (r, r, | ) in unit cell coordinates. Given a set of particle diameters, and assuming close-packing, we can use the geometric relations implied in Figure 5.13(b) to give the relations among a, c, and r. Choosing a as the independent parameter we can have  Chapter 5. Calculation of Phase Coexistence  128  Figure 5.13: Crystal structure of LiI-3H 2 0 (a) Projection in the a, b plane of the positions of atoms. The circles are, in order of increasing size, L i + , H 2 0 , and I - . Numbers give the vertical position in c-axis coordinates, (b) Geometry of three particles at close packing. The vertical axis is c, and the horizontal axis is the long diagonal of Figure (a), (c) Demonstration of how replication of the primitive cell gives six-fold symmetry.  Chapter 5. Calculation of Phase Coexistence  129  (a)  (b) v/3a  (c)  129  Chapter 5.  Calculation of Phase  Coexistence  where we recall that dij = (d,- + dj)/2.  130  These relations are useful in helping us find the  minimized single particle density. The primitive cell of the hexagonal space group does not possess all the symmetry elements, but from Figure 5.13(c) it is clear how replication of the cell generates the sixfold symmetry. The protons are in the plane of the figure, and are oriented towards two anions at the same vertical level, while the lone pairs are oriented towards two cations situated above and below the page on the central axis. Configurations of lower symmetry are possible, especially at low temperature [96], but preliminary calculations have shown that they are not predicted by this model at 25°C.  5.4.3  R e s u l t s and D i s c u s s i o n  To apply the density functional theory to NaCl-2H20 and LH-3H20 , we require liquid correlation functions for aqueous solutions of NaCl and Lil. All RHNC solutions are obtained at 25°C. For NaCl, CHS model solutions are obtained at a number of concentrations up to stability (these are the same as for the anhydrous case), and we also make use of CHSR9 results at 2.0M for various values of A*. Lil solutions present a greater challenge, and we were unable to obtain any CHS or CHSR9 correlations at all. This appeared to be due to the extremely strong cation-anion and cation-solvent interactions. In any case, we were able to obtain RHNC solutions by adding repulsive terms to both the ion-ion and the cation-solvent interactions. This ad hoc adjustment was exploratory in nature, and we did not make any effort to fine tune this model to match experimental properties. Three sets of results are shown in Figure 5.14.  In each plot, one value is obtained  by minimization of the functional, and all other points are obtained by evaluation of the functional for parameter values obtained at the minimized point. The justification for this procedure is straightforward.  In treating molecular systems, one may readily  Chapter 5.  Calculation of Phase  Coexistence  131  Figure 5.14: Functional values for alkali halide hydrates at 25°C. (a) CHS models of NaCl at various concentrations, (b) CHSR9 models of 2.0M NaCl for various values of the repulsion parameter, A*, (c) Modified CHSR9 model of Lil for various concentrations. The parameters values for which the functional is evaluated are discussed in the text.  Chapter 5. Calculation of Phase  Coexistence  132  NaCl-2H 2 0, A*=0 -35  (a)  -40-  -45-  -50 0  2  n  r  4  6  8  10  mol/L  NaCl-2H 2 0, 2.0M "KF.  -38^  1  (b) -40—* >  Q  A")  c»  1  2  1  4  1 A  .e  i  8  1)  Chapter 5.  Calculation of Phase  Coexistence  133  encounter situations in which the number of minimization variables becomes unwieldy for typical optimization routines. It is therefore reasonable, that, having obtained minimizing parameters for one point, these should then be applied to the other points of interest, generating an upper bound to the true minima of these other points. If the values thus obtained are positive, one would then proceed to minimize them more carefully. However, if all values are negative, as is the case of Figure 5.14, there is no interest in obtaining the true minima.  It is clear that, for these models, we will not obtain any crystal  hydrates in equilibrium with solution at this temperature. There are two motivations for not exploring other repulsion parameters for Lil other than the one used to obtain the results shown in Figure 5.14(c). First, by considering Figures 5.11(c) and 5.14(b), it appears that increasing the parameter beyond its present value of A* = 6.5 will not likely result in significant change, while liquid correlations for much smaller parameters are difficult to obtain numerically, as discussed above. Secondly, the evidence from Figure 5.14(b) suggests that we will not obtain as dramatically different a value by varying the repulsion as we did with anhydrous alkali halides and with ice. We can understand at least some of the reason for the relative insensitivity to model exhibited by  NaCl-2H20 . In the anhydrous crystals of Section 5.2, the only nearest  neighbours were cation-anion pairs. Therefore, when increased repulsion lowered the contact value of the cation-anion direct correlation function, the functional was made less negative. In the case of NaCl, this apparently made the crystal become unstable for some value of A* less than 6.5. Hydrate crystals, on the other hand, have both cationanion neighbour pairs, as well as ion-solvent neighbour pairs, both of which give negative contributions to the functional for the CHS model. When repulsion is added, the cationanion contribution again becomes less negative, but the contact value of the ion-solvent correlation functions is enhanced, making its contribution more negative. Thus there is at least a partial cancellation of the effect of adding repulsion in the CHSR9 model,  Chapter 5.  Calculation of Phase  Coexistence  134  with the ion-solvent contributions apparently being more important, according to Figure 5.14(b). We have indicated that most of the points shown in Figure 5.14 are upper bounds to the true minimized functional values. This is true of all the points for  NaCl-2H 2 0  in another sense as well. We have not performed these minimizations with respect to all variables, but have set a number of the variables equal to their experimental values. These variables are the positions and orientations of the particles, the ratios of cell edge lengths, and the angle between the a and b axes. These values were given in Table 5.6. This leaves as properly minimized parameters the average density, the Gaussian widths of translational oscillation, and the widths of orientational oscillation, a reduction from thirty variables to nine. We attempted some minimizations in which other quantities were allowed to vary, and found that substantially lower values were obtained, as low as /3AQ./V  = —73. Even so, the true minima may be considerably lower yet. In Ta-  ble 5.7, we give the surface-to-surface distance of nearest neighbours for the minimized experimental structure, as well as for the more fully minimized crystal. In both cases,  Neighbour Pair  Experimental  Partially Minimized  Na-Cl  .061,.060  .0067, .056  Na-H 2 0(l)  .0054,.026  .00025,.00028  Na-H 2 0(2)  .00019,.0026  .0025, .00053  C1-H20(1)  .13, .18  .30 .16  Cl-H 2 0(2)  .12, .12  .29, .12  Table 5.7: Nearest neighbour surface-to-surface distances in NaCl-2H 2 0 . Distances are in solvent diameters. the minimized Gaussian widths for a given particle are of the same order as the smallest  Chapter 5. Calculation of Phase  Coexistence  135  surface-to-surface distance with any neighbour of that particle. Thus, for the experimental structure, Na and H 2 0 ( l ) have similar Gaussians to what we have seen in other systems, while H20(2) has a somewhat larger value, and CI has a considerably larger value. In the minimized crystal, the Gaussians are somewhat closer to what we expect, and the minimization seems essentially to have moved the water molecules closer to the cations. This is a logical change, but may simply represent a local minimum. For instance, the true minimum for this model may be one in which CI is packed more tightly. At these low temperatures, a serious search for the true minimum, such as would be required if the values in Figures 5.14(a) and (b) were positive, would benefit from being preceded by a search for the minimum energy configuration of the crystal at absolute zero. This would be a much simpler and quicker procedure than the finite temperature minimization of the grand canonical potential we have conducted here.  Chapter 5.  5.5  Calculation of Phase  Coexistence  136  P r a c t i c a l I m p l e m e n t a t i o n of t h e D e n s i t y Functional P r o g r a m  5.5.1  Programming Tests  A number of tests were made during implementation of the density functional theory in order to insure its accuracy. T h e two areas most susceptible to error are the sums over various contributions in the interactive term (Sections 5.1.3) and the Ewald sums (Section 5.1.4). The spherical terms were tested by comparison with results obtained by Method I in Section 5.2.1 for hard sphere freezing and anhydrous alkali halide crystallization. The dipolar contributions from the Ewald sum were tested using configurations and results in a paper by Onsager [97]. Higher-order multipolar contributions were also tested, as we explain below. Essential to the higher-order contributions is the quantity W^*1 M, which was tested against the requirement that it approach -ft™/M(HiP) in the limit 6, ( —> oo. The rotation matrices and Euler angle definitions are a potential source of error in the Ewald sums, and these were tested by allowing the angles of NaCl-2H20 to minimize freely among all possible angles, while other variables were held fixed, to see that they returned to nearly experimental values. Another valuable procedure is to test r-space and Ewald sums simultaneously against each other. In the r-space routine we calculate electrostatic contributions for a number of lattice vectors using harmonic expansion expressions. The Ewald sum, when supplied with a very small Ewald constant, -few, can be made to calculate the same quantity, but by a Cartesian coordinate method. Coulson and Eisenberg [98] performed calculations on electrostatic interactions in ice which are useful to us in testing our calculation of electrostatic contributions for ice, including quadrupolar terms. In Ref. [98], the various contributions to the electrostatic interactions of nearest neighbours are calculated, using an average over the various possible arrangements. To reproduce these numbers, we also calculated nearest-neighbour interactions in our model. From Figure 5.10(b) we can see that each molecule has one  Chapter 5.  Calculation  of Phase  Coexistence  137  neighbour in an eclipsed configuration and the remaining three are staggered. Of each of these there are two asymmetric configurations and one symmetric configuration. Therefore, in a proton-disordered crystal where all configurations have equal probability, we would expect the symmetric and asymmetric versions of staggered and eclipsed configurations to occur with the frequency indicated in Table 5.8. In this table we have also given  Frequency Eclipsed, symmetric  i i _ i 4*3 12  Eclipsed, asymmetric  1.1 — 1  Staggered, symmetric Staggered, asymmetric  4 ' 3 6 3 1 _ 1 4*3 4 3 2 _ 1 4 ' 3 2  U  dd 0 3^  ~f -¥2  u  uqq  dq  ~ ^ Q T -^SHQT  -3><?T -J75"QT  _  56^)2 27VT 9 VT  80^2 ~27L*T  -±£f) 2 9 VT  Table 5.8: Nearest neighbour pair configurations and energies in ice Ih.  the contributions to the electrostatic energy for each such pair in a close-packed crystal. Combining these contributions with the appropriate weights, we compare our result with those adapted from Ref. [98], as shown in Table 5.9. We also compare the magnitude of the total dipole moment (i.e., permanent plus induced) and of the quadrupole moment employed, which allow us to conclude that the two methods are reasonably consistent. Since the quantities in Table 5.8 have been checked by calculation using the Ewald sum and r-space routines, this is a confirmation of the integrity of these two components. A less quantitative but good qualitative test is to consider some of the contributions to the electrostatic energies of the crystal hydrates, as shown in Table 5.10. In NaCl-2H20 each ion interacts with two counterions, so that we would expect the charge-charge contribution to be negative. In both crystals the water is surrounded by an approximate tetrahedron of ions, the protons pointing towards two anions and the lone pairs towards cations. The water's tetrahedral charge distribution is a combination of the dipole and  Chapter 5. Calculation of Phase Coexistence  138  Ref. [98]  Present Model  2.60 x 10~19  2.77 x 10- 19  -1.13 x 10" 26  2.57 x 10- 26  esu cm2)  .945 x 10" 26  -2.57 x 10" 26  Qyy (esu cm2)  .135 x 10" 26  0  Q dipole-dipole zz (esu cm )  -10.40  -10.22  dipole-quadrupole  -5.24  -12.38  quadrupole-quadrupole  -1.08  -5.00  polarized dipole (esu cm) tyxx \  2  Table 5.9: Comparison of nearest-neighbour energies with Coulson and Eisenberg. The energies are in units of kT.  c-c  c-d  c-q  d-d  d-q  q-q  NaCl-2H 2 0  -47.1  -21.8  -6.30  2.27  .601  .046  LiL3H 2 0  -34.1  -62.9  -12.7  12.9  6.26  -.572  Table 5.10: Contributions to the electrostatic interaction of hydrates, "c" refers to charge, "d" refers to dipole, and "q" refers to quadrupole. The energies are per volume, in units of solvent diameter and kT. square quadrupole, so it is reasonable that both charge-dipole and charge-quadrupole contributions are negative. Considering Figure 5.13(c) of LiI-3H20 , we observe that all dipoles are oriented towards the central axis. It is obvious that this should give a positive dipole-dipole contribution. Remembering as well that alternate H 2 0 ' s are staggered in their vertical position, we can imagine taking the circle that surrounds the Li + ions and laying it flat, as in Figure 5.15. From this picture it is clear that we would expect an attractive quadrupole-quadrupole interaction, as is the case.  Chapter 5. Calculation of Phase  Figure 5.15: LiI-3H 2 0 .  139  Two-dimensional representation of the arrangement of quadrupoles in  +  + +  5.5.2  Coexistence  +  + +  +  +  + +  +  +  Series C o n v e r g e n c e  As previously discussed, the question of convergence of the perturbative density functional theory is very much open. In this section we discuss two convergences of importance within the second-order perturbation theory, namely, convergence of the lattice vector expansion, and convergence of the rotational invariant expansion. In Table 5.11 we give an evaluation of the functional for various radial cutoff lengths (re) of the lattice vectors for ice Ih and  NaCl-2H20 . We give results both with and  without the Ewald sum technique. Clearly, when using the Ewald sum technique, convergence is quite rapid. In ice, where we have dipole-dipole interactions, the Ewald method is seen to be helpful. In systems containing ions it is absolutely necessary. It is not as easy to verify convergence for the rotational invariant expansion. Because of the effort involved in obtaining solutions for large basis sets, and because we have tried to conduct these exploratory investigations using as much as possible pre-existing liquid d a t a [8], we have not performed the density functional theory for rimax > 4. As discussed in Section 4.2, thermodynamic properties of liquid solutions have been shown to be reasonably converged by imposing n m a x = 4. This would not be obvious from calculations  Chapter 5. Calculation of Phase Coexistence  Icelh  140  NaCl-2H 2 0  re  Ewald  No Ewald  Ewald  No Ewald  1.2  -185.93  -181.61  -37.719  -71.219  1.3  55  55  -37.688  -80.469  1.4  55  55  -37.805  -36.872  1.5  55  r>  -38.163  -27.811  1.6  55  55  -38.275  -31.610  1.7  -186.04  -184.39  -38.389  -0.860  1.8  55  55  -38.553  -46.963  1.9  55  55  -38.670  -108.27  2.0  -185.99  -183.48  -38.666  -96.332  2.5  -185.95  -184.14  -38.511  4.848  2.9  -185.92  -184.09  -38.450  18.052  3.0  55  55  -38.438  4.582  Table 5.11: Convergence of lattice vector expansion for ice Ih and NaCl-2H 2 0 . Values of the functional are given for various lattice vector length cutoffs ( r e ) . Results are shown both with and without treatment of long-range correlations by the Ewald method. only of Umax = 2,3,4, as the convergence is, for some properties, somewhat abrupt [61], which lends weight to the possibility that our density functional results are also reasonably converged.  There is an argument to the contrary as well. T h e thermodynamic and  dielectric properties of these liquid systems depend essentially on projections for which m,n  < 2, which are fairly well converged for n m a x = 4. Density functional theory,  though, makes use of the infinite expansion, which, if truncated before contributions become negligible, would affect the results. Convergence for this expansion, then, is still an open question, and it would be well to investigate it in, say, freezing of dipoles. There  Chapter 5. Calculation of Phase  Coexistence  141  have been, as yet, no reliable published results which positively demonstrate the freezing of dipolar hard spheres [31,81]. Since all studies to date have used only the  first-order  perturbation theory for the liquid structure, this may be an indication of the importance of the question of invariant expansion convergence.  5.5.3  Computing Resources  The resource requirements for the calculation of liquid structures are discussed in Ref.  [7]. Two factors contribute strongly to the time requirements of the functional  evaluation in this molecular density functional theory. One is the numerical integration required to calculate W^*1  M  when // / 0, and the other is the execution of a series of  loops for performing the various sums in Section 5.1.3. In Table 5.12 we show, for ice Ih, NaCl-2H 2 0 , and LiI-3H 2 0 , the number of independent W^  ** terms to be calculated  by numerical integration, and, for a number of different lattice vector length cutoffs, the number of times a command is executed if it is at the center of the loops. In performing the sum over R, i,j,p,  q, we use a minimum image convention and the radial cutoff length  mentioned above, as this gives the most efficient use of resources. Using this method, we find that evaluating the functional for re = 1.2 requires on the order of a few seconds on a 5 megaflop SGI workstation, while using re = 3.0 requires on the order of a minute. For small re, calculation of the W™*1 M, which is done outside the loops, accounts for a significant proportion of the execution time, while for large re, it becomes unimportant. As mentioned in Section 5.1, the only other numerical integration is that over r\i, and this makes a negligible contribution. It is also the only numerical integral to be done when the system is composed of spherical and/or linear particles.  Chapter 5. Calculation of Phase Coexistence  142  Ice  NaCl-2H 2 0  LiI-3H 2 0  120  240  360  loop iterations for re = 1.2  39,824  1504  60,696  loop iterations for r e = 2.0  256,920  191,556  240,240  loop iterations for re = 3 . 0  503,620  912,132  989,220  Number of independent W^  M  , fi ^ 0  Table 5.12: Resource requirements for the molecular density functional program. Given the order of execution time for evaluation of the functional, this is then multiplied by a hundred or even a thousand for the minimization, depending on the number of variables involved. This functional is a difficult one to minimize, for two reasons. One is that large numbers of variables are often required (a full minimization of NaCl-2H20 would require thirty varibles, even for our relatively simple parametrization of the density), which rapidly increases the number of function evaluations, and the other is the "flatness" of the functional with respect to Gaussian widths and angular oscillation widths. Thus, values very close to the minimum of the functional can be obtained even without being particularly close to proper solid parameters. This is good for obtaining a value for the minimum, but frustrates determination of the parameters which give the solid structure. Very flat functions are generally difficult for minimization routines to deal with, and a number of different library routines were tried, none being completely reliable or satisfactory. Eventually, a simple "shotgun" method was developed, and was checked for certain cases against other routines. In this brute force method, a range for each minimization variable is specified, and a random number generator is used to select a number of points for evaluation. These are compared, and the ranges are gradually reduced. Minimization is considered to have been achieved when changing any variable  Chapter 5.  Calculation of Phase  Coexistence  143  Xi by a small quantity ±&r,- always yields a larger value of the functional. Thus, derivatives are not calculated, as these are deceptive indicators for a very flat function. This "shotgun" routine would be dangerous to use on an unknown function, especially if there are multiple local minima, but in this particular application it was found to be quite adequate.  Chapter 6  Calculation of S t a b i l i t y Limits  In Section 2.2 we introduced thermodynamic stability theory, and showed that the vanishing of certain derivatives signaled the onset of instability. We distinguished between the strictest stability criteria, whose collective positivity is necessary and sufficient for a system to be thermodynamically stable, and the weaker criteria, whose positivity is necessary but not sufficient for stability. The relevant derivatives for simple binary mixtures were given in terms of the integrals of direct correlation functions in Eqs (2.25)(2.31) and the importance of using strictest criteria was illustrated. In this chapter we obtain analogous quantities for electrolytes, in terms of both total and direct correlation function integrals. We then apply stability theory to a simple binary mixture of hard spheres in water and to three different electrolyte solutions.  6.1  D e r i v a t i o n of S t a b i l i t y Criteria for E l e c t r o l y t e s Kirkwood-Buff theory relates the integrals of spherically averaged total correlation  functions to the thermodynamic quantities of mixtures. After establishing some notation, we use the Kirkwood-Buff expressions to obtain a stability theory in terms of total correlation function integrals. This is convenient for relating instability to fluctuations in the liquid. We then use the Ornstein-Zernike equation to obtain equivalent stability expressions in terms of direct correlation function integrals. These allow us to relate instability to the loss of numerical solution in integral equation calculations.  144  Chapter 6. Calculation of Stability  6.1.1  Limits  145  Total Correlation Function E x p r e s s i o n s  Let us suppose that we have an electrolyte solution consisting of a compound of the form A" + B"~ in solvent. For this system we define g  =  G/N = x+fi+ + x-n-  /  =  F/N = g-pv,  (6.183)  Pi  =  p+/v+ = p-/v-,  (6.184)  x2  =  —, P  (6.185)  (*» =  /^solvent*  (6.186)  H2 =  v+p+ + v-P-,  (6.187)  v  =  + xsfx3,  V+ + V-,  (6.182)  (6.188)  where G is the Gibbs free energy, N the number of particles in volume V, X{ and //, the mole fraction and chemical potential of species i, F the Helmholtz free energy, p the pressure, and v — V/N.  Using Eqs (6.182)-(6.188) and the standard expression for dF  we can also show that  * - <£)  dF FdN N ~ N2 -pdV - SdT + £,- VidNi N  =  —pdv — sdT + (fi2 — vps)dx2.  (-pV  + £,•  fiiNi)dN  2  N  (6.189)  If we apply the stability theory derivation of Section 2.2 to Eq.(6.189) we could obtain  Chapter 6. Calculation of Stability  Limits  146  Eq.(2.18) modified for ionic systems, or  (d2r\ dv2  T,fi-2-ufi,  _  (6.190) ,dxl/T,v  T,x2  To obtain explicit expressions for the four quantities in Eq.(6.190) we also use the expression equivalent to Eq.(2.17), which is  (d*f  (6.191)  \dx2dv)  M)T,V  M)T,V  y T,x2  \dv2  Kirkwood-BufF theory supplies us with expressions for the quantities (Q^) and ( f ^ j  , ( dx gv) ,  , which are sufficient to fully determine Eqs (6.190) and (6.191). It is con-  venient to first express these three derivatives in terms of more familiar thermodynamic quantities. It can be shown that [44] X2  (6.192)  Xs  \dx2)Tp  \dxi' 2 / T . p  and we can employ Eq.(6.192) to show that  (d2g\  d  2  dJX 2 / T , p  dx2  dps  f>2 - Vfls + X2  = (p) \dx2JT,v  dx2i  + xs dx2 T,v  , JT,P'JT,p  'dps dx2, T, P (6.193)  xs\dx2)Tp' From Eq.(6.189) we have  '(Pf 2  dv ,  T,x2  'dp dv  T,x2  U— v\dPj XT,X2 '  T,x2  (6.194)  Chapter 6. Calculation of Stability  Limits  147  Eq.(6.189) and a familiar identity give us  - (dPX  d2/ y  dx2dv)  \dx*)T,v  Jdp\2  (dvV  \dvJT,X2\dx^T,P The partial molar volumes, {u,}, satisfy v = X2V2 + xavs.  (6.195)  A similar relation to Eq.(6.192)  is obtained by replacing //,- with r;,-, and this allows us to write  Eqs (6.193), (6.194), and (6.196) allow us to apply the results of Kirkwood-BufF theory to Eqs (6.190) and (6.191). By means of Kirkwood-BufF theory [51] modified for ionic solutions, Kusalik and Patey [9] showed that  _ V.  =  1 + ps(G+- + Gss — 2G+S) 1 + ps{Gss - G+s) (*[! + Ps(Gss + G + _ - 2G+S)Y G+- — G+s 1 + ps(Gss + G+- — 2G+S)'  \dPiJT,v  (6.197)  lD,iy8j  (6.199)  Pi{G+--G+s)  where Gij = J dr_f dQhij(r_,Q). One may similarly show that  \dp2jTtP  Ps  \dp2)TiP  and  ^w)T,r---^X,r^J,.-o,,Y To convert the { I ——- 1  >dPiJT,P)  > into derivatives with respect to x2 we write  (6 20i)  -  Chapter 6. Calculation of Stability Limits  dp*.  148  'dp2 dp 2 / T . p  T,p]  ' dx 2 X2  +p  9fi2  T,P  Rearranging this result and combining it with Eqs (6.200) and (6.201) gives  (dx2 X c  2  /T,p  0  dps  'dp2 ,dn2  x2\ dp,2 /T,p  T,p  dpA Kdx2JTp  1 + G+--  x2 p2ps(ys+Gss  2G+S)'  (6.202)  We can now use Eqs (6.193), (6.194), and (6.196), and the Kirkwood-BufF expressions, Eqs (6.197)-(6.199), and (6.202), to solve Eqs (6.190) and (6.191). Noting that p = l/v we obtain <Q2£ dv2, T,x  2  1 SG f3v2DG'  (6.203)  s  2 fl dJJ_ 2  dv  1 2  ' T,H2-VV3  Kdx2)Tp  PV 1  SQ/DG +  X2XSV3/DG'  1 1 /3x2xs SG ' 1 1+ fix2xs  (6.204) (6.205)  x2xsV2/DG SG  (6.206)  where DG  — p2G+-(l  + psGss) -  p2psG+s,  SG  =  [p2(l + psGss) + P2ps(G+- -  VG  =  [l +  (6.207) 2G+s)]/p,  (6.208)  pt(Gaa-G+,)]-uP2[G+.-G+s].  (6.209)  The results above depend upon the three quantities G+-, G+s, and Gss, but the other Gij do not appear explicitly. This is because electroneutrality requires that G+- — G++ + — — G— -\ P+  , P-  (6.210)  Chapter 6. Calculation of Stability  Limits  149  G+s = G-s.  (6.211)  Results for a non-ionic mixture are as in Eqs (6.203)-(6.206) with v = s — 1 and + P2G22)-piP2G212,  DG  =  (l + piGn)(l  SG  =  l + pxix2(Gn  VG  =  Pi(G11-G12)-P2(G22-Gl2).  (6.212)  + G22-2G12),  (6.213) (6.214)  The quantities Gij can be readily obtained after any integral equation calculation of the liquid structure. We will return to the question of stability, but first it is useful to consider quantities in terms of the direct correlation function.  6.1.2  D i r e c t Correlation Function E x p r e s s i o n s  For nonspherical particles the Ornstein-Zernike (OZ) equation takes the form of a matrix equation [58] giving exact relations between the expansion coefficients of total and direct correlation functions. Relationships for the {Gij} are obtained in the limit k —* 0 of the transformed OZ equation, and these expressions can be simplified by noting that f£%(k  = 0) = 0 for all f™%(r) which have / ^ 0 and which decay faster than r" 3 .  Using such methods Kusalik and Patey [9] have shown that for a solvent with dipole and tetrahedral quadrupole, as is the case here, the following relations hold G+S = G_S  n  Gss  = =  {U+ s  + V C s)  (P+G+S  + P-C-S)G+S + Css  ^ (1 -  - ; P2G+^ psCss)  ji T7-, (1 — psUss)  (6.215) ,  /fi01fiv  (6.216)  where C,j = 47r/0°° CQQ°J(r)r2c?r is analogous to the definition of Gij. From Eqs (6.215) and (6.216) it is clear that if we can obtain G+- in terms of the {C,j}, that we will be able to give the stability expressions in Eqs (6.203)-(6.206) in terms of direct correlation function integrals.  Chapter 6. Calculation of Stability  Limits  150  From the matrix OZ equation of Ref. [58] we extract the expression  ::}(:::) j ^ - r O ^ , ie.m)  C-C-I>Z:£JT{: a  where i,j, k represent ions, I  b  c  > is a Q-j symbol [75], and the tilde represents the  Hankel transform defined in Section 5.1.10. For the case of spherical ions in a tetrahedral quadrupole solvent it can be shown [9] that in the limit k —• 0 Eq.(6.217) can be simplified to i.  V1  ~  1  ~  hij - Cij = l^pkfiikCkj  1000 -000  i  + pshQQtiscmjs  k  *  1011 -Oil  - -psli00tisc00js.  ic 01 o\  (6.218;  6  The small k expansions of charge-charge and charge-dipole quantities are given by [99]  ^  = -inpqiqjk-2  c°olls =  + 4 0 ) + •" •  (6.219)  -i(-^Pqi»k-1+c°0Q1}s1)k  + ..-)  (6.220)  h\f + h\fk2 + ---  (6.221)  h°ol% = - » $ S 1 ? ) * + •••),  ( 6 - 222 )  hij  =  which may be substituted into Eq.(6.219). Equating coefficients of k~2 gives the electroneutrality conditions of Eq.(6.210), and equating constant terms gives  Ga ~ ~cf = EPk[G*$ ~ *2 Wfltfi)] + PsGisCjs - \PshZV(^PWs)6  k  (6-223)  We multiply by pj and sum over the index j . The neutrality condition J2j PjQj — 0 simplifies this result to  E PiGa = Eft* + PkGik] E Pi*8 + PsG,s E P&..  (6.224)  Using Eqs (6.210) and (6.215) to substitute for G„, G t s we can obtain the final result p2G+_  =  l ^ f i , L>c  (6.225)  Chapter 6. Calculation of Stability Limits  151  where Dc  = (1 - p.Cu){y  - p2Cu) - p2p.C]„  (6.226)  (6-227)  C„ = 5>.^-c£\ ij  Cu  = X>C,S.  (6.228)  i  From Eqs (6.215) and (6.225) it follows that G+s = G.S = ^  (6.229)  and, with Eq.(6.216), l+PsGss=V-=-^-.  (6.230) *->C  Substitution of Eqs (6.225), (6.229), and (6.230) into Eqs (6.207)-(6.209) yields DG  = l/Dc,  (6.231)  SG  = Sc/Dc,  (6.232)  VG  = Vc/Dc,  (6.233)  where Sc  = l-p(x2Cn  + x2sCss + 2x2xsCIs)  (6.234)  Vc  = vpsCss - p2Cu + {vpv - ps)Cn  (6.235)  We are now able to obtain stability indicators in terms of the {C,j}. Substitution of Eqs (6.231)-(6.233) into Eqs (6.203-6.206) gives us the desired result  =  (6236)  (0L £**  (d2f'\  1 (d2g\ <^XVT,P  'd2f\ 9xlJT„  1  -  Sc DC  fix2xsSc' 1 Dc (l+x2xsV£/Dc). f3x2xs Sc  (6.237) (6.238) (6.239)  Chapter 6. Calculation of Stability  Limits  152  In displaying numerical results on stability, it is convenient to divide the quantities in Eqs (6.236)-(6.239) by their ideal counterparts. Ideal systems have no phase transitions, and hence, no instabilities, so the qualitative nature of the figures is unchanged. We have for the ideal gas the familiar expressions p = kT/v,  (6.240)  Hi = H°i(T,p) + kTlnxu  (6.241)  where /x°(T,p) is the chemical potential of the pure component, and depends only on T and p. From Eqs (6.240) and (6.241) one may readily derive the relations  \dv2)T \  \dv2)T / l,tl2 — VHs 2 id  /^2'  ^  2  / 1 ,X2 id  (d g\  fd f\  \OX2JTp  A) \OX2J  "  1 pXiXs  Tv  We therefore define the following quantities for plotting: J X  '  =  Sc  =  'd2f'Y dv2  Sc l  J r,„-w  + W.VS/Dc  (aPgY  £c_  \dxl)Tp  Sc  ^2  =  -z-V+W.VS  7T,  (6-242)  SG  1  DG 1 +  x2xsV3/DG'  (6.243)  j _  (6.244)  SG' Dc)  =  =  •  The asterisk indicates that a quantity is divided by its ideal counterpart.  6.245) These ex-  pressions are the same as those obtained for a simple binary mixture, when appropriate definitions of Do, SG, VG, DG, SG, and VG are used, and v = s = 1. Eqs (6.242) and (6.243) are indicators of mechanical stability and Eqs (6.244) and (6.245) are indicators of compositional stability. When a phase is stable or metastable, the  Chapter 6. Calculation of Stability  Limits  153  derivatives in Eqs (6.242)-(6.245) will all be positive. If any of the derivatives vanishes, this indicates that the system has become unstable. However, such behaviour will always be subject to the restriction that Eq.(6.242) is greater than or equal to Eq.(6.243), and that Eq.(6.245) is greater than or equal to Eq.(6.244).  Chapter 6. Calculation of Stability  6.2  Limits  154  H a r d S p h e r e and C s l S o l u t i o n s In this section we consider two systems. The first is composed of hard spheres (HS)  in water at 25° C, where the water is represented by the model given in Section 4.2. The hard sphere diameter is equal to the water diameter, and the cross-term in the pair potential is also a hard sphere interaction, so that the hard spheres are expected to be hydrophobic. The density is fixed at p* — .732, very nearly the density of real water at 25°C and 1 atm, and the mole fraction of the hard spheres is varied up to ~ .018 (or ~ 1 mol/L), where an instability occurs. The other system is an aqueous solution of Csl at 25°C, and we use two different models, CHS and CHSR9, which were also introduced in Section 4.2. It is convenient to discuss the HS and Csl systems together, as their behaviour is similar in certain respects, and we shall discuss KC1 and NaCl in the following section. In Figures 6.16 and 6.17 we plot strictest mechanical and compositional stability indicators, as well as the quantities DG, SG, DC, and Sc for hard spheres in water and for CHS and CHSR9 models of Csl in water.  According to Figures 6.16(a),(b) and  6.17(a),(b) both mechanical and compositional stability indicators are signaling instabilities. We see this from the fact that at high concentration they are small and rapidly decreasing. Of course, because of the numerical difficulties which occur as one approaches an instability, we do not know for sure whether these quantities do actually vanish, but we assume that the displayed behaviour will be continued. The choice of interionic potential affects the location of the Csl instability, though not as strongly as the location of the coexistence predicted by density functional theory. The essential behaviour of the two models appears to be the same, and we continue with only the CHS model in the remainder of this section. From Figures 6.16(c),(d) and 6.17(c),(d), we see that the quantities DG and SG are  Chapter 6. Calculation of Stability  Limits  155  Figure 6.16: Stability indicators of hard spheres in water. The points are shown explicitly for the first graph, (a) Defined in Eq.(6.243). (b) Defined in Eq.(6.244). (c) Defined in Eq.(6.212). (d) Defined in Eq.(6.213). (e) Defined in Eq.(2.25). (f) Defined in Eq.(2.28).  Chapter 6. Calculation of Stability  Limits  Chapter 6. Calculation of Stability  Limits  157  Figure 6.17: Stability indicators of Csl in water for both the CHS model and the CHSR9 model with A* = .7. The points are shown explicitly for the first graph, (a) Defined in Eq.(6.243). (b) Defined in Eq.(6.244). (c) Defined in Eq.(6.207). (d) Defined in Eq.(6.208). (e) Defined in Eq.(6.226). (f) Defined in Eq.(6.234).  Chapter 6. Calculation of Stability Limits  d'f'V  (a) S*  158  (b) (  v  '  x  iV  CHSR9  2-  /^~~S. 2-  1-  n_  u —|  10  )  (  1  5 mol/L  10  SG  (d)  CHS  CHSR9  321 1 1 ,  1n un ()  sc  (f)  (e) CHS  CHSR9  1 *? 1Z  ZU  1  5 mol/L  CHS  10  CHSR9  10i  i  /  /  \  .  '1  10-  86-  u "1 0  I 5 mol/L  I  10  0  5 mol/L  10  Chapter 6. Calculation of Stability  Limits  159  both divergent, although, according to Eq.(6.243), the divergence of DQ must be faster than that of SG- DG and SG are composites of the {Gij}, which diverge when the total correlation functions become long-range. Looking at Eq. (6.208) we see that divergence of SG can be an indication of a tendency to demixing, as it can result from the development of long-range attractive correlations between like components and long-range repulsive correlations between unlike components, as discussed in Section 2.2. It is important that we refer to "components" rathers than "particles". The discussion in Section 2.2 is illustrated by a simple binary mixture, but these concepts apply as well to ionic systems in which the electrolyte component consists of two types of particle, which because of electroneutrality cannot become macroscopically demixed. Therefore, we are seeing increasing correlation between all particles of the electrolyte, and decreasing correlation between each particle of the electrolyte and the solvent particles. The onset of instability is straightforwardly related to the interruption of numerical procedures since both are produced by the vanishing of Dc, for both components. We can see this from Eqs (2.22)(2.24), (6.225), and (6.243)-(6.244). Also, it is clear from Eqs (6.243) and (6.244) that Dc must vanish more rapidly than Sc- In fact, Sc need not vanish at all as appears to be the case in the hard sphere system. This is important, since, according to Eq.(6.242), Sc is equal to a weaker criterion of mechanical stability. Thus, if we considered only the weaker criteria displayed in Figure 6.18, we would conclude that the HS system is becoming compositionally but not mechanically unstable, and that the Csl system is becoming mechanically but not compositionally unstable.  These examples underline  the importance of referring to strictest stability criteria. Returning to the discussion of liquid correlations, we can show how these become long-range near the instabilities. In Figures 6.19-6.21 we plot the angle-averaged total correlation functions multiplied by r 2 . The r 2 factor arises from the spherical symmetry of the function, and helps to amplify any long-range effects present. In each graph there  Chapter 6. Calculation of Stability  Limits  160  Figure 6.18: Weak stability indicators of hard spheres in water and Csl (CHS model) in water, (a) Defined in Eq.(6.242). (b) Defined in Eq.(6.245). (c) Defined in Eq.(6.242). (d) Defined in Eq.(6.245).  Chapter 6. Calculation of Stability  Limits  (c) 1n  8-  6-  A —  4 -)  0  |  5 mol/L  1)  161  Chapter 6. Calculation  of Stability  Limits  162  are two lines, one far from instability and the other close to it.  The principal  feature to be observed in Figures 6.19-6.21 is the upward shift of correlations between like components, and the downward shift of correlations between unlike components. In the HS system this is most evident in the HS-HS and HS-H2O correlations, but there is a perceptible upward shift of the entire H2O-H2O function as well. The total correlation function is asymptotically equal to the negative of the potential of mean force [47], so that what is illustrated in these figures is in fact a long-range hydrophobic interaction between HS and water, such as would be experienced by nonpolar particles in aqueous solution. Although the HS system is somewhat aside from the central occupation of this study, it is included here, both because of its intrinsic interest, as well as because its interpretation is straightforward and thus guides us in the Csl analysis. In the case of the Csl solution, the signals of demixing are again more readily perceived in some projections than in others. Therefore, the three functions involving solvent, shown in Figure 6.21, are clearly oscillating about a non-zero value at the higher concentration. However, from the ion-ion correlations in Figure 6.20 we are unable to clearly ascertain such behaviour, which, if it exists, is obscured by significant structure extending to considerable distances. It is instructive to consider a linear combination of these quantities which is given in Figure 6.22(a).  This quantity is related to the  probability of finding one ion (of any charge) a given distance from another ion (of any charge), and, near a compositional instability, would be expected to show a long-range attraction of ions for themselves. It does appear to do this, at least more convincingly than the individual correlations in Figure 6.20. In Figure 6.22(b) we plot a linear combination related to the probability of finding an ion (of any charge) the given distance from a solvent particle. We see here the long-range repulsion expected near a compositional instability. In Figure 6.23(a) we combine all of the correlation functions of Csl to give a quantity related to the probability of finding any particle a given distance from any  Chapter 6. Calculation of Stability  Limits  163  Figure 6.19: r 2 times the spherically averaged total correlation functions for the hard spheres in water. 1 denotes hard spheres and 2 denotes the solvent. The higher concentration is near the stability limit, (a) The HS-HS correlation, (b) The H2O-H2O correlation, (c) The HS-H2O correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  Limits  164  r2h22{r)  1.0 M  n A u. 4  0.2-  0>)  1.032 M  /%  n n u. u -0.2n  A  U. 4  2  I  i  3  4  5  r/d}i3o  r% 2 (r)  1.0 M  1.032 M  4  1  \  0.5-  (c)  f\  u  -  // -0.5— i "i  0  i  i  2  4 r/dH,o  6  Chapter 6. Calculation of Stability Limits  165  Figure 6.20: r 2 times the ionic total correlation functions for Csl (CHS model) in water. The higher concentration is near the stability limit, (a) The cation-cation total correlation, (b) The anion-anion total correlation, (c) The cation-anion total correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  Limits  166  Chapter 6. Calculation of Stability  Limits  167  Figure 6.21: r 2 times the total correlation functions involving solvent for Csl (CHS model) in water. The higher concentration is near the stability limit, (a) The cation-solvent total correlation, (b) The anion-solvent total correlation, (c) The solvent-solvent total correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  Limits  168  .2 LOOO  '  8.8486 M  1.0 M  o _  Z  1"  (c)  'A  / \  n u  1 -  I  •  2  1  4 r/dH7o  Chapter 6. Calculation of Stability  Limits  169  Figure 6.22: r 2 times composite total correlation functions for Csl (CHS model) in water. The higher concentration is near the stability limit, (a) The ion-ion total correlation, (b) The ion-solvent total correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  (a)  Limits  170  r= [ft++(r) + 2/i+_(r) + /i__(r)] 1.0 M  8.8486 M  4-  2-  -2-  -4  I 4 r/duto  2  r 5  i  6  r 2 [CW + C W ]  ^  1.0 M  8.8486 M  "7  3  2-  #  / 1  M i I I i  / \ 1 \ / I  ''  1 -  /I  t  \  I  \  '  A/  1  0.V 1  -1 -  1  '  '  >  1  /  \  \ '' 1/  I  •<M  1  /  2  I  I  I  3  4 r/dHio  5  I  6  7  Chapter 6. Calculation of Stability  other particle.  Limits  171  We argue that the long-range nature of this function is analogous to the  long-range nature of a pure fluid near an instability, and therefore is consistent with and a signal of the onset of mechanical instability. To introduce the final graph in Figure 6.23(b), we combine Eqs (6.208), (6.210), and (6.211) to give  SG  =  v- + i: + PI(G++ + G— +  +-)  /92(r + psGss) + ps—— i/+ • + v-' Jz) . P2Ps to + tti  =  2G  IG++ + 2G+. + g - _  p l G++ + 2G+_ +G__ _ -V + px2xs 4  {G+s  + Q_s) +  {G+s  Qs  p2ps(G+s - G-s) IP + G_a) + Gg  if v+ — V-.  In Figure 6.23(b) we plot r 2 times the function whose integral is equal to  (SG -  ^l^)l{pX2Xs).  This shows explicitly how the divergence of SG arises from the long-range correlations. Given that the more concentrated Csl solution is well above the experimental saturation concentration, it is tempting to investigate the structures for signs of crystal nuclei. In Table 6.13 we gather information relevant to this question. In the first column we give the locations of various peaks in Figure 6.20, along with a description of their shape. In the second column we list the distances at which we could expect to find peaks associated with certain configurations in the liquid. These are all comprised of the two correlated particles separated by a third particle, with the bridged configurations being non-linear triples, and the others being linear triples. The third column contains distances at which we could expect to find peaks associated with certain configurations in a Csl crystal. Values in the three columns are arranged so that values at the same horizontal level are comparable. The location of most peaks is clearly best described by liquid configurations. Not only are the distances more accurate, but the smoothness of the first peak of Figure 6.20(a) and (b) is better explained by a continuum of bridged configurations  Chapter 6. Calculation of Stability  Limits  172  Figure 6.23: r 2 times composite total correlation functions for Csl (CHS model) in water. The higher concentration is near the stability limit, (a) The particle-particle total correlation, (b) The function whose integral is equal to (SG — 1 /'v)I'(px2%s). Distance units are in solvent diameters.  Chapter 6. Calculation of Stability Limits  173  (a) r* {*i[fc ++ (r) + 2 i + . ( r ) + * - W ] + 2x 2 x s [Ag? +5 (r) + C ? - s M ] + * s £ s ( r ) } 1.0 M  1 _ 1  8.8486 M  1 • 1 *  0.5-  «  i * 1 f •  A # 1 / \ /s  i  / W i  •  1*1 * * 1  /  1 |  n U  _  \ \  l  1 f 1 t \ t 1 $ \ 9 \ 1 \ 1 \ t  I'  /  /i ** ** /1 . ' /  \/  /  \  *  /  \  ^ s *-* •  ,  —'  ™~  /'''  vv  \ \ y/ y  ••  -0.5-  4  1 ~"|  i  i  i  i  2  3  4  5  6  »"/^HjO  (b)  r 2 {[/> ++ (r) + 2A + -(r) + h..(r))/4  - [h™+s(r) + h™_s(r)] + 1.0 M  h™s(r)}  8.8486 M  z 11 11  1-  i  I*/  \  /.  Mi  /I  r 1 \ '1 •  \''*-A 0-  • I f ' i f  V  /  v  '  -1 '•  Z ~]  i  i  i  i  2  3  4  5  r/dH,o  6  Chapter 6. Calculation of Stability  Limits  174  Peak Positions  Liquids Configurations  Crystal Configurations  1.69 round 2.31 sharp 2.77 sharp  bridged +s+ 2.28  + - + 2.72  1.57 2.22 2.72  h._  1.77 round 2.49 sharp 2.74 sharp  bridged -s- 2.44 - + - 2.72  1.57 2.22 2.72  h+-  2.37 sharp 2.64 sharp  h++  +s-  2.36 —  —  2.60  Table 6.13: Peak positions and shapes refer to Fig. 6.17. Bridged liquid configurations refer to the two correlated particles in a continuum of rolling contacts with a third particle. Non-bridged liquid configurations refer to a linear arrangement of the three particles shown. Crystal configurations refer to possible distances between the two correlated particles in a close-packed Csl crystal. than by a fixed configuration in the crystal. The 2.72 distance appears in both liquid and crystal, and is the same configuration in both. This alone then cannot be evidence of crystal nuclei if peaks are absent at 1.57 and 2.22. We conclude that the formation of crystal nuclei, even if it is occurring undetected, is not the explanation for the long-range behaviour of these solutions. Instead we are seeing the Csl interacting hydrophobically with water to induce a compositional instability, while at the same time undergoing a mechanical instability. In this sense the mechanism is similar to the demixing of hard spheres and water, and may be termed "hydrophobic".  Chapter 6. Calculation of Stability  6.3  Limits  175  KC1 and N a C l S o l u t i o n s In this section we examine the stability of model NaCl and KC1 solutions at 25° C. We  employ the CHS model for both compounds, and also CHSR9 model for KC1. Looking at Figures 6.24 and 6.25 we immediately see behaviour very different from that in Section 6.2.  From the rate of decrease at high concentration of the quantities plotted in Fig-  ures 6.24(a) and 6.25(a) we would conclude that the systems are becoming mechanically unstable, but the actual values of the indicators at loss of numerical solution are much larger than in Figures 6.16(a) and 6.17(a). In Figures 6.24(b) and 6.25(b) the indicators of compositional stability actually appear to diverge at high concentration, rather than go to zero. The system does not at all appear to be heading towards compositional instability. Because there is no compositional instability, we would not expect to see a divergence of SQ, and indeed there is none, nor for Do-  Instead, from Eqs (6.243) and (6.244)  we would expect to see SG —• 0, and this is consistent with the Figures 6.24(c),(d) and 6.25(c),(d). Instead of the growth of long-range demixing correlations, we see a balancing and cancellation of the terms comprising SG (see Eq.(6.208)). Looking at Eqs (6.243) and (6.244) we can deduce that, in the absence of divergent Dc, we must have Sc —> 0 while Dc either remains finite, or vanishes more slowly than Sc- By reference to Figures 6.24 and 6.25 we confirm that Sc vanishes while Dc appears to remain finite. There is yet another important aspect to this picture, though. include the detail of Dc and ( g ^ J  near the instability.  In Figure 6.26 we  This figure suggests that  certain changes may be happening very quickly near the stability limit which allow Dc to vanish and these systems to become compositionally unstable. We hesitate to base any firm conclusions on these results, though they may be significant, since, as is evident from Figure 6.26, numerical convergence of the RHNC equations near these points becomes  Chapter 6. Calculation of Stability  Limits  176  Figure 6.24: Stability indicators of KCl (CHS) in water and KCl (CHSR9) in water. Points are shown explicitly for the first graph, (a) Defined in Eq.(6.243). (b) Defined in Eq.(6.244). (c) Defined in Eq.(6.207). (d) Defined in Eq.(6.208). (e) Defined in Eq.(6.226). (f) Defined in Eq.(6.234).  Chapter 6. Calculation of Stability  w m,„-  Vlis  °  Limits  CHS  8  (c)  DQ CHS  0.15  CHSR9  0. 10-  0.05  0.00  5 mol/L  (e)  DC  10  CHSR9  CHS  /  20-  /  10•  n  ^  _  ^  ^  I  0  5 mol/L  10  177  Chapter 6. Calculation of Stability Limits  178  Figure 6.25: Stability indicators of NaCl (CHS) in water. Points are shown explicitly in the first graph, (a) Defined in Eq.(6.243). (b) Defined in Eq.(6.244). (c) Defined in Eq.(6.207). (d) Defined in Eq.(6.208). (e) Defined in Eq.(6.226). (f) Defined in Eq.(6.234).  Chapter 6. Calculation of Stability Limits  U2)T>2-^5  (a) A  A  4.2-  4-  0  10  SG  (c) n  5 mol/L  1  =L  U . ID  0. 1 0 -  0  i  5 mol/L  10  179  Chapter 6. Calculation of Stability Limits  Figure 6.26: Detail of Dc and ( g - ) * Figure 6.24 and 6.25.  180  near the instability for the solutions shown in  Chapter 6. Calculation of Stability  (a)  DC  (KC1,CHS)  20-1  1  <i  ( b ) ©r,„ (KC1 - CHS ' 3. 1  e-e^  e—  3.05  19.5-  19-  1 ft ^ -  10.31  1  9  (c)  181  Limits  1  9.005 9.01 mol/L  DC  1  2.95 9  9.015  9.005 9.01 mol/L  9.015  (KC1,CHSR9)  28.9  28.8-  28.7-  28.6 7  (e)  7.05  7.1 mol/L  7.15  ( f ) (S)r,, (NaC1' CHS)  DC (NaCl, CHS)  2.39  16.616.55-1 16.5(> 1 fi A.*,  i  9.2  9.3 mol/L  9.4  2.382 9.35  1  1  9.36 9.37 mol/L  9.38  Chapter 6. Calculation of Stability  Limits  182  increasingly unreliable. Next we consider the correlations of the KCl (CHS) solution (for brevity we do not include KCl (CHSR9) and NaCl (CHS)) and we see in Figures 6.27 and 6.28 that all the correlation functions oscillate about zero to all distances shown. do the combinations given in Figures 6.29 and 6.30(b) as well.  Not surprisingly, so This is consistent with  the non-divergence of SG and the absence of a compositional instability. In Section 6.2 we argued that mechanical instability is related to long-range behaviour of the function plotted in Figure 6.30(a).  The fact that Figure 6.30(a) does not display long-range  character, though, could be related to the fact that the mechanical stability indicator in Figure 6.24(a) is still relatively large. There are other differences as well between these graphs and those of the Csl system. The unusual structure which persists to considerable distances in Csl is largely absent here, and correlations show a regular, periodic decay after the first few peaks. The assignment of the initial ion-ion peaks follows a similar pattern to those in Csl (i.e., they are well described by common liquid configurations), but the solvent-separated peak in /i+_ is much depressed in the Csl solution relative to the KCl (this peak is at r = 2.37 in Figure 6.20(c) and r = 2.15 in Figure 6.27(c)). This suggests that the hydrophobicity of the previous system is not present here. Indeed, from the correlations of the KCl system, we find no evidence for the existence of anything like a long-range hydrophobic force. The behaviour of the NaCl and KCl systems is clearly very different from that of the Csl system, and we can perhaps explain the origin of this difference. We have modelled the ions as hard spheres with point charges at the centre, and these charges interact with the dipole and quadrupole of the solvent. The strength of interaction at contact decreases as the radius increases, so that the overall ion-solvent interaction is stronger for NaCl and KCl than for Csl. By way of analogy we can consider a simple liquid binary mixture of spherical particles. If the pair potential between unlike species is less attractive than  Chapter 6. Calculation of Stability  Limits  183  Figure 6.27: r 2 times the ionic total correlation functions for KCl (CHS) in water. The higher concentration is near the instability, (a) The cation-cation total correlation, (b) The anion-anion total correlation, (c) The cation-anion total correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  Limits  184  Chapter 6. Calculation of Stability  Limits  185  Figure 6.28: r2 times the total correlation functions involving solvent for KC1 (CHS) in water. The higher concentration is near the instability, (a) The cation-solvent total correlation, (b) The anion-solvent total correlation, (c) The solvent-solvent total correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  r.21  Limits  186  000  V + » M ,0M  (a)  9.013 M  1 -  -1  r h00tSS(r)  M  9.013 M  1 -  (c) -1  2  4 r/dH3o  Chapter 6. Calculation of Stability  Limits  187  Figure 6.29: r 2 times composite total correlation functions for KCl (CHS) in water. The higher concentration is near the instability, (a) The ion-ion total correlation, (b) The ion-solvent total correlation. Distance units are in solvent diameters.  Chapter 6. Calculation of Stability  (a)  Limits  188  r2[ft++(r) + 2/i+_(r) + /l_.(r)] 1.0 M  9.013 M  8  (b)  i  i  6  7  r2lh^+,(r) + h™_,(r)] 1.0 M  9.013 M  2-  1 -  -1 -  -2  -r 2  i  4 r/dHlo  6  Chapter 6. Calculation of Stability  Limits  189  Figure 6.30: r 2 times composite total correlation functions for KCl (CHS) in water. The higher concentration is near the instability, (a) The particle-particle total correlation, (b) The function whose integral is equal to (SG — l/v)/(pX2Xs). Distance units are in solvent diameters.  Chapter 6. Calculation of Stability Limits  (a)  190  r2 {xl[h++(r) + 2h+.(r) + h..(r)] + 2x2xs[h™+s(r)  + h™_s(r)] +  1.0 M  h™s(r)}  9.013 M  0.5-  -0.5-  -1  T  5  2  (b)  r/dHio  r 2 {[h++(r) + 2& + _(r) + A_.(r)]/4 - [Ag» s(r) + C ° - s ( 0 ] + C ° s s ( r ) } 1.0 M  9.013 M  1.5  1-  0.5-  -0.5-  -1  -i  r  3  4  r/du2o  5  Chapter 6. Calculation of Stability  Limits  191  that between like species, it is possible that at some temperature the liquid will demix into two phases of different concentrations. If on the other hand, the attraction between unlike species is greater than that between like species, then the liquid is much less likely to become compositionally unstable. In aqueous electrolytes, the interactions are composites of various contributions, and we cannot easily say how the strength of the ion-solvent interaction compares to that of the ion-ion and solvent-solvent interactions. We can only say that the ion-solvent interaction of systems with small ions is stronger than that of systems with large ions, so that if hydrophobic behaviour is going to be displayed in ionic systems, then large ion salts are more likely candidates than small ion salts. As mentioned in Section 5.4.2, alkali halides form crystal hydrates only when at least one ion is small. This may be interpreted as an indication of the difference in ionsolvent interaction strengths of various real salts, so it is possible that the hydrophobic and non-hydrophobic behaviour displayed by our models near instability do occur in real solutions as well. While it is easy to attribute the source of the difference between the Csl system and the KCl, NaCl systems to difference in ion size, it is somewhat more difficult to say what the exact nature of this difference is. To answer this question would require a better knowledge of KCl and NaCl very near the instability, but we can state some of the possible answers. If the trends shown in the detail in Figure 6.26 are real, then we can imagine a very rapid onset of compositional instability in the KCl and NaCl systems. In this case, the difference between these systems and the Csl system would be one of degree rather than a fundamental one. All systems would become compositionally unstable, but some would do so very rapidly and others more gradually. Another possibility is that the trends in Figures 6.24 and 6.25 are correct, and that the KCl and NaCl solutions become mechanically unstable, but not compositionally unstable. The correlation functions would not necessarily become long-range in this case. Still another possibility is  Chapter 6. Calculation of Stability  Limits  192  that no stability limit is reached, but that the exact RHNC solutions terminate while the solution is still stable. This type of behaviour has been observed in other integral equation calculations near instability [52]. Such a result leaves unanswered the separate question of how an exact description of the model (as opposed to the approximate RHNC description) would behave. The principal result of this chapter is that we have used the strictest criteria of stability to show that there are at least two very different types of behaviour near the limit of supersaturation for certain models of aqueous alkali halide solutions. One kind of behaviour, of which hard sphere solutes and Csl are examples, can be characterized as hydrophobic. Such systems become mechanically and compositionally unstable, and show long-range attractions and repulsions characteristic of demixing. The other kind of behaviour, of which NaCl and KCl are examples, has a very narrow region near the limit which we have not managed to characterize completely.  Prior to entering this  small region such systems cannot be termed hydrophobic, as they do not develop longrange correlations and do not become compositionally unstable. Instead, they appear to approach mechanical instability, and this in such a way that solvent and solute are still closely associated.  Chapter 7  Conclusions  We have developed a molecular density functional theory of crystallization and a thermodynamic stability theory for electrolytes, and have applied them to solutions of alkali halides in water. These applications were motivated by the availability of solutions of accurate integral equations [7,8] for two different models of these mixtures. In the simpler model we represented a molecule in solution by a hard sphere with an embedded multipole, so that the ions were charged hard spheres (CHS). The other model considered was the same, except that the ion-ion potential, in addition to the charged hard sphere components, contained an r - 9 repulsion (CHSR9). The magnitude of the repulsive term was controlled by an adjustable parameter, the manipulation of which produced significant changes in thermodynamic properties and ion correlations in the liquid. Activity coefficients, for instance, could easily be made to agree with experiment, so that in this way at least the CHSR9 model was closer to real potentials. In order to apply density functional theory to these liquids, it was necessary to extend previous work, such as the r-space results for spherical particles [22], and the fc-space results for linear molecules [31,32,81]. The theory derived is the first r-space molecular theory, and is the first molecular theory for non-linear molecules based on spherical harmonic expansions. Two principal contributions were important to produce this real space theory. The first was to recognize the usefulness of identities involving the spherical modified Bessel function in performing a number of integrals analytically. This resulted in a derivation which is in many ways simpler than that originally given for spherical  193  Chapter 7.  Conclusions  194  particles. The second contribution was the implementation of the Ewald sum technique to perform lattice sums over long-range direct correlation functions.  This was shown  to be especially vital for systems containing charged particles. A number of other issues were discussed as well, such as the elimination of identical terms from sums in the functional, the effect of solvent polarizability on the crystallization transition, and the parametrization of angular distributions to give tractable integrals. We showed that a particular angular distribution, as well as being useful for linear molecules [32], can be applied to certain non-linear molecules as well, including H2O. Our density functional calculations lead to a number of conclusions.  One of the  questions posed in approaching the study of crystallization by a basic second-order density functional theory was whether there is sufficient resemblance between correlations in the liquid and the solid to allow any prediction of a transition at all. We conclude that the theory is capable of predicting equilibrium transitions to reasonable crystal structures for these systems, as demonstrated by the anhydrous alkali halides and by ice Ih. In the case of the crystal hydrates, although no coexistence was obtained, the functional at least possessed minima corresponding to crystal structures. This is especially important for the CHSR9 model of NaCl with a large repulsion. A physical minimum was obtained for the dihydrate, but not for the anhydrous salt, which makes sense for a solution dominated by solvent-separated cation-anion pairs. This example indicates one of the limits of second-order density functional theory. A striking feature of many of our calculations was the sensitivity to model of the functional values. Is this an intrinsic feature of the potentials themselves, as would be confirmed by a simulation of crystallization for the different models, or a consequence of the truncated density functional theory? We argued in Section 5.3.2 that such dependence can be understood in terms of the influence of model on the liquid pressure and on the near-contact portion of the direct correlation functions. This suggests that the sensitivity  Chapter 7.  Conclusions  195  is a real property of the potentials. It is certainly possible that this effect can be obscured by a truncated &-space calculation, which would not represent as accurately the nearcontact portion of the correlation, and could give fortuitously insensitive results. Our conclusion is that to calculate the equilibrium coexistence of a real system, it is important to choose the potential very carefully. This choice will also affect the quality of solid structures obtained. We have argued that the extremely small Gaussian widths and high densities obtained in this study are again due to the contact-cusp of the direct correlation function, which results from the hard core potential. Thus the agreement with typical widths of real solids was poor. When small repulsions were used, as with the anhydrous alkali halides, the agreement improved slightly. When large repulsions were used, as in ice Ih, it improved significantly.  In contrast to the positional Gaussians, the roughly  Gaussian angular distributions would not be expected to show the same dependence on the contact values of the direct correlation function, and these were, in fact, predicted to have values much closer to those of real systems. In attempting to account for all important contributions to the second-order functional, we tested the hypothesis that inclusion of a residual entropy term would be sufficient to account for the proton disorder of ice Ih. Our hypothesis was found to be false, and we conclude that the various configurations can give significantly different interactive contributions. Taking these into account would require some effort, but it is easier to consider doing this in r-space than in &-space. In r-space we could appropriately average over functional contributions from a variety of different configurations within the first few shells of a given particle, even if these do not correspond to a regular order. In fc-space, however, only ordered configurations could be averaged over, which in the case of ice Ih would leave many of even the nearest neighbour configurations unrepresented. One final practical lesson. From preliminary minimizations of NaCl-2H20 , we found that allowing particles to vary from their experimental positions had a considerable effect  Chapter 7.  Conclusions  196  on the minimum of the functional. This is of importance in studying the minimization of crystals in which particle locations in the unit cell are not all fixed by symmetry. In exploring the stability of the liquid solutions from which crystallization occurs, it soon became apparent that the approach taken in the literature to the stability of liquid mixtures was inadequate to explain the results obtained for a system such as KCl in water.  This prompted a more rigorous examination of the theory, which resulted  in the derivation of the strictest criteria of stability, and the relationships between the various criteria. These relationships have been previously given [42,43], but we consider their application to this field to be a useful contribution, particularly as they are in direct contradiction to a commonly held notion [44,45] that a mixture must become compositionally unstable before it becomes mechanically unstable. We have shown in the first example of Chapter 6 (hard spheres in water) exactly how consideration of only the weaker criteria, which are necessary but not sufficient conditions for stability, can result in this erroneous conclusion. The system of neutral spheres in water was shown to approach a stability limit near a particular concentration. As this stability limit was approached, the system developed what may rightly be termed long-range hydrophobic forces. We showed that very similar behaviour occurred as the Csl solution approached its high concentration limit. In both the hard sphere and Csl systems, the compositional instability related to hydrophobicity was attended by a mechanical instability as well. Plots of various combinations of radial distribution functions showed clearly that fluctuations in the solution were forming regions of higher and lower concentration, and of higher and lower density. All of these effects for the hard spheres in water and Csl in water were shown to be related to the breakdown of the numerical solution of the integral equations. The approach to instability of NaCl and KCl solutions was in considerable contrast to the hydrophobic cases. The compositional stability indicators of the chlorides appeared  Chapter 7.  Conclusions  197  to be diverging, rather than tending towards zero as in the case of demixing. The mechanical stability indicators were decreasing rapidly, suggesting an approach to mechanical instability, but their actual value near the loss of solution was still much larger than for Csl and the hard spheres. Thus, even if instability was imminent, it was not apparent in the distribution functions, which showed no evidence of any long-range behaviour in any of their combinations. Similarly, their was no indication of numerical failure, at least in the gross features of the curves. Examination of the detail near the instability suggested that a demixing mechanism similar to that of the hydrophobic systems could come into effect in a very small region before the instability, but this could not be ascertained clearly from the numerical results. We have then two very different approaches to instability. The obvious difference between the two classes of electrolytes is the size of ions, especially cations in this case. N a + and K + bind H2O more tightly than C s + does, and this affects the overall salt-water interaction, suggesting that the two classes of behaviour are fundamentally different. On the other hand, there is, as we have noted, an indication that demixing may occur in some extremely narrow range just before the stability limit of NaCl and KCl, suggesting that the two behaviours are fundamentally related. The question of whether the instabilities are fundamentally different, or just well separated on a spectrum of behaviours, could not be positively answered by our calculations. The results of this research suggest several avenues for future work. Building on the extensive development of hard sphere models for aqueous alkali halides [9], it would be appropriate to apply the liquid theory in a careful way to a consistent set of more realistic potentials. Our results suggest that these would give reasonable density functional predictions of real phase transitions. There is much that can be done with simpler potentials as well. A timely study that would help assess the usefulness of the molecular density functional theory is the freezing of dipolar hard spheres. This would be the easiest  Chapter 7.  Conclusions  198  system with which to study the effect of truncation of the rotational invariant expansion discussed in Section 5.5. So far freezing of dipolar hard spheres has proved elusive [31,81], but we expect an r-space method would be substantially more accurate for this system than previously employed &-space calculations. The use of more accurate liquid correlations may also be important, and would be consistent with the object of studying convergence of the rotational invariant expansion. Unfortunately, computer simulation of freezing for this system has not yet been performed, although such a simulation does not presently seem unreasonable to attempt. Modifications to the present theory are also a direction which may be pursued. Thirdorder corrections using triplet correlations are not a practical alternative, but extension with one of the nonperturbative techniques discussed in Chapter 2 is certainly possible. A completely different approach to improving the predictions of density functional theory has quite recently been suggested by Baus [100]. In this method, the results of hard sphere freezing, which have been successfully calculated by a number of theories, are used as a reference for the freezing of other systems, analogous to the methods of liquid perturbation theory [47]. This recent technique is successful in prediction of soft sphere freezing, and would perhaps lend itself to transitions involving multipolar hard spheres. Turning to stability studies, it would be of interest to study the strictest criteria of stability for simple binary mixtures of spherical particles with different strengths of cross-interaction. A number of such systems have been investigated [46], but not with the strictest stability criterion for mechanical stability. The advantage of these systems is that the region close to the stability limit could be investigated in much more detail than was possible here.  Bibliography  [1] Yang, C.N., Phase Transitions and Critical Phenomena vol. 1, eds Domb, C , Green, M.S., Academic Press, London, 1972; Dresden, M., Physics Today 41, 26 (Sept.1988). [2] Yang, C.N., Lee, T.D., Phys. Rev. 87, 404 (1952). [3] Lee, T.D., Yang, C.N., Phys. Rev. 87, 410 (1952). [4] Stewart, I., New Scientist  131, 29 (July 13,1991).  [5] Uhlenbeck, G.E., Fundamental Problems in Statistical Mechanics, vol.2, ed. Cohen, E.G.D., North-Holland, Amsterdam (1968). [6] Widom, B., Science, 157. 375 (1967); Longuet-Higgins, H.C., Widom, B., Molec. Phys., 8, 549 (1964). [7] Kusalik, P.G., "The Structural, Thermodynamic and Dielectric Properties of Electrolyte Solutions: A Theoretical Study", Ph.D. Thesis, University of British Columbia, 1987. [8] Kusalik, P.G., Patey, G.N., J. Chem. Phys. 88, 7715 (1988). [9] Kusalik, P.G., Patey, G.N., J. Chem. Phys. 86, 5110 (1987). [10] Haymet, A.D.J., Ann. Rev. Phys. Chem. 38, 89 (1987). [11] Baus, M., J. Stat. Phys. 48, H29 (1987). [12] Hohenberg, P., Kohn, W., Phys. Rev. 136, 864 (1964). [13] Mermin, N.D., Phys. Rev. 137, 1441 (1965). [14] Evans, R. Adv. Phys. 28, 143 (1979). [15] Evans, R. Liquids at interfaces (Les Houches, Session XLVIII, 1988), eds Charvolin, J., Joanny, J.F., Zinn-Justin, J., Elsevier, New York, (1989). [16] Kirkwood, J.G., Monroe, E., J. Chem. Phys. 9, 514 (1941). [17] Haymet, A.D.J., Oxtoby, D.W., J. Chem. Phys. 84, 1769 (1986). [18] Ramakrishnan, T.V., Yussouff, M., Phys. Rev. B19, 2775 (1979). 199  Bibliography  200  [19] Haymet, A.D.J., Oxtoby, D.W., J. Chem. Phys. 74, 2559 (1981). [20] Haymet, A.D.J., J. Chem. Phys. 78, 4641 (1983). [21] Laird, B.B., McCoy, J.D., Haymet, A.D.J., J. Chem. Phys. 87, 5449 (1987). [22] Jones, G.L., Mohanty, U., Molec. Phys. 54, 1241 (1985). [23] Hoover, W.G., Young, D.A., Grover, R., J. Chem. Phys. 56_, 2207 (1972). [24] Barrat, J.L., Hansen, J.P., Pastore, G., Waisman, E.M., J. Chem. Phys. 86, 6360 (1987). [25] Laird, B.B., Kroll, D.M., Phys. Rev. A 42, 4810 (1990). [26] Barrat, J.L., Hansen, J.P., Pastore, G., Molec. Phys. 63, 747 (1988); Barrat, J.L., Europhys. Lett. 3, 523 (1987). [27] Curtin, W.A. J. Chem. Phys. 88, 7050 (1988). [28] Rick, S.W., Haymet, A.D.J., J. Phys. Chem. 94, 5212 (1990). [29] Barrat, J.L., Baus, M., Hansen, J.P., J. Phys. C20, 1413 (1987); Zeng, X.C., Oxtoby, D.W., J. Chem. Phys. 93, 4357 (1990); Denton, A.R., Ashcroft, N.W., Phys. Rev. A 42, 7312 (1990). [30] Brami B.R., Joly F., Barrat, J.L., Hansen, J.P., Phys. Lett. A 132, 187 (1988). [31] Smithline, S.J., Rick, S.W., Haymet, A.D.J., J. Chem. Phys. 88, 2004 (1988). [32] Singh, U.P., Mohanty, U., Singh, Y., Physica A 158, 817 (1989); Marko, J.F., Phys. Rev. Lett. 60, 325 (1988); McCoy, J.D., Singer, S.J., Chandler, D., J. Chem. Phys. SI, 4853 (1987). [33] Ding, K., Chandler, D., Smithline, S.J., Haymet, A.D.J., Phys. Rev. Lett. 59, 1698 (1987). [34] Chandler, D., McCoy, J.D., Singer, S.J., J. Chem. Phys. 85, 5971 (1986). [35] Lovett, R., J. Chem. Phys. 88, 7739, (1988). [36] Lovett, R., Stillenger, F.H., J. Chem. Phys. 94, 7353 (1991). [37] Lutsko, J.F., Baus, M., Phys. Rev. A41, 6647 (1990). [38] Rosenfeld, Y., Phys. Rev. A 43, 5424 (1991).  Bibliography  201  [39] Verlet, L., Phys. Rev. 184, 150 (1969). [40] Ramakrishnan, T.V., Materials Science Forum 3, 57 (1985); Yussouff M., Phys. Rev. £23_, 5871 (1981); March, N.H., Tosi, M.P., Phys. Chem. Liq. 11_, 79, (1981); Mahato, M.C., Lakshimi, M.R., Pandit, R., Krishnamurthy, H.R., Phys. Rev. A 38, 1049 (1988). [41] Ramakrishnan, T.V., J. of Non-Crystalline Solids 117/118, 852 (1990); Sengupta, S., Sood, A.K., Phys. Rev. A 44, 1233 (1991). [42] Tisza, L., Generalized Thermodynamics M.I.T. Press, Cambridge, Mass., 1966. [43] Casas-Vazquez, J., Stability of Thermodynamic Systems eds Casas-Vazquez, J., Lebon, G., Springer-Verlag, 1982. [44] Prigogine, I., Defay, R., Chemical Thermodynamics, ch. 15-16, University Press, Glasgow, 1954. [45] Rowlinson, J.S., Swinton, F.L., Liquids and Liquid Mixtures, 3rd ed., ch. 6, Butterworth, UK, 1982. [46] Malescio, G., J. Chem. Phys. 95, 1198, 1202 (1991); Caccamo, C , Giunta, G., Hoheisel, C , Phys. Lett. A 158, 325 (1991). [47] Hansen, J.P., McDonald, I.R., Theory of Simple Liquids, 2nd ed., Academic Press, London, 1986. [48] Klein, M., Green, M.S., J. Chem. Phys. 39, 1367 (1963); Foiles, S.M., Ashcroft, N.W., Phys. Rev. A 24, 424 (1981); i Timoneda, J.J., Haymet, A.D.J., Phys. Rev. A 40, 5979 (1989). [49] Belloni, L., Phys. Rev. Lett. 57, 2026 (1986); Caccamo, C , Malescio, G., Phys. Rev. A 40, 6384 (1989); Malescio, G., Phys. Rev. A 42, 2211 (1990); [50] Abramo, M.C., Caccamo, C , Giunta, G., Phys. Rev. A 34, 3279 (1986); Caccamo, C , Malescio, G., J. Chem. Phys. 90, 1091 (1989); Caccamo, C , J. Chem. Phys. 91, 4902 (1989); Caccamo, C , Giacoppo, A., Phys. Rev. A 42, 6285 (1990). [51] Kirkwood, J.G., Buff, F.P., J. Chem. Phys. 19, 774 (1951). [52] Kerins, J., Scriven, L.E., Davis, H.T., Adv. Chem. Phys. 65, 215 (1986); Poll, P.D., Ashcroft, N.W., Phys. Rev. A 35, 5167 (1987). [53] Stillinger, F.H., Buff, F.P., J. Chem. Phys. 37, 1 (1962). [54] Lebowitz, J.L., Percus, J.K., J. Math. Phys. 4, 116 (1963).  202  Bibliography  [55] Ursenbach, C.P., Wei, D.Q., Patey, G.N., J. Chem. Phys. 94, 6782 (1991). [56] Verlet, L., Weis, J.J., Phys. Rev. A 5, 939 (1972). [57] Stell, G., The Equilibrium Theory of Classical Fluids, eds Frisch, H.L., Lebowitz, J.L., Benjamin, New York, 1964. [58] (a) Blum, L., Torruella, A.J., J. Chem. Phys. 56, 303 (1972); (b) Blum, L., ibid. £7, 1862 (1972); 58, 3295 (1973). [59] Messiah, A., QuantumMechanics,  vol.11, Benjamin, New York, 1973.  [60] Fries, P.H., Patey, G.N., J. Chem. Phys. 82, 429 (1985). [61] Kusalik, P.G., Patey, G.N., Molec. Phys. 65, 1105 (1988). [62] Kusalik, P.G., Patey, G.N., / . Chem. Phys. £9, 5843 (1988). [63] Carnie, S.L., Patey, G.N., Molec. Phys. 47, 1129 (1982). [64] Dyke, T.R., Muenter, J.S., J. Chem. Phys. 59, 3125 (1973). [65] Verhoeven, J., Dymanus, A., J. Chem. Phys. 52, 3222 (1970). [66] Morris D.F.C., Struct. Bonding 4, 63 (1968). [67] McMillan, W.G., Mayer J.E., J. Chem. Phys. 13, 276 (1945). [68] Friedman, H.L., Dale W.D.T., Modern Theoretical Chemistry, Vol. 5, Chap. 3. ed. Berne, B.J., Plenum, New York, 1977, [69] Ramanathan, P.S., Friedman H.L., J. Chem. Phys. 54, 1086 (1971). [70] Krishnan, C.V., Friedman H.L., J. Solution Chem. 3, 727 (1974). [71] Pettitt, B.M., Rossky, P.J., J. Chem. Phys. 84, 5836 (1986). [72] Hirata, F., Rossky, P.J., Chem. Phys. Lett. 83, 329 (1981); Hirata, F., Rossky, P.J., Pettitt, B.M., J. Chem. Phys. 78, 4133 (1983). [73] Pauling L., The Nature of the Chemical Bond. Cornell U.P., Ithaca, New York, 1960. [74] Ursenbach, C.P., Patey, G.N., J. Chem. Phys. 95, 485 (1991). [75] Gray, C.G., Gubbins, K.E., Theory of Molecular Liquids, vol. 1, Clarendon Press, Oxford, 1984.  203  Bibliography  [76] Gradshteyn, I.S., Ryzhik, I.M., Tables of Integrals, Series, and Products, Academic Press, New York, 1980. [77] Arfken, G., Mathematical Methods for Physicists Academic Press, Florida, 1985. [78] de Leeuw, S.W., Perram, J.W., Smith, E.R., Roy. Soc. London, (1980).  Proc. A 373, 27  [79] Kittel, C., Introduction to Solid State Physics 3rd ed., Wiley, New York, 1966. [80] Levesque, D., Weis, J.J., Patey, G.N., Molec. Phys. 51, 333 (1984). [81] McMullen, W.E., Oxtoby, D.W., J. Chem. Phys. 86, 4146 (1987); ibid, 88, 1476 (1988). [82] Jacucci, G., Klein, M.L., McDonald, I.R., J. Physique Lett. 36, L-97 (1975). [83] Tosi, M.P., Fumi, F.G., J. Phys. Chem. Solids 25, 45 (1964). [84] (a) Baus, M., J. Phys.: Condens. Matter 2, 2111 (1990) (b) Baus, M., Colot, J.L., Molec. Phys. 55, 653 (1985). [85] Franks, F., Water - A Comprehensive Treatise, vol. 1, ed. Franks, F., Plenum Press, New York, 1972. [86] Hobbs, P.V., Ice Physics Oxford University Press, England, 1974. [87] Burns, G., Glazer, A.M., Space Groups for Solid State Scientists, 2nd Ed., Academic Press, San Diego, 1990. [88] Pauling, L., J. Am. Chem. Soc. 57, 2680 (1935). [89] Hoover, W.G., Ree, F.H., J. Chem. Phys. 49, 3609 (1968). [90] Lindemann, F.A., Z. Phys. U , 609 (1910). [91] Eriksson, A., Hermansson, K., Acta Cryst. £ 3 9 , 703 (1983). [92] Leadbetter, A.J., Roy. Soc. London, Proc. A 287, 403 (1965). [93] Klewe, B., Pederson, B., Acta Cryst. £ 3 0 , 2363 (1974). [94] Ford, T.A., Falk, M., J. Mol. Struct. 3, 445 (1969). [95] West, C D . , Z. Kristallogr. 88, 198 (1934). [96] Brink, G., Falk, M., Can. J. Chem. 49, 347 (1971).  Bibliography  [97] Onsager, L., J. Phys. Chem. 43, 189 (1939). [98] Coulson, C.A., Eisenberg, D., Roy. Soc. London, Proc. A 291, 445, 454 (1966). [99] Levesque, D., Weis, J.J., Patey, G.N., J. Chem. Phys. 72, 1887 (1980). [100] Lutsko, J.F., Baus, M., J. Phys. C: Condens. Matter 3_, 6547 (1991).  204  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items