A T H E O R E T I C A L S T U D Y O F A M M O N I A - S A L T M I X T U R E S I N B U L K S O L U T I O N S A N D I N T E R F A C I A L R E G I O N S By John S. Perkyns B. A. (Chemistry), B. Sc. (Mathematics) Dalhousie University M. Sc. (Chemistry) University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF CHEMISTRY We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA March 1990 © John S. Perkyns, 1990 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. Department of CHff M I fTT The University of British Columbia Vancouver, Canada Date MfrY 2*1 MO DE-6 (2/88) Abstract Five models, ranging from a full molecular polar/polarizable model with C 3 v symme-try, to the primitive model of ions in a dielectric continuum, have been used to study the properties of ammonia both as a pure liquid and as a solvent. Ammonia is modelled as a multipolar polaxizable hard sphere and ions as charged hard spheres. Using these models, in conduction with the Reference Hypernetted-chain integral equation theory, ammonia has been studied as a pure liquid, as a solvent near charged and uncharged surfaces, and in electrolyte solutions of finite concentration. The formalism of Kirkwood and Buff was used to obtain thermodynamic quantities of ionic solutions from calculated distribution functions. Structural, thermodynamic and. dielectric properties were calculated for pure am-monia and were compared with experiments where possible. Values for the dielectric constant and the configurational energy were found to compare favorably with experi-ment. Ammonia formed a relatively dense, highly structured layer within two solvent diameters of an uncharged surface. This structure was analyzed in terms of angular probability distributions of the molecular dipole vector and the NH-bond direction, and was found to be similar to that of frozen ammonia. The extreme asymmetry of solvation of unlike charges in ammonia was also investigated. Small cations were found to be more favorably solvated than small anions, but as the ion size was increased, the situation reversed. Estimates of the number of ion pairs in liquid ammonia and their effects on the be-havior of mean molar activity coefficients were examined. Large differences between ex-perimental activity coefficients and the Debye-Hiickel hmiting law could not be explained n by the usual ideas of ion pairing. It was found that the integral equation theories appear to have no solution between ionic concentrations of about 2 x 1 0 _ 4 M and 2 x 1 0 _ 2 M for either molecular or continuum models. Using rigorous stability criteria, this behavior was shown to be consistent with the onset of a phase change. It is proposed that such phase separation phenomena might explain the unusual behavior of the experimental activity coefficients. in Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement xi 1 Introduction 1 2 Theory 8 2.1 Statistical Mechanical Theory. 8 2.1.1 Considerations in Choosing Models. 8 2.1.2 Model Definitions 10 2.1.3 The Reference Hypernetted-chain (RHNC) Theory 17 2.1.4 Reduction of the OZ Equation 23 2.1.5 Reduction of the Closure 29 2.1.6 The Static Dielectric Constant and Average Quantities 32 2.1.7 The Self-Consistent Mean Field Treatment of Polarizability Effects. 36 2.1.8 Angular Probability Distributions 43 2.2 Theory for Ionic Solutions 46 2.2.1 Kirkwood-Buff Theory for Electrolyte Solutions 47 2.2.2 Stability Criteria for Ionic Solutions 51 2.2.3 Limiting Behavior 53 iv 2.2.4 Individual Ion Partial Molar Volumes 56 2.3 Ideas of Ion Pairing 59 2.3.1 The Bjerrum Theory. . 59 2.3.2 Other Definitions 61 2.3.3 Ion Pair Definitions for the Present Study 62 2.4 Computing Considerations 64 2.4.1 Method of Numerical Solution 64 2.4.2 Choice of Basis Set . 64 2.4.3 Handling Long-Ranged Potential Tails 66 3 Results and Discussion 77 3.1 Results for Models of Pure Ammonia 77 3.2 Results for Macrospheres at Infinite Dilution. . 96 3.2.1 Reference Hard-Sphere Bridge Functions for Large Particles. . . . 96 3.2.2 Model Dependence in Large Particle calculations 120 3.2.3 Dependence of Radial Distribution on Macrosphere Size 134 3.2.4 Dependence of Angular Structure on Macrosphere Size 137 3.2.5 Combined Positional and Orientational Distributions. . . . . . . . 149 3.3 Results for Ions at Infinite Dilution 165 3.3.1 Dependence of Solvation on Ion Size. 165 3.3.2 Solvent Distribution About Ions 185 3.4 Dependence of Ionic Behavior on Concentration 204 3.4.1 Examination of Experimental Results 207 3.4.2 Results for Electrolytes in Liquid Ammonia 214 3.4.3 Indication of a Molecular Model Phase Change 226 4 Summary and Conclusions 242 v Bibliography 249 vi List of Tables 2.1 Function index restrictions due to particle symmetry 13 2.2 Non-zero multipole elements for each model 15 2.3 Computing resources needed for each basis set 65 3.4 Input parameters for polarizable pure ammonia models 77 3.5 SCMF effective dipoles and renormalized polarizabilities 78 3.6 Thermodynamic quantities, dielectric constants and coordination numbers for pure ammonia . . . 79 3.7 Proportion of Free Ions a Calculated From Model III. 218 vii List of Figures 2.1 The molecular axis system for the ammonia molecule 18 2.2 Comparison of absolute errors for a test function 70 2.3 Result of logarithmic transform of a step function into t-space 73 3.4 Pure ammonia radial distribution functions 80 3.5 ^/x1 • fi^j vs. r for pure ammonia 83 3.6 hH°(r) for pure ammonia 85 3.7 h-QQ2(r) for pure ammonia 88 3.8 ^oo3(r) f° r P u r e ammonia. . 90 3.9 hH4(r) for pure ammonia 92 3.10 hll6(r) for pure ammonia 94 3.11 Hard sphere (reference) radial distribution functions 97 3.12 Hard sphere bridge functions 101 3.13 Model I 9ooms(r) functions using different bji(r) approximations 104 3.14 Model I / i o"m«( r ) functions using different bR(r) approximations 106 3.15 Model I h°l^ms(r) functions using different bn(r) approximations 108 3.16 Model I h°ofma(r) functions using different 6i?(r) approximations 110 3.17 Probability density for the dipole-surface normal angle at contact. : . 113 3.18 Probability density for the dipole surface normal angle at r = 0.bds. . . . 115 3.19 Probability density for the dipole surface normal angle at r = 1.0d„. . . . 117 3.20 Model dependence of g^ms{r) 121 3.21 Model dependence of C U ( r ) 123 vm 3.22 Model dependence of h°™ms(r) 125 3.23 Model dependence of K™me(r) 127 3.24 Model dependence of PfZ(6) at contact 129 3.25 Model dependence of PNH{@) at contact 132 3.26 Dependence of 5fSo^n«(r) o n hard sphere size 135 3.27 Dependence of ^ocwnst7") o n n a r c l sphere size 138 3.28 p^(d,r) for lda diameter hard spheres in Model I N H 3 140 3.29 Pti(0,r) for ZdB diameter hard spheres in Model I N H 3 143 3.30 Pn(8,r) for 5ds diameter hard spheres in Model I N H 3 145 3.31 Pv(8,r) for 20ds diameter macrospheres in Model I N H 3 147 3.32 PNH{8,T) for \ds diameter hard spheres in Model I N H 3 150 3.33 PNH(8,T) for 20<is diameter macrospheres in Model I N H 3 152 3.34 The dominant orientation of ammonia at a macrosphere surface. . . . . . 154 3.35 goo%is{r)P^^r) f° r 20<fs diameter macrospheres in Model I N H 3 . . . . . . 157 3.36 gm?ms(r)Ptl(8,r) for Id, diameter hard spheres in Model I N H 3 . . . 159 3.37 9ocrmii(r)PNH{@,r) f° r 20ds diameter macrospheres in Model I N H 3 . . . . . 161 3.38 gm°ms{r)PNH{6,r) for lda diameter hard spheres in Model I N H 3 163 3.39 g^is(die) for hard-sphere ions in Model I N H 3 . 166 3.40 UjsT/NkT for hard-sphere ions in Model I N H 3 169 3.41 UID/NkT for hard-sphere ions in Model I N H 3 171 3.42 Uio/NkT for hard-sphere ions in Model I N H 3 174 3.43 Uw/NkT and UIH/NkT for hard-sphere ions in Model I NH 3 . 176 3.44 Coordination numbers for Model I N H 3 molecules around hard-sphere ions. 178 3.45 Coordination sphere density for Model I N H 3 molecules around hard-sphere ions 181 3.46 Partial molar volumes at infinite dilution 183 ix 3.47 g$%(r) for ions in Model I N H 3 186 3.48 h%£it(r) for ions in Model I N H 3 188 3.49 h™s(r) for ions in Model I N H 3 190 3.50 h°™is(r) for ions in Model I N H 3 192 3.51 Pfi(0,r) for Model I NH 3 molecules near an anion : 195 3.52 PJVK(^) 7 ') ^ o r Model I N H 3 molecules near an anion 197 3.53 Pn(9,r) for Model I N H 3 molecules near a cation 199 3.54 PNH(8,r) for Model I N H 3 molecules near a cation 201 3.55 Primitive Model Spinodal Decomposition Line Predicted by the HNC The-ory 205 3.56 Experimental mean molar activity coefficients for 1:1 Electrolytes. . . . . 208 3.57 Experimental mean molar activity coefficients for 1:2 Electrolytes. . . . . 210 3.58 G+- vs. p2 for 1.52<is diameter ions in NH 3 . . 216 3.59 vp2G^ vs. p2 for 1.52<is diameter ions in NH 3 . • 220 3.60 y/p^(d\a.y±ldp2)Ttp vs. ^J~p~2 for 1.52 ,^ diameter ions in N H 3 . . . . . . . . 222 3.61 sfp-tipIny±/dp2)x,P v s - yfpi f° r l-52<is diameter ions in N H 3 224 3.62 lny± vs. ^fp2 for 1.52d„ diameter ions in N H 3 . 227 3.63 l / p 2 G + _ vs. concentration for Models III, I V and V 230 3.64 l/p2G+- vs. concentration for Models I V and V 233 3.65 Mechanical and Compositional Stability Criteria for Model III 235 3.66 lny± vs. ^/p~2 for 1.52d„ assuming a phase change 239 x Acknowledgement I wish to thank my academic advisor, Dr. G. N. Patey, for his guidance and encour-agement over the last six years, and in particular the time he has expended on my behalf over the last few months. I would also like to thank the University of British Columbia Chemistry Department and the Natural Sciences and Engineering Research Council of Canada (NSERC) for their financial support, Dr. P. G. Kusalik and Dr. G. M. Torrie for the use of some computer programs, and Dr. J. P. Valleau, Dr. P. Attard and C. Ursenbach for the use of unpublished results. I would also like to thank my wife, Jane, who managed to be supportive while com-pleting a doctoral degree of her own. xi Chapter 1 Introduction There is a vast wealth of information in the literature, both experimental and theoret-ical, dealing with many aspects of water both as a pure liquid and as a solvent, reflecting the importance of water to human endeavor. By comparison, non-aqueous solvents, and ammonia in particular, have not been as intensely studied and are therefore perhaps less well understood. This lack of understanding comes in large part because water is much more the exception than the rule; it has physical properties unlike most other solvents. Yet, theories of solvation and ionic interaction, such as the Debye-Huckel (DH) approach [1] and Onsager's conductivity theory [2], while applying to electrolyte solutions in all solvents, were developed primarily for the purpose of understanding aqueous systems. There has been much success in improving the basic ideas of Debye and Hiickel to obtain theories which accurately describe water solutions [3], so attempts to understand the behavior of ammonia have also grown out of these ideas, with only modest success. Thus there are two problems for any theoretical study of liquid ammonia. There is relatively little experimental or theoretical work on which to build, and the theories and methods which are needed to understand such solutions were developed largely for systems which have markedly different physical properties. Experimentally determined mean molal activity coefficients, 7-t, for electrolytes in non-aqueous solvents usually deviate significantly from the Debye-Hiickel limiting law at ionic concentrations well below those where differences can be seen for corresponding aqueous systems [3,4]. The activity coefficients measured for electrolytes in ammonia 1 Chapter 1. Introduction 2 deviate in the opposite direction to those for most other electrolyte solutions, and they do not appear to follow the Debye-Htickel limiting law anywhere [5-8]. Despite this difference, attempts to explain the behavior of activity coefficients in ammonia have been centered around the idea of ion pairing, even though this contention is not supported by Raman spectroscopy [9]. The concept of a pair of ions held together only by electrostatic forces was first treated qualitatively by Bjerrum [10,11] who, by building on the equations of Debye and Hiickel, developed a method for calculating proportions of "free" (a) and "paired" (1 — a) ions. The newly defined species, the ion pair, is seen to be in equilibrium with free ions and so, in the usual manner, we can write an expression for the related association constant in terms of electrolyte concentration c, 1 — a , The Bjerrum theory is based on a somewhat arbitrary definition of an ion pair: Ions at separations such that their mutual potential energy is greater than twice the average thermal energy of a solvent molecule are considered paired. Fuoss objected to Bjerrum's definition of an ion pair and developed two theories of his own [12,13] in an attempt to remove some of the arbitrariness. For any theory of ion pairing, however, some definition of an ion pair must be used. Since there is no distinct potential barrier, such as that of a covalent bond, preventing the separation of electrostatically attracted ions, some measure of arbitrariness is inevitable. There is also a more serious objection to the use of ion pairing notions to explain the experimental activity coefficients for ions in ammonia. Unfortunately, any simple method of accounting for deviations from the DH limiting law by including ion pairing ideas predicts that the activity coefficients will deviate in the opposite direction from that which is actually found experimentally. Another property of interest in the study of electrolyte-ammonia solutions is the Chapter 1. Introduction 3 asymmetric solvation of positive and negative charges. This is in turn due to the partic-ular symmetry of the charge distribution around a single ammonia molecule. The three hydrogen centers spread the local positive charge over a wide region around the molecule, whereas the single lone pair of electrons create a "negative end" where the negative lo-caL charge is more focussed. Consequently the interaction between a positive ion and an ammonia molecule is characterized by a different total angular dependent potential than is the corresponding interaction between a negative ion and an ammonia molecule. The effects of this highly asymmetric solvation on thermodynamic quantities have not been examined for ammonia, and neither has a possible connection between the activity coefficients and solvation properties been investigated. One last subject of interest discussed in the literature concerns the so-called "primitive model" [14-16], or ions immersed in a dielectric continuum. There has been much interest in recent years in an apparent phase change predicted by this model [17-23] using many approximate theories. Unfortunately, the phase changes predicted by these theories give critical temperatures which differ by a factor of two and critical densities which differ by a factor of more than twenty. Recent Monte Carlo calculations [24,25] place a lower bound on the critical density at a somewhat higher value than the approximate theories, but accurate values have not yet been determined. Since it is not known where the correct critical point for the primitive model lies, it is not easy to judge the accuracy of any of the theories. All theories agree, however, that the primitive model predicts phase transitions at relatively low densities. How well the primitive model approximates real electrolyte-ammonia solutions is also not known. It has been known since the beginning of this century [26,27] that dilute solutions of alkali metals in ammonia separate into two immiscible liquid phases. About forty years ago, it was noticed that the addition of small amounts of salts had a marked effect on the temperature of phase separation [28,29]. Since then more data on phase separations in alkali metal/ammonia solutions, Chapter 1. Introduction 4 with and without electrolytes, have become available [30-33]. Such phase instabilities are not always visible and so are confirmed by indirect means. In light of these ideas, one motivation of this study is to examine the concepts of ion pairing and asymmetric solvation, with a view to explaining the large deviations from the Debye-Hiickel limiting law. Because ammonia has not been studied extensively, another motivation is to perform a relatively comprehensive survey of ammonia in situations where it has received little attention. In order to make these goals tenable, we will focus on the structural and thermodynamic properties of pure ammonia and the effects of ion size on solvation in ammonia. We will also examine the structure of ammonia near a macrosphere or macroion surface. Finally, we will investigate the possibility of a phase change in low concentration electrolyte-ammonia solutions, and examine what qualitative effect a phase change would have on the mean molar activity coefficient. Five models covering a wide range of complexity will be used here. The first three models define ammonia as a molecular species, characterized by a hard-sphere core with a set of electric point multipoles embedded at the centers. Each model represents a different truncation in the multipole series. Model I includes all multipoles up to hexade-capole, Model II up to octupole and for Model III the multipole series is truncated at the quadrupole. Since the quadrupole of ammonia is axially symmetric, Model III ammonia also has this symmetry. Models I and II, however, conform to the true C3„ symmetry of the ammonia molecule. The asymmetric solvation of oppositely charged ions in ammo-nia will be apparent in Models I, II and III. Unlike water [34], the combination of the dipole and quadrupole in ammonia are sufficient to induce a high degree of asymmetric solvation. The solvent is also treated as polarizable using the Self-Consistent Mean Field (SCMF) approach [35]. The polar-polarizable solvent is reduced to an effective system characterized by an effective permanent molecular dipole. Molecular models of this type have been used with some success for water-like solvents [34, 36-38], and a model very Chapter 1. Introduction 5 similar to Model III has recently been applied to pure ammonia [39]. The results obtained were in good agreement with experimentally determined radial distribution functions [40] and dielectric properties [41], as well as the same quantities determined by computer sim-ulations for the same model. Computer simulations can be regarded as essentially exact for a given model, and as such provide a test of the approximations used in the integral equation theories applied here. McMillan and Mayer [42] showed that the solvent in an ionic solution need not be considered as a molecular species in the low concentration limit if solvent effects are included as a short ranged effective ion-ion potential. Model IV considers just the ions as distinct species, with the dependence of ions on solvent behavior included in this manner. The potential for a pair of ions which defines Model IV is simply the infinite dilution potential of mean force Wij(r) for a molecular solvent. For Model IV the total ion-ion potential can be written in the form Uij(r) = w{j(r) = < ( r ) + ^ (1-2) where w°j(r) is the short ranged effective ion-ion interaction due to the molecular nature of the solvent, r is the distance between the ions i and j with charges qi and qj, and e is the dielectric constant of the pure solvent. The potential of mean force Wij(r) which is used here is taken from a calculation for ions at infinite dilution in Model III ammonia. In the infinite dilution limit, Models III and IV become equivalent. Approaches of the McMillan-Mayer (MM) type have been widely used in the literature [16, 43-45]. Model V, the primitive model, consists simply of hard-sphere ions in a dielectric continuum. It is equivalent to Model IV with w*j(r) set to zero. Integral equation theories have been used extensively in conjunction with primitive model calculations [14-16], molecular multipolar fluids [46-48], and aqueous electrolytes Chapter 1. Introduction 6 [34, 36-38]. Calculations using these theories incorporate a set of experimentally mea-sured physical properties (electric multipoles, polarizabilities, densities, effective diame-ters) and using no freely adjustable parameters derive probability distribution functions for pairs of components. Integral equation methods are based on the equation of Orn-stein and Zernike [49] which relates the pair distribution function, 5(12), to the direct pair correlation function c(12). The direct correlation function was originally expected to depend only on the pair potential, but now it is known to depend on other functions in a manner which as yet has no simple exact expression. Each integral equation theory is characterized by a further relationship between g(12) and c(12) known as a closure, and the approximations inherent in each theory are due to closures being inexact. Closures such as the Mean Spherical Approximation (MSA) [50], the Percus-Yevick (PY) theory [51] and further approximations to the Hypernetted-chain (HNC) theory [52,53] in which the logarithmic term was expanded and truncated have been readily solvable for angu-lar dependent potentials for some time, but a method of solving the full HNC equation has been found only recently [48]. The HNC closure together with the Ornstein-Zernike (OZ) equation is known as the HNC integral equation theory. We will be using the Ref-erence Hypernetted-chain (RHNC) theory for most calculations. The RHNC theory uses a perturbation technique introduced by Lado [54,55], which writes all the pair functions in terms of deviations from a reference system which is known exactly. The reference systems we will be using are simple hard sphere fluids. The RHNC theory is now being used extensively [36,37,46,48,56]. For Models I, II and III the RHNC theory is applied and for Models IV and V the HNC theory is used. In the limit of low density and small interaction potential, the solution to the HNC theory approaches an exact solution. The accuracy of the HNC theory has been questioned for interaction potentials and densities similar to those used for Model V in this study [18,57] so we will examine such results very carefully. Chapter 1. Introduction 7 This thesis is divided into four Chapters. Chapter II is devoted to theory, and we begin by defining the models that will be used. The mathematical machinery of the integral equation methods is then presented followed by the SCMF theory and the method for calculating thermodynamic quantities of interest. Some results from the Kirkwood-Buff theory [58] for electrolyte solutions are then given, as well as a presentation of some ideas of ion pairing. A section on computing considerations follows, giving an estimate of the computing resources used, and a discussion on the method of logarithmic fourier transforms. Chapter III contains a discussion of the results obtained, and begins with a section on pure ammonia. The macrosphere/macroion interface is then examined and finally a section comparing ionic solutions of finite concentration to some experimental results. Chapter IV summarizes the main conclusions. Chapter 2 Theory 2.1 Statistical Mechanical Theory. 2.1.1 Considerations in Choosing Models. This study will not primarily be concerned with theoretically duplicating real sys-tems, partly because the enormous complexity of real systems requires correspondingly complex models, and partly because there is relatively little experimental data for com-parison. Instead, we wish to isolate a few particular properties for investigation. For a model to be useful it must unambiguously display those phenomena of a physical system which will answer the questions that originally motivated the study. The first question to be considered here is, which properties of a single ammonia molecule are responsible for the properties of bulk liquid ammonia and ammonia at a surface, and in particular, what are the causes and consequences of asymmetric solvation? This is essentially a restatement of the fundamental goal of statistical mechanics, to form a bridge between microscopic and macroscopic properties of matter. The second question to be considered is, what is the nature of the phenomena observed for ionic solutions with low dielectric constant solvents? The phenomenon of particular interest here is the enormous difference between the low concentration behavior of the mean molar activity coefficients and what is predicted by the Debye-Hiickel limiting law. Generally the simplest model which will display these desired phenomena will yield the most insight into the nature of the physical system itself. As mentioned in the 8 Chapter 2. Theory 9 Introduction, a sequence of five distinct models of decreasing complexity is used in this study with the aim of showing the minimum requirements for a given phenomenon to be present in a system. It is also our intention to show that when the model complexity is increased to where it more closely describes a real system, some thermodynamic quantities which result from the model, such as configurational energy and dielectric constant, will approach reality. We choose models with essentially the same features as those which have been previously used for both ammonia and water with some success [34,37,38,43,59], in an attempt to show which model parameters lead to which experimentally observed phenomena. The bulk properties of liquids are determined from interparticle potentials using sta-tistical mechanical methods. These interparticle potentials describe the interactions be-tween molecules and are electrostatic in nature, arising from the charge distributions of molecular species. The potentials have two basic features, a short ranged repulsion caused by the overlap of electron clouds and nuclear repulsion, and longer ranged forces due to coulombic interaction between net asymmetries in positive and negative charge distributions. The charge distributions can be described simply by the multipole series, and the effect of long-range forces can be written as interactions between the multipoles. The models are chosen to obey several conditions. They are purely classical models, and the machinery of classical statistical mechanics is used. Quantum effects are expected to be small over the range of temperatures used [60]. The molecules are assumed to be rigid, with the effects of intramolecular vibrations reflected in an average way by the multipole moments and the experimentally determined polarizability. It is also assumed that contributions to the potential between particles due to classical polarizability effects can be written in terms of effective two-body interactions only. These conditions could be considered as approximations to a real system introduced for mathematical and com-putational simplicity, but they do not compromise the quahtative validity of the results Chapter 2. Theory 10 if the models used display the phenomena being investigated. They are quite distinct from the approximations used in the theory itself, which will be discussed below. 2.1.2 Model Definitions introduce the mathematical language in which they and other functions will be expressed. The space of angles for two particles in an arbitrary laboratory frame is exactly spanned by an orthogonal basis set of functions known as rotational invariants, where R™^{Q) is a Wigner generalized spherical harmonic [61], fi,- — (xi,6i,(f>i) are the Euler angles [60,61] for the orientation of particle z, (:::) is the usual 3-j symbol [62], and m,n,l,p,i/,p',v' and A' are integers taking the values for which the 3-j symbols are non-zero, and for which the Wigner generalized spherical harmonics are defined. The notation r represents the polar and azimuthal angles of the vector from particle 1 to particle 2, whose magnitude, the distance between particle centers, will be denoted r, The relations between the indices necessar3r for non-zero 3-j symbols are In order to give expressions for the two body potentials, w(12), it is necessary to (2.3) m — n\ < I < m + ra, (2.4a) p\ < m and \p'\ < m, (2.4b) 0 v\ < n and < n. (2.4c) The coefficients fm are arbitrary non-zero constants, and for mathematical conve-nience the definition fmnl = [(2m + l)(2n + l ) p , (2.5a) Chapter 2. Theory 11 will be used to simplify equations, and is the form which is actually used to solve the equations. The alternate form rmnl m n I v 0 0 0 j will be used when the r-dependent projections (see definition below) are given, in order to be consistent with the literature. Since the pair potentials and all other pair distribution functions used in this study depend only on the relative particle orientations and the distance between their centers, they can be written as expansions in rotational invariants, in the form fap(12) = ZfZit(r)*^in1,&,£), ' (2-6) mnl where the coefficient functions, f™ap{r), often called projections, are dependent only on r, and the notation implies that particle 1 is of species a and particle 2 is of species (3. The pair potential uap(12) can be considered as a special case of the general pair function /<rf(12). All the pair functions must obey the following properties: (i) Invariance on rotation of the lab frame. (ii) Invariance on translation of the lab frame. (iii) Invariance on interchange of identical particles. (iv) The functions must be real. (v) Invariance under the symmetry operations of particles 1 and 2. (vi) Invariance under positive and negative charge interchange. This choice of basis is a convenient one because it makes it possible to enforce these properties easily. Properties (i) and (ii) are true for homogeneous isotropic fluids such Chapter 2. Theory 12 as those considere here, and have been assumed in order to write the present expansion. Property (iii) can be expressed as A ( r ) = (-r+B/^L(')- (2-7) Property (iv), together with the general principle that the product of two real quantities must be real,, generates the condition, f™Lp(r) = {-T+M"+v f-;l-^{r). (2.8) Equations 2.7 and 2.8 have the effect of decreasing the number of unique projections. An-other condition, which applies to the potential only, arises from the multipole expansion (equation 2.13, below), and can be written as m + n = I. (2.9) The most restrictive conditions are those due to the symmetries of the pair of species described by a given pair function (property (v)). For the combinations of symmetries used in this study, the restrictions on the function indices are given in Table 2.1. For all symmetries used here, the symmetry restrictions in Table 2.1 in conjunction with equation 2.8 also imply the condition, m + n + I — even, (2-10) for all functions. There are in general, for all distribution functions, an infinite number of projections contributing to each expansion. It has been found, however, that these expansions converge relatively quickly [47,59,63,64], and so for all calculations in this study the restrictions m < 4, n < 4 and I < 8 are used. In previous studies with similar models the thermodynamic quantities of interest were found to have converged to within less than 1% when these restrictions were used. Chapter 2. Theory 13 Symmetry of species a Symmetry of species 8 spherical linear C31; spherical m, n,l,p,v =0 m,/x, 1/ = 0 n = I m,// = 0 n = / 1/ = 0 ,±3,±6,.. . /o"o^(r) = (~ Y fo-v;a.i3{r) linear n, p, v = 0 m = I fi,u = 0 not used n,u = 0 m = / / i = 0,±3,±6,... fjiO;a0 (r) — {~ Yfl-!iO;aP (r ) not used p,v = 0,'±3, ± 6 , . . . Table 2.1: Function index restrictions due to particle symmetry. Using this basis expansion for distribution functions in the equations which relate them is also useful because the relationships between individual invariants is well estab-lished. In particular, the orthogonality of the basis elements, expressed by [64], J$mi , n i 1 , 1 (fia, a,, t )$m 2 '"2 h (fij, a,, 1 )dn1dn2 = ( fmmili\2/o 2\2 .(2m 1- r-l)(2n 1 + l)(2Z1 + l ) 1 J < W 2 < W 2 ^ 2 , ^ . n j where 6 m i m 2 is Kronecker's delta, ensures that two invariant expansion expressions are equal if and only if the r-dependent coefficients are individually equal. This is a powerful tool in integral equations involving distribution functions, leading to relatively simple expressions for thermodynamic properties of interest which depend on only a very small Chapter 2. Theory 14 subset of projections. This in turn makes it possible to trace the elements of the pair potentials directly to their effects on the thermodynamics of the system. As stated above, the potentials have a short-ranged repulsive contribution. For all species in this study a simple hard-sphere core will be used to model this feature, which contributes to only one projection, u^0ai3(r), 0 for r > dan, ~ (2.12) oo for r < d^, where da@ = a —- and da,dp are the hard-sphere diameters of species a and j3. Long-ranged electrostatic interactions are modelled by a truncated series of point multipoles embedded at the centers of the hard spheres. The interaction of a pair of real molecules would have this limit at large separations, but for the present models we assume the whole charge distribution for all species to be inside the hard sphere so the multipole series will converge everywhere outside the sphere. The interaction potential of two non-overlapping charge distributions for spherical, linear or C3„ symmetrjr can conveniently be written as an invariant expansion with each term corresponding to a multipole moment in spherical tensor form [65,66] -Wi2) = E ( - r ( ( 2 m { 2 ^ ! + 1 ) , ) 5 c n ^ a , (2.13) with I = m + n in all terms and the superscript e denoting the electrostatic contribution. Here Qf^.a is defined for a discrete charge distribution by [65,66] 0£ ; Q = E [eiC^o"(a)] , (2-14) i where the discrete charges e; contributing to the charge distribution of species a are lo-cated at the point (r^ , fi;) in the molecular coordinate frame, and * indicates the complex conjugate. This quantity represents an element in the electric m-pole moment of species Chapter 2. Theory 15 a. The potential expansion will have fewer non-zero terms than the other distribution functions because the multipoles have few non-zero elements for the symmetries used here. Table 2.2 gives the non-zero multipole elements for each model. If we now simply Model Symmetry Multipoles included Multipole elements I ions spherical charge solvent dipole quadrupole octupole hexadecapole QI QI Q3 J Qzi Q3 QAZ,Q1,Q\ II ions spherical charge Qo solvent dipole quadrupole octupole Q°2 Qz >QziQ% III ions spherical charge QI solvent linear dipole quadrupole Q\ Qi IV ions spherical charge Q°o V ions spherical charge Q°o Table 2.2: Non-zero multipole elements for each model. combine the hard-sphere core and electrostatic parts of the potential we have the full expression for the potential, uafi(12)=uHS.^+ueal3(12), (2.15) Chapter 2. Theory 16 which with equations 2.12 and 2.13 and Table 2.2 defines Models I, II and III. Model IV is a McMillan-Mayer (MM) [42] level model, where only the ionic species are explicitly defined, and they are the same as the previously defined models. At the M M level the effect of the molecular solvent is included as a constant, short-ranged contribution to the potential of mean force, wsR;a.p(r), which is added to the ionic interaction to give, where a and j3 are ionic species only, and the potential of mean force is taken from a full molecular solvent calculation in the limit of infinite ion dilution. In this study we obtain wSR;a.p{r) from a Model III calculation. .The defining expression for wsR-a/3(r) is given by equation 2.138 in section 2.2.3 below. In the infinite dilution limit Models III and IV become equivalent. Model V is the primitive model, with the ions as in all the other models, but with the solvent reduced to a dielectric continuum. Model V can be seen as a.model IV calculation with wsR-api'r) set to zero. It should be noted that Models I through V are a series of simplifications in the model for the solvent, liquid ammonia, with the first three representing progressively severe truncations in the multipole series. This natural sequence of models will be used to elucidate the dependence of thermodynamic quantities on solvent properties. It should also be noted that the quadrupole tensor of ammonia obeys axial symmetry, and the well-known tripod structure of ammonia is not represented in the multipole series until the octupole is included. The multipole coefficients are used in spherical tensor form, but can be expressed in the more familiar cartesian tensor notation. Conversion between the two systems has been outlined in the literature [67,68] and for the symmetries used here the following relations hold: ua/3{12) = uHS]a/3 + uea0(12) + wSR.<ap(r), (2.16) Qo=q, (2.17a) Chapter 2. Theory 17 Q\=pz, (2.17b) Q°2 = ®zz, (2.17c) Q°3 = nzzz, (2.17d) Ql = -Os 3 = - (^ ) 2 (fix-x - 3 0 ™ ) , (2.17e) Ql = *zz*z, (2-17f) Q 3 = -Q4 3 = - ( ^ ) 2 - 3$ x y s / 2 ), (2.17g) where the usual notations oiq,p, Q,Q and $ represent charge, dipole, quadrupole, oc-tupole and hexadecapole moments, respectively. The cartesian molecular axis system corresponding to the definitions of the multipole moments given above is shown in Fig-ure 2.1. 2.1.3 The Reference Hypernetted-chain (RHNC) Theory In general the problem we wish to solve is the procuring of n-particle distribution functions, gffl(Xi, X 2 , X n ) where N is the number of particles in the system and the Xi denote the positions and orientations of the N particles. These distribution functions are related to the interparticle potentials uj^(X\, X2,Xj^) by [14], („) p~nN\ J exp[-8uN(X1, X 2 , X N)} dXn+1dXn+2...dXN gN ^uA2,...,A.n) (N_n)iSexp[_pUNfx1}X2,...,XN)]dX1dX2...dXN ' 1 j where 8 = (kT)-1, k is the Boltzmann constant, T the absolute temperature, and p = N/V is the number density, with V the volume of the system. If we are only concerned with two-body and effective two-body interactions, the problem becomes one of finding Chapter 2. Theory 18 Figure 2.1: The molecular axis system for the ammonia molecule. The molecular dipole is aligned with the z-axis and one hydrogen center is located in-the xz-plane. Chapter 2. Theory 20 <7(12) = g KN'(Xi,X2), which only explicitly depends on the positions and orientations of two particles so we write, / p~2N(N ~ 1)1 exp [Et<j -Bu(ij)\ dX3dX4...dXN g(l2) = = 1 / • 1 ]- . (2.19) /exp p i K j -Bu(ij)\ dX-i.dX2-.dXN The method used here to obtain ^(12) is one of a group of integral equation theories and is based on an equation first used in the study of density fluctuations near critical points known as the Ornstein-Zernike (OZ) equation [49], h(\2) = c(12) + ~ [ Ml3)c(32)<LY3, (2.20) where h{12) = g(12) - 1 (2.21) is called the total correlation function. The OZ equation in this form can be considered as a defining relation for c(12), the direct correlation function. It was at first anticipated that c(12) would depend on it(12) only, but it is now known [14] to depend on7i(12) as well. For the hard sphere potentials used in this study, it is numerically convenient to define a function 77(12) = h(12) — c(12) which, for the theories used here, is continuous everywhere, whereas both h(12) and c(12) are discontinuous. If we generalize equation 2.20 to represent a multicomponent system, and rewrite it in terms of 77(12), the OZ equation becomes 7^(12) = Y, J fo*r(13) + ca7(13)] c 7 / 3(32)dX 3 l (2.22) where a, 6 and 7 can represent any components in the system. The problem now becomes one of finding another independent relation between cQ|g(12) and hap(12). This is done by introducing another independent relation between cap{12) and the other previously defined functions. Such a relation is known as a closure. Rushbrooke and Scoins [69] Chapter 2. Theory 21 showed that cap(\2) can be written as a density expansion in multidimensional integrals known as clusters or diagrams, as can hap(12) and wap(12) = — In{gap(l2)} [70,71], from which exact functional relationships can be deduced. Using these relationships, and introducing the bridge function, bap(12), the exact closure can be expressed in the form c^(12) = M 1 2 ) -ln{<7«g(12)} - + ba0(12). (2-23) There is one class of diagrams in the cap(12) expansion, known as bridge diagrams, which cannot simply be expressed in terms of ca(g(12) or hap(12). While each individual bridge diagram integral could be calculated separately, the computing resources needed for only a few of them would be prohibitive, and many of these integrals give significant contributions. As is written in equation 2.23 the sum of these integrals can be considered as a single function, fc0/g(12). While fca/g(12) is not known analytically, it has been shown by computer simulations [14,72,73] to be short ranged for many systems, and since the most significant diagram in its expansion varies with p2, with all other diagrams varying with higher powers of density, it becomes insignificant in the limit of low density. If c0/g(12) is approximated by setting 6ajg(12) to zero one has the hypernetted-chain (HNC) equation, c ^ ( 1 2 ) = f e a / 3 ( 1 2 ) - l n { 5 ^ ( 1 2 ) } - ^ | ^ , (2.24) or if the system being studied can be regarded as a perturbation to a similar "reference" system for which fcajg(12).is known exactly, one can obtain the reference hypernetted-chain (RHNC) equation [48], cal3(12)= ha0{l2) - ln{ga/9.(12)} - U°f^ + bR.a^l2), (2.25) where the subscript R denotes the reference system. The RHNC equation is usually written [47,48,59,63,64] in terms of deviations of each of the functions from the reference system. This way of writing the RHNC equation is exactly equivalent, however, and Chapter 2. Theory 22 although the perturbation nature of the equation is not as apparent, it is a more conve-nient form for numerical calculation. This perturbation approach was first introduced by Lado [54,55], and it ensures that the solution is exact in the limit of ua/g(12) approaching UR;ap(l-2)- Equations 2.24 and 2.25, together with the OZ equation, are known as the HNC integral equation theory, and the RHNC integral equation theory, respectively. For hard sphere systems the exact closure is reduced to cQ/3(12) = -1^(12) - 1 for T < da/3, (2.26) so it replaces either equation 2.24 or equation 2.25 in this region. Other similar integral equation theories can be regarded simply as different approximations for bap(l2). In this study, Models I, II and III use the RHNC theory, with uncharged hard spheres as the reference system for all particles, and the exact functions 0R.ap(12) used are those of the Lee-Levesque extention [73] of the Verlet and Weis [72] fit to Monte Carlo calculation data. Model IV uses the potential of mean force from a pure solvent calculation using Model III, so at the Macmillan-Mayer level, the RHNC theory is used for the solvent. Model V uses the HNC theory, but the density of the ions in calculations using models IV and V is small, and the corresponding hard-sphere bridge functions would not contribute significantly. Chapter 2. Theory 23 2.1.4 R e d u c t i o n of the O Z E q u a t i o n It is possible to expand all functions of interest in rotational invariants. As was shown by Blum and Toruella [74] and by Blum [75,76], such an expansion applied to the func-tions of the OZ equation can reduce the equation to a set of algebraic relations between Hankel transforms of the r-dependent coefficients only. What follows is a summary of this reduction, with particular application to linear and C3v symmetry. The method of solution is one of finding r\ap (12) from a given ca/3(12), so the equations for the trans-forms used are given in terms of the function cap(l2) until the point where rjap(12) is determined, and in terms of 7/a/g(12) from that point on. Analogous transforms for the functions not given are implied. If we write the OZ equation explicitly in terms of its dependence on interparticle vectors, ^ ( 0 1 , ^ 2 , 1 1 1 2 ) = J^PfJ IVcc-y(Hulls. His) + Ca7(0l,H3)Ill3)] xar0(Q.3,Q2,r13 - r12)dr13dQ3, (2.27) noting that r 3 2 = r 1 3 — r 1 2 we can see the integral is clearly a convolution, suggesting that we Fourier transform the equation to obtain ^(•12) = g ^ X X / ftry(13) +ca7(13)]c7/3(32)dH3, (2.28) where caf3(12) = J cap(12)exp(ik -r12)drt2. (2.29) If the integrand of equation 2.29 is expanded in rotational invariants as was done for the interparticle potentials above, then the result of the transformation is also a rotational invariant expansion [74], 5^(12) = £ S X v ' W ^ & i ' & M (2:30) mnl Chapter 2. Theory 24 where $ is defined analogously to that in equation 2.3, and the fc-dependent coefficients are Hankel transforms of the r-dependent coefficients which are given by ,p(k) = Jo r23i{kr)c^{r)dr. (2.31) The ji(kr) are spherical bessel functions [77] of order I. The calculation of a large set of these projections is computationally time consuming unless fast Fourier transform (FFT) techniques are used, which is only possible to do directly for transforms of order I = 0 and / = 1. It is possible to use FFT integration for all other I values if a two step procedure is followed. As is discussed below, a two step procedure is needed for all orders of I if the logarithmic FFT method is used. The original r-space projections are first "hat" transformed, to yield (r) „mni co C nmnl Ft " mnl lUs]pofr_ where Pf(x) a n d Pf(x) a r e polynomials defined by, \t-i„2i ds for I even, ds for I odd, P2t+2ix) — ^y~i i=0 i=0 i\(t - + \)\ (-y-'x^jt + i + iy. -*)!(* + 1)! PS(z) = P°(x) = 0, for t > 0, for t > 0, (2.32a) (2.32b) (2.33a) (2.33b) (2.33c) and the usual non-integer factorial [77] tj! = T(q +1) is used. Now the Hankel transforms can be completed using FFT integration on the hat transforms [76], /•oo ClU f c) = 4 7 r / r2Jo(kr)c™Ur)dr for 1 even> (2-34a) J 0 Chapter 2. Theory 25 ~mnl CfJ.V /•oo %{k) = 4TTZ J T2j^)c^{r)dr for / odd, (2.34b) w. here Jo{kr) = —~J-, (2.35a) and . sin(fcr) — kr cos(kr) , , , *(*»•)= 1 J ( f c r ) 2 1 j - (2.35b) The rotational invariant expansions for (12) and cap(12) are then inserted into equation 2.28, and after angular integration and much simplification [74,75] the OZ equa-tion becomes, *l™Uk) = E P T E (~r ^ n t te7£r(*) - ^ ; 7 / 3 ( H (2-36) where the coefficients are given by mnni V / 21 + 1 2raj + 1 m n n^ (2.37) with {•:::} representing the usual 6-j symbol [62]. If the form of fmnl given by equation 2.5a is used, the coefficients reduce to \ zl±L = - ( - ) m 4 " + n i ( 2 / + i ) 1 / 2 / 1 1 m n ni (2.38) In this form of the OZ equation, a, 3 and 7 can represent any component in the system, all the correlation functions can be fully represented by restricting m,n,n\,l,l\ and 72 to non-negative integers, and the indices are restricted to those which give non-zero values for . These restrictions are that the triangle inequalities must be satisfied for the triplets (m,ni,/i), (ni,n,Z2), (m,n,l) and (Z2,/i,Z), and we must also have \p\ < m, \pi\ 5: ni and \v\ < n. Now using an identity [78] of the 6-j symbols we can rewrite the Chapter 2. Theory 26 Zmnni c o e r n c i e n t s as ^ =(2/ + i)B-)" ni n I2 \ X -x 0 ) m n\ li ^ ( X X 0 ) (2.39) the form of which suggests the use of the transform / <w(*) = E 1 m n l (2.40) V x - x 0 which is known as Blum's %-transform [75], with the sum over the I values that give a non-zero 3-j symbol. N™™£p(k) is defined in a similar way and these transforms are introduced into equation 2.36 to give = E ^ E ( - ) x +^ tasw + csn&w (2.41) This equation in turn suggests the use of matrix algebra, and following Blum [75,76] we let each C™™£p{h) be the (i,j)-th element of the matrix C ^ . For general symmetry, the positions of the elements are given by i = m(m + 1) + p + 1, (2.42a) j = n(n + 1) + v + 1. (2.42b) When the symmetry is such that this general assignment must be used, there can be many rows and columns of zeros in matrices representing terms for which m,n < %, and these rows and columns are removed for numerical ease. For most symmetries, however, many more of the elements will be zero, so even more compact assignments are possible. In particular for Models I and II of this study, a convenient way of assigning elements to matrix positions for symmetry, so that all matrix elements are non-zero, is m2 + m + 1 + .3. X 2 + 2 + 1, (2.43a) Chapter 2. Theory 27 n2 + n + 1 + v .3. X 2 + 2" + 1, (2.43b) where the square brackets indicate integer division. An even more compact assignment for Model III (linear symmetry) is, i = m - x + 1, (2.44a) j = n - X + 1- (2.44b) In order to write equation 2.41 in terms of matrices it is necessary to define two additional matrices. The first is simply /97 = p7I, where I is the identity matrix of the same order as the other matrices. The other is the matrix P with elements P{j all zero except i = m(m + 1) — p + 1 P{j = (-)" when j — m(m + 1) + p + 1 for general symmetry, (2.45a) = (-)" when < m2 + m + 1 > ' 3 .3. m2 -f m + 1 >" 3 .3. J = X 2 + 2 X 2 + 2 + 1 + 1 and Pn = 1 for Unear symmetry. Equation 2.41 is now written sL = E(-)x[&+£Jrf£x •y/3-for symmetry, (2.45b) (2.45c) (2.46) In general for multicomponent systems, these matrices are elements in larger matrices such as C X with (a.,8) elements £ ! ^ g - Corresponding definitions of the matrices R = diag(g,Zg,---), (2.47a) Chapter 2. Theory 28 g = diag(o±,pj, •••,o^), (2.47b) are also used, and substitutions into equation 2.46 give E x = (-)* [fi* + £ * ] ££fi x > (2-48) which is rearranged to give the form of the OZ equation which is actually solved, Rx =QxpRQx{(-)xI-PRO*}'1. (2.49) Thus, given the functions c™^(r) , the various transforms are performed to get the functions C^^(k) which are then used build the matrices The OZ equation is solved for £f which is dismantled into the functions N™£p(k) which are back %-transformed using ^ * ( * ) = M 2 i + l ) X ; f m n ZW"£(fc), (2.50) * V x -x o J where % ranges between ± min(m, ri): The fj^fap(k) a r e * n e n back Hankel transformed using 1 r°° t™Us) = ^ 2-J0 k2Uks)^(k)dk for /even, (2.51a) C;U*) = ^51 k2ji(ks)ij™lp(k)dk for /odd. (2.51b) Finally the real space r-dependent coefficients of the invariant expansion of 77(12) are obtained from the inverse hat transforms, CiflO") = V^U(r) - -1 jT s2P{ ( i ) for / even, (2.52a) V™Ur) = vXeW - h i ' S 3 P ° (;) iZltW ' ^ d , (2.52b) where P[(x) and P°(x) are the same polynomials used in the forward transforms (equation 2.33a-c). Chapter 2. Theory 29 2.1.5 Reduction of the Closure The closure equation is now written in the same rotational invariant language as the OZ equation. It is not possible to expand the logarithm in equation 2.23 in rotational invariants, so we must rewrite the equation to make this possible [48]. We begin by taking the partial derivative of the exact closure relationship given in equation 2.23, with respect to interparticle separation, r, giving the result, dca0(12) _ dha0(l2) 1 dgaP{12) 1 dua0{\2) + 56^(12) dr dr gap{!2) dr kT dr From equation 2.21 we can easily see that dga0(l2) _ dha0(l2) dr dr dr and if we define the quantity (2.53) (2.54) •«>c0(12) = ca/3(12) - ha0{12) + -7^/3(12), we can rearrange equation 2.53 to give dca0(\2) dr . = -ha0(12) dwa0(\2) dba0(12) dr dr 1 dua0{\2) dba0(12) kT dr dr We now can integrate the equation over r which gives dca0(12) dr dr — J°°ha0(l2) r kT Jr dwa0(12) dba0(12) dr dua0(12) dr dr dba0{12) '-dr + I kT Jr dr Jr dr dr. (2.55) (2.56) (2.57) The functions ca0(12), ua0{12) and ba0(\2) all go to zero as r —> 00, so the integration can be easily performed for those functions yielding, /oo ^ ( 1 2 ) dwa0{12) dba0(12) dr dr dr-^ua0{12)+ba0{12). (2.58) Chapter 2. Theory 30 The logarithm in equation 2.23 has thus been replaced by a product followed by an integration, both of which are processes that can be carried out with invariant expansions. Since the only parts of the invariant expansions that depend on r are the projections, the derivatives of pair functions can be simply written as expansions such as r ", Or ^ mnl (2.59) In general the product of two rotational invariants can be written $ ™ ^ ( 1 2 ) $ ^ ( 1 2 ) = EP M " ' $ ^ ( 1 2 ) > (2.60) mnl where the p coefficients are scalars given by [48] _mnl x < •Fmnl mi ni li m2 n2 h m n I (2m + l)(2n + 1)(2Z + m! m2 m n i n2 .71 ^ ^ y u1 v2 -v 0 0 0 (2.61) and | : : : | is the usual 9-j symbol [61]. The closure equation can now be expanded in rotational invariants and the r-dependent coefficients of equal invariants equated to give the following set of closure relations, mnl (r) = E TT\\Tl\l\lL \ U\ TTI2TI2/2M2 ^2 fC P - Jr dwm2n2hJr) Or at TT12T12/2 (r) dr dr umnl kT (2.62) for all allowed index values. The number of computations can be greatly reduced and the expression for the p coefficients can be simplified if we use the real space equivalent Chapter 2. Theory 31 of the ^-transform, 1 V x - x o (2.63) where the sum is over the allowed I values. Using similar expressions for dw™*™22.lJp(r) / dr and db™^}£p(r) I dr and substituting them in equation 2.62 gives, mnl /•oo m i n i x i "12712X2 <9r <9r (2.64) Now we use the choice of fmnl as given in equation 2.5a, and the expression for the scalar coefficients becomes, P™* = [(2mx + l)(2m2 + l)(2m + l ) ^ + l)(2n 2 + l)(2n + 1)]^ (-)"+" / m „. _ \ / _ _ _ \ / x mi m 2 m \ Xi X2 - x mi m 2 m rt\ n2 n \ " X i - X 2 X rii n2 vi v2 n \ (2.65) This is useful because the number of non-zero coefficients is reduced enough to more than offset the computation time incurred by forward and back ^-transforming, as well as increasing the accuracy of the result. The back transforms are given by an expression analogous to equation 2.50. Chapter 2. Theory 32 2.1.6 The Static Dielectric Constant and Average Quantities. For pure ammonia we make use of three different expressions in order to obtain the static, or equilibrium, dielectric constant, e. The first is the Kirkwood relationship [79,80] ra = (^ ifi±il (2.66a) where [81,82] 9 = 1 + ^ j f : r X Z s ( r ) d r , (2.66b) is known as the Kirkwood g-factor and y = -MT, (2- 6 6 c) where p represents the dipole moment of the pure solvent molecules. The second route to e is via the limit [83] KlUr) -> ( £ 1 ) 2 3 , asr.-^oo, (2.67) with y given by equation 2.66c. The dielectric constant can also be given by the expression [84] l - ^ ( l i m f c ^ o £ S ? „ ( f c ) + 21im^o^.(fc)) flo ^ e = (2.68a) 1 - * ( H m ^ 0 c™.(k) - Hm f c_ 0 <&?„(*)) where the first limit can be found by [OO Jim <&?„(*) = 4TT / r*clZs(r)dr, (2.68b) «—»0 Jo and if we note that there are no short range contributions to the k —» 0 limit we can write = < 2 - 6 8 c > It should be noted that A; in the product kT on the right side of equation 2.68c is the Boltzmann constant, and should not be confused with the Fourier space variable k on the Chapter 2. Theory 33 left side. These three routes should yield the same result if the OZ equation is satisfied, and so can be used as a test of convergence of the numerical solution. The value reported in the next chapter is taken from equation 2.66a using the trapezoidal rule for integration. The integration was also performed using Simpson's rule and the discrepancy between these methods was found to be negligible. For electrolyte solutions the first two routes (equations 2.66a and 2.67) are no longer correct because of the effect of Debye screening on the projections involved. Equation 2.68a is still valid however as is the equation [84] e. = Zyg + l (2.69) where e„ denotes the dielectric constant of the solvent, and y and g are the same as above. The Stillinger-Lovett [85,86] second moment condition also furnishes a route to e„ for ionic solution, and it can be expressed by the equations Air —* - (2) e* = -fiji .2s PiWiWj (2.70a) and W = ~r-f hZ+rydr, (2.70b) where i and j represent ionic species. For electrolyte solutions is a well defined quantity [34,63,81], however it can not as yet be unambiguously determined from experiments due to dynamical contributions that are not easily calculated theoretically [87-89]. The total instantaneous configuration al energy for a mixture can be expressed by utot = ^zZzZzZvaeiiJ), (2-71) ot/3 i i where Na is the number of particles of species a and ua@(ij) is as defined in section 2.1.2. If we now use the general expression [14,90,91] for the average of any mechanical Chapter 2. Theory 34 quantity map(l2) corresponding to the species pair a/3, (2.72) with V representing the sample volume, we can write the total average configurational energy as ^ 1 ^ 7 ^ J g^i^u^^d^d^dr, (2.73) U> TOT N in which pi = N/V denotes the total number density and xa the mole fraction of species a. If we expand gap(12) and uap(l2) in rotational invariants, and simplify using orthog-onality of the invariants (equation 2.11) we get UTOT _ 27T „ mnl (2m + (2.74) In some cases average energies of specific component pairs Uap are given, and in these cases it is convenient to use N„ 2w(2 - M P / s X mnl ^-tmnl^2 too L(2m + l)(.2n + l)(2Z + l)./o r* <"UrK?Ur)dr , (2-75) where Kroenecker's 8 has been used. The only non-zero contributions to the sum are those for which m+n = I, because the only non-zero projections of the pair potential have this restriction (see section 2.1.2). The non-zero energy contributions can be attributed to the interaction between specific multipoles in the expression for each u™£?ap(r), so individual interaction terms are expressed as a|3;nra ~N7~ 27r(2 — Sap^m +n,lPP {fm Y f°° (2m + l)(2n + 1)(2Z + 1) Jo r*9™UrK?Ur)dr (2.76) where the subscript a.8; nm indicates the interaction between the 2"-pole of species a and the 2m-pole of species 8. For ionic systems, the integrals representing ion-ion energies in Chapter 2. Theory 35 equations 2.74-2.76 diverge individually, but in charge neutral systems the divergences cancel exactly and the remaining contribution is the quantity of interest. Whenever these quantities are discussed it is understood that it is the remaining non-cancelling contribution to the total energy which is given. The compressibility factor Z = pV/'NkT is found in a similar way [14] from * = 1 - m ^ ^ W Y Ir^(l2)^^dnidn2dr, (2.77) which can be derived from the virial equation. The same situation with divergent integrals occurs here as for the energies, but once again the divergences cancel. We also note the derivative of the potential function at the hard-sphere discontinuity is a (negative) dirac 8 function. The coordination numbers around ions, solvent molecules and macrospheres will be needed for the discussion below. We define the coordination numbers as Cn = 4*p. [Rr2g™B(r)dr, (2.78) where j represents s, + or —, and R denotes the separation between the species concerned where <7oo°js(r) n a s ^ s first minimum. The coordination density is defined simply as Chapter 2. Theory 36 2.1.7 The Self-Consistent Mean Field Treatment of Polarizability Effects. The dependence of electrolyte properties on dielectric constant is extremely important for this study, and it has been shown [92] that neglecting polarization effects leads to a significant underestimation of the dielectric constant of pure ammonia. Similar conclu-sions can be drawn from the recent literature on models of water and aqueous electrolyte solutions [34,38,59]. A recently published simulation [39] of a model of ammonia1 has shown that remarkable agreement between calculated and experimental [40] distribution functions can be achieved by including molecular polarizabilities. Clearly polarization effects on equilibrium properties can not be ignored at liquid densities. Calculating polarization effects exactly is a difficult problem, where many-body in-teractions have a significant influence. It is not possible to treat these effects exactly in the present integral equation formulation, so we must reduce the many-body interactions to pairwise additive effective potentials. The SCMF theory is an approximate theory [35], but it has been shown to be a relatively accurate way of reducing the many-body problem into an effective pair potential, for models of both water [34,38,59] and ammonia [39,47,64,92]. The SCMF theory also has a physical basis behind its mathematical formu-lation with approximations corresponding to distinct physical concepts. The polarizing effect of placing a molecule in an electric field due to the close proximity of other molecules is a time dependent quantity, but the effect on the system as a whole is approximated by replacing all the changing dipoles with an effective permanent dipole which is consistent with the average electric field felt by a molecule when all other molecules in the system are treated in a similar way. This is done by. ignoring fluctuations in the local electric fields experienced by the molecules. All multipoles in all species, other than molecular dipoles, are assumed to be constant, but their contributions to the average electric field 1The model is similar to Model III with a Lennard-Jones potential in place of the hard-sphere repulsion. Chapter 2. Theory 37 are included. We note that the SCMF theory is formally equivalent to Wertheim's 1-R theory [93] when only dipoles are used. Since the solvent is the only species considered to be polarizable, the subscript i throughout will refer to the solvent only. A polar-polarizable particle will have an in-stantaneous dipole moment m ; , which can be written as, ULi^p.+p., (2.80) where p. is the contribution from the permanent (constant) dipole moment of particle i, and p. is the instantaneous induced contribution. The induced part is caused by (J?;),-, the local electric field which is felt by particle i , so we can write, £ = (2.81) where Q is the polarizability tensor. If we now define the average total dipole moment m', we have, m'-=(m) = (/£) + ( a - ^ ) = M . + a-C&>- (2-82) For C 3 „ as well as linear symmetry, the average local electric field will be non-zero only in the direction of p, which in turn implies that m' will be in this same direction so we can write, (Et) = C(m')m', (2.83) with C(ra') being just a scalar. If we substitute equation 2.83 into equation 2.82 and then iterate the result, rn = p +e • C(m')m', (2.84) we get the infinite expansion, m' = p -f Q • C(m')p + Qi • Qi • C2(m')p + Cz(rn)p + ••• , (2.85) Chapter 2. Theory 38 which suggests the introduction of the renormalized polarizability a' into equation 2.85 giving, m' = p + C{m')^ • p, (2.86) where, a' = £ + C(m')a' • ct. (2.87) The renormalized polarizability a' of an ammonia molecule in a polarizable fluid describes the fluctuations in the total dipole about its mean in the same way as the original gas phase polarizability a did for an isolated ammonia molecule, so we can write = (m')2 + 3a'kT (2.88) m 2 N where a' = ^Tra/ (2.89) and (m2) is the mean square total dipole. We now wish to evaluate C(rn') in terms of a_' and m' so all three quantities can be found. Comparison of equation 2.86 with equation 2.84 shows that a particle with polarizability a which experiences a local field due to the surrounding particles behaves as if it had the renormalized polarizability ct'. Therefore, we define an effective system analogously to that described in section 2.1.2, in terms of interparticle pair potentials containing products of pairs of multipole moments, with the total system potential given by ^ = X ^ ( 1 2 ) + ^-EPn, (2-90) where the sums are over all particles of all species i and j , and polarizable species n. In the second sum, representing the polarization energy, we ignore fluctuations in the induced dipole, substituting p1, by (p2). The sum becomes N (p2) /2a, a constant, and so Chapter 2. Theory 39 it can be dropped. All multipoles except dipoles are assumed to be constant so for terms not involving dipoles the interparticle potentials remain unchanged. The dipole-dipole potential term for a single pair of particles i and j contains the product m,mj, which for a system where fluctuations in the instantaneous dipole moments are ignored can be replaced with an effective permanent dipole moment, m e , defined for all solvent particle pairs as rrnmj = (m 2 ) = m 2 . (2.91) The potential terms for interactions between dipoles and other multipoles involve prod-ucts such as 77i; 0 or m^O and again when fluctuations are ignored we can replace m; by m ' . In order to simplify the expression for the potential of the effective system, we also make a further approximation by assuming m ' — m e , so we can now describe the effective system simply as a system with an effective permanent dipole. We note that this approximation is not used in any of the expressions below, but only in the calculation for the effective system once m e and m' are known, and the maximum difference between these quantities for any calculation used here is 3.8%. We find explicit expressions for C(m') by writing equations for the total average configurational energy for both effective and polarizable systems and equating separate terms. The instantaneous configurational energy for the polarizable system is given by, U N — U H S + U U + U I D + UJQ + uIO + uIH + uDD + UDQ + UD0 + UDH + UQQ + UQO + UQJJ •+ u0o + u0H+uHH+uP, (2.92) where UHS is the hard-sphere energy, up is the energy of polarization and all other terms are the energy of interaction between multipoles of a pair of particles with I,D,Q,0 and H representing ionic charge, dipole, quadrupole, octupole and hexadecapole, respectively. In species where multipoles are not defined, the corresponding terms in this expansion Chapter 2. Theory 40 are zero. Since only up and the dipole terms are affected by polarization we write them in terms of their dependence on local electric field, and replace the sum of the remaining terms by to obtain, uN = - E & j • (&i)i - \zZV±i • (RID)j ~ E ^ ' (—«?)i i i j - !>;• (Mw)i -Y,Bkr i&H)i + \ EE,- • (2-93) 3 3 3 where the sums are over the j solvent particles and (£,); = (£,/),- + {R^ + ( £ / < ? )j + (Elo)i + ( l u r ) ^ (2.94) with (En )j being the contribution to the total local electric field from all the ions in the system, and all other terms on the right hand side of equation 2.94 representing analogous quantities. The average total energy is then UTOT = ^ t - i V J ( ^ . ^ / ) - i ^ ( ^ . ^ ) - i V s ( / £ . ^ Q ) - i V 5 ( / x . E i ( - N. (/* • ElH) - ±NS (p • Eu) - ±N. (p-Eiq) -\Ne(p-Elo)-±Ne(p.ElH), (2:95) where = £ 7 * . Now, ignoring fluctuations in the local electric field we write {P-E1) = (P)-(E1), (2.96) which enables us to write equation 2.95 in the form UTOT = £ / f - \N.p - ±N.{m' + p) [(Eu) + (lhQ) + ( ^ 0 ) + CEur)] • (2:97) The equivalent expression to equation 2.97 for the effective system can be written, UTOT = U*- ]-N.me (E,D)e - Nsme [ ( ^ 7 ) e + + ( E i G ) e + (EjH)e} . (2.98) Chapter 2. Theory 41 We now write C(m') = dim') + CD(m') + CQ(m') + C0(m') + CH(m% (2.99) and from equation 2.94 we have (EJJ) = (7,(m')m', (2.100a) (Ew)=CD(m')m', (2.100b) EiQ) = CQ(m')m', (2.100c) (Ejo) = Co(m')m', (2.100d) i&H) =CH(m')m'. (2.100e) A n d finally, because both polarizable and effective systems are assumed to have the same structure, it follows that ( £ a - ) e = ( £ z x > , . (2.101a) for X = I,Q,0 and H, and (ElDy = ^ (Kto) • (2.101b) m We just have to identify the individual terms in equation 2.98 with terms in equation 2.76 giving UEDD=-\NBmE(EW)\ (2.102a) and UEDX =-NsmE(E_LX)\ (2.102b) again for X = I,Q,0 and H, and apply equations 2.100a-e and 2.101a,b to get the final result, Chapter 2. Theory 42 — 1TJe Cr>(m') = — ^ , (2.103b) Nsml Co(m') = - ^ - r (2.103d) I\smem' CH{rn') = -^™-r (2.103e) Nsmem' Now it is only necessary to find m e , m', a/ and C(m') self consistently. Chapter 2. Theory 43 2.1.8 Angular Probability Distributions. The interpretation of angular distributions near a macrosphere or hard wall will be considered in subsequent sections, so the definition and method of calculating these quantities is given here. The two vectors in the molecular coordinate frame which will be used to interpret the orientation of an' ammonia molecule with respect to a surface normal vector are the molecular dipole vector and the molecular NH-bond. First we write the expression ( r ) \ = /^( i2)^% p ( i2)dn 1 rfn 2 (/mnT (2.104) (2m + l)(2n + l)(2* + l) [g°0^(r) which is the average value for a rotational invariant at a separation r for a given distri-bution gap(l2). We require the quantity pT(0), the probability density function for finding a molecule with a fixed vector v at an angle 6 to the surface normal vector, at some fixed separation r. We can use the rotational invariant language given above to introduce an angular distribution function, ^(12), for the general vector v fixed in the molecular frame. The angle 6 with which we are concerned is the polar angle for this vector, as well as one of the Euler angles for the distribution, so we can simply write ,n^n SfdxSfd* [g™, (r) + £ h^(r)*™1 (X,0,</>)] dS PrMde = — s^v(r) ' ( 2 - 1 0 5 ) where (%, 6,<j)) are the Euler angles of the distribution, the sum is over all (m, n,Z, p,u) except (0, 0, 0, 0, 0) but restricted by the relations given in section 2.1.2, and the subscript v distinguishes this distribution from that of the ammonia molecule itself. Since the orientation of the vector does not depend on % or (/), this expression reduces to PrA0) = l+Y^(KL^))Kol(S), (2-106) Chapter 2. Theory 44 where we have used equation 2.104. If we now use [94] (Pl(cos6)r) = ($°00l(6)T), (2.107) equation 2.106 becomes (2.108) where (P;(cos #)r) is the average of the Ith Legendre polynomial of cos 8 at a distance r. This probability density function satisfies the nomalization condition and in a random system with no correlations pT(6) = 1/2. It will be interpreted as the average orientation per particle, which implicitly assumes there is little correlation between the average orientation and average position of a particle. We can calculate the angular distribution for the molecular dipole in a straightforward manner by equating the molecular Euler angles (a , /3,7) , to the required vector Euler angles, (%, 6,•</>). This also equates #(12) for ammonia to <7„(12) for the vector, and equation 2.106 can be evaluated directly. For the case of the NH-bond distribution, or in general any intramolecular direction for which 9 is not equal to one of the Euler angles of the molecular distribution, we must calculate each of the (P;(cos 6)r-v) in the following manner. The rotation matrix which carries any vector in the molecular reference frame xyz into the reference frame XYZ of the vector v is (2.109) cos {3 cos 7 — cos (3 sin 7 sin /3 R = sin 7 cos 7 0 (2.110) — sin (3 cos 7 sin (3 sin 7 cos (3 If we now denote the cartesian coordinates of the vector in the molecular reference frame of an ammonia molecule as (xa, ya, za), and if we note that we wish to rotate the molecular Chapter 2. Theory 45 NH-bond into the Z-axis, the relationship we require is given by the bottom row of the rotation matrix, cos 6 = — xa sin3 cos 7 -f ya s m 0 sin 7 + za cos 8. (2.111) Now that cos 6(6, f) is known we can evaluate 1 /-2TT r7r hmnl (r) (P,{COB9U) = - d1\ Bin/9P,(coBgro(A7))4<9 + E , " 0 0 0 A 47r Jo Jo 47T5r 0 0 ; m (r) x Jo d~f jo rin/9P,(cosfim(/3>7))*™nI(i8l7m (2.112) where once again the usual subscript restrictions apply and the sum is over all values except when they are all zero. One final relation, the average cosine of the angle between the dipoles of particles 1 and 2, (cos^) = ( £ 1 . / a 2 ) , (2.113) where p is a unit vector in the dipole direction, will be required. This quantity can be identified [59,81] with {$oo°(r)) a n < ^ s o *s g i y e n directly by equation 2.104. o Chapter 2. Theory 46 2.2 Theory for Ionic Solutions. This section will be concerned with the theory of electrolyte solutions in the form given by Kusalik and Patey [95] as extended from the original formalism of Kirkwood and Buff [58] for mixtures of uncharged particles. The need for this extension of Kirkwood-Buff theory stems from the fact that in electrolyte solutions the neutrality requirement means the concentrations of individual ions can not be varied independently. When the Kirkwood-Buff theory is directly applied, in conjunction with charge neutrality condi-tions, the expressions for all thermodynamic properties are indeterminate. The nature of the extension is one of finding well-defined thermodynamic quantities for general molecu-lar solvent systems as was done by Friedman and Ramanathan [96] for continuum solvent models. The results given here will be for a single ionic solute of the general form Av+Bv_ , but could be extended to systems with more components. Throughout this section the parameters q+ and g_, related by q+ = —q_v_/v+ will be used to denote the charges of positive and negative ions, although in all calculations only monovalent ions are consid-ered. The number density of the solute will be denoted p2 which is related by p+ = v+p2 and p_ = V-P2 to the individual ion number densities. Of particular interest in this study are the mean ionic activity coefficients. These are defined, using the notation of Harned and Owen [97] in terms of the chemical potential of species a as Poc = P°a + kThxaa = /^°;m + kT'ln(7Qma) = vl;c + kT\n(yaca) = p°a.tP + kTln(yaPa), (2.114) where aa is the activity of species a, and ma and ca are concentrations in units of Chapter 2. Theory 47 molality and molarity, respectively. The quantities 7 a and ya are activity coefficients which depend on the concentration scale used. In order that both the activities and the chemical potentials of a given species be independent of concentration scale it is necessary to define standard states (superscript 0) which conform to this convention, so each choice of concentration scale has its own unique standard state. The standard states are defined in such a way as to allow all the activity coefficients, regardless of scale, to approach a value of 1 in the limit of low concentration. We point out that the same activity coefficient, ya, is used for both molarity and number density concentration scales. Single ion activity coefficients can be formally defined using u-2 = v+P+ + v-P-, (2.115) but in order to compare with experimental data, the mean ionic activity coefficient, y±, is defined (for the molarity scale) as y, = (iCy- )•': (2.H6) where v = v+ + •!/_, and from which we have p2= plP + kT\n(u^u1-) + iykTHy±p2). (2.117) 2.2.1 Kirkwood-Buff Theory for Electrolyte Solutions. The formalism of Kirkwood and Buff [58] gives exact expressions for the thermody-namic properties of an n-component system in terms of a matrix B with elements defined by Bap = PaSap + PaPpGafi,. (2.118) where 8a@ is the Kroenecker delta and the quantities Gap, known as Kirkwood G's, are given by roo GaP = 4TT / r2ha0(r)dr, (2.119) Jo Chapter 2. Theory 48 where hap(r) = h^al3(r) is the radial correlation function £r™°a/3(r) — 1. Charge neutrality in electrolyte solutions can be expressed in terms these Kirkwood G's by the equations YPi9iGii = -qj1. (2,120a) i and X > 9 ; G i s = 0. (2.120b) i If we now label the solvent as component 1 and the other species as components 2,... ,771, and introduce the quantity m S = X PaP /s lS let f , (2-121) a,/3 = l where |B | a jg is the cofactor of matrix element Bapy we can write the quantities of interest as fcTXT = | l l l , (2.122a) m ^ 7 = 1^-) = 4 X / M f i U (2.122b) U1*y/ T,P,Na ° /3=1 V / cty*a \ 1 • / 5/i a\ |B la/3 and kT\dNp)TytKn kT\dP,JTpi HI V:fdua\ \R\cf, Va% (2.122c) (2.12'2d) kT \dNpJTpN^ |BJ kTXT These quantities apply to ionic solutions only in a formal sense because single ion quan-tites are not determinable thermodynamically, but they are nevertheless useful because they lead to physically meaningful quantities for electrically neutral solutions. If these quantities are evaluated in terms of the Kirkwood G's and the charge neutrality condi-tions are applied directly, the results are indeterminate. It is necessary to find a way of expressing matrix quantities such that when the charge neutral limit is applied the result Chapter 2. Theory 49 is well defined. Following Kusalik and Patey [95] this can be done by rewriting the ma-trix elements in terms of their related corresponding fc-space quantities. If we recognize that Gap = hmfc_0 hap(k) = h^p, which comes directly from the Hankel transform as in equation 2.31, and that at finite ion concentration hap{r) decays exponentially at large r (see section 2.2.3), for small k we can write (2.123a) where hKa0 is given by equation 2.70b, and we can thus introduce the matrix elements Bap{k) = pa6a/3 + pappha0(k), (2.123b) which in the limit as k —> 0 become the matrix elements defined by equation 2.118. It is now only necessary to find the correct limits of these equations and apply them in equations 2.122a-d. For the present system of three constituents, as k —> 0 we have [95,59], \R(k)\is -» 0, fori = + , - , 5 , (2.124a) which are then used to show that ~ 2 „ 2 P + P-Ps Dk2 + (2.124b) and, if S(k) is defined analogously to equation 2.121 with terms of the form of equation 2.123b, we find that S(k)- « 2 n 2 „ 2 P + P-Ps Ps Dk2 + (2.124c) where D = h(+l + h{2)_ - 2h{2l. (2.124d) These relations are used in equations 2.122a-d and when the ratios are simplified before the —• 0 limit is taken, the quantities become determinate. Thus we have, 1 + P> (Gts - G+s) P2 [1+/>.(<?„+ G + - - 2 G + < ) ] : (2.125a) Chapter 2. Theory 50 and V. = G+- — G+a l + P . (G„ + G + _ - 2 G + . ) ' G+.+P. (G+_Gee-Gle) = l + Pa{Gaa + G+^-2G+ay kT [dp. V ( dp2 1 + psGa T,p. pl [G+-+P.{G..G+--Gl.)Y kT \3N2 J T p N a pl [1 + P a (GBa + G+- - 2G+S)} (2.125b) (2.125c) (2.125d) (2.125e) In order to obtain mean molar activity coefficients for comparison with experiment we take the partial derivative of equation 2.117 1 dlny±\ = _ L _ fdpA _ _1_ dp2 )TJ ukT\dp2)T. p2 (2.126) where j can be either ps or the pressure p. Using equation 2.125d we have immediately 1 + psGaa din y±\ _ 1 dP2 )TiP> ~ p2 vp2 \G+-+pa{GsaG+--G\a Evaluation of the constant pressure derivative requires the ratio V ! dfx2 \ kT\9p2)TiP~ V and the identity kT\dp2 L v(^) dp2 \ 1 dN,l " VpsVe' Now using equations 2.128a,b, 2.126 and 2.125e, we have 'd\ny±\ = _1 dp2 )Tj> P2 1 upa(G+.-G+.) (2.127) (2.128a) (2.128b) (2.129) Chapter 2. Theory 51 A final well known relationship [70,71] follows immediately from equations 2.118, 2.121 and 2.122a if a pure solvent system is considered. Since the only element of the relevant matrix is Bas = p„ + plGts, S reduces to p2B and we write kTXT = - + GBS. (2.130) P» 2.2.2 Stability Criteria for Ionic Solutions. The criteria whose behavior mark the onset of absolute instability for ionic solutions have been derived2 by Ursenbach and Patey [100] and the results as applied to the models in this work are stated here in terms of the Kirkwood G's defined in the previous section. The mechanical stability criteria are ^ = T T ^ T T - > 0, (2.131a) and > 0, (2.131b) and the compositional stability criteria are 'T,p and (dfia\ \ d x a ) T y DG+VC ^Mf = J - > 0, (2.132a) G G > 0, (2.132b) DGSG \ d x a J T y where xa is the mole fraction of species a, Ap = p2 — Us, the difference between the chemical potential of the solute and the solvent, and id denotes the quantity for a system behaving ideally. The three parameters DG, SG and VG depend in this case on the model in question. For Models I, II and III they are given by DG = vp2 + pB (G+_GSS - G2+B)} , (2.133a) 2 For a discussion of instability see reference [98] or [99]. Chapter 2. Theory 52 SG = 1 / 3 2 [1 + p, (G+_ - 2G+. + G„)] , • (2.133b) and VG = i/a^jc. [1 + p . (G„ - G+.) - (G+_ - G + s ) ] 2 . (2.133c) For Models IV and V they reduce to DG = vp2G+- (2.134a) SG = u (2.134b) and Vb = 0. (2.134c) Chapter 2. Theory 53 2.2.3 Limiting Behavior In the Debye-Hiickel theory for ions in continuum solvents the ion-ion potential of mean force, it;+_(r), is equated to the electrostatic potential and is obtained by solving the Poisson-Boltzmann equation [1,11]. In order to obtain an analytical solution, and in order that the principle of superposition of charges not be violated, the Boltzmann distribution is linearized. This is an approximation that becomes exact in the limit of large r and small ionic concentration. The solution to the differential equation is Q w + _ ( r ) = — exp(—/tr), (2.135a) where the Debye screening parameter K for a single salt solution is given by 4.7T \q+q_\vp2 (2.135b) .ekT1 and e is the dielectric constant of the pure solvent. The constant C is evaluated using boundary values. One method of finding C is to assume that at long range the electro-static potential due to the central ion i will approach that of a point charge, qi/ekTr. This will also become exact in the zero concentration limit, and the result which follows is one expression of the Debye-Hiickel limiting law, w+^DH(r) = Z | ± i i e - » - . (2.136a) It has been indicated that this potential of mean force is correct in the limit as r —> oo and as p2 —• 0 for both continuum and molecular solvents [81]. Another boundary value which has been used to evaluate C is to consider the extent of the ions. If we use the condition that the total charge outside a central ion must be equal and opposite to the charge of the central ion it follows that -a, a e~^T~d±^ w = X ( i + ^ r (2136b) Chapter 2. Theory 54 where d± is the distance of closest approach of a cation and an anion. It can be seen that as the concentration (and therefore K) becomes small, equation 2.136b approaches the same limit as equation 2.136a. We can use the form of w +_(r) in equation 2.136a to define the short-ranged effective ion-ion interaction, wsR-apif) (see equation 2.16) which is used in the description of Model IV. In the infinite dilution limit K —> 0, so we have - ( r ) - 1 ^ asr-oo , (2.137) ' ekTr which is the contribution from direct ion-ion interactions only. We therefore define w8Rx(r) = w™ij(r)-^, (2.138) where w^?-(r) = — In [^ So°ij(7')] *s taken from a full molecular model at infinite dilution (in this case Model III), and i and j represent (possibly the same) ionic species. Sub-tracting the direct ion-ion interaction leaves the short-ranged ion-solvent contributions as an effective ion-ion term. The function wsR;ij{f) is now used in equation 2.16 to define Model IV. We can use the potential of mean force given in equations 2.136a to obtain the ion-ion total correlation function = exp[—w+_(r)/kT] — 1 exp lg+g-l e-itr 1, (2.139) ekTr for the DH expression, and use the Kirkwood-Buff formalism to get thermodynamic quantities at low concentration. If we expand the outer exponential in equation 2.139, keep terms up to second order and apply the resulting expression in equation 2.119 we get 1 A G+-.j)H = + - = + •••, (2.140a) VP* JVp2 C h a p t e r 2. Theory 55 and the corresponding result G+-4 = — + VP2 • y/vp~2~(l + Kd±) 2 for the d-dependent relation (equation 2.136b), where + (2.140b) (2.140c) 2 \ekT J In order to obtain expressions for G+a and Gea in the infinite dilution limit we must introduce analogous quantities to the Kirkwood G's. We define K—»0 = 4TT / Jo (2.141) and appeal to the OZ equation (2.36), which we rewrite in terms of /i(12)aig and c(12)Q|g for the relevant projections, giving ^00°a/3 CO ^00°a/3 ( k ) ~ • E P 7 E ^ r E (-rfe( f c )^-«( f c ) ifj,~—m (2.142) where the Zomm coefficients are given by equation 2.38. Since all the h^™p(r) (m > 0) decay exponentially for electrolyte solutions at finite concentration, we have l i m / C ^ ( f c ) = 0, for all m ^ 0 which simplifies equation 2.142 to (2.143) Gaa Gag ^ ^ Pn/ Gay C*yg , 7 which, together with the charge neutrahty conditions (equations 2.120a,b) gives 'v+C+. + v_C. (2.144) G+a — G_„ 1 - paC„ P*G+~. (2.145) Chapter 2. Theory 56 Since G+e does not diverge in the infinite dilution limit, one can always find a low enough non-zero concentration such that equations 2.127 and 2.129 are dominated by which leads to the limiting law for the derivative of the activity coefficient. In the limit as p2 —> 0, d In y+ \ i d In y+ \ y ± 4 ' y ± 1 (2.146a) and (d\ny±\ _^ -Ay/i for the DH hmiting law, and fd\ny±\ _^ -Ayfi \ 9P* 1T,p ~* y / P * ( 1 + K d ± y QP2 JT.O. \ 9P'< _ (2.146b) (2.146c) for the corresponding d-dependent result. To obtain the usual limiting law expressions in terms of \o.y± equations 2.146b and c are integrated to give -AJvp~2 l n y ± = Y (2.147a) for the DH hmiting law, and — A^/vpn ^ = WT±) ( " 4 7 b ) when the d-dependence is included. 2.2.4 Individual Ion Partial Molar Volumes. By combining equations 2.125a, 2.130 and 2.145 the infinite dilution limit of the partial molar volume of the solute can be shown to be V2° = u+kTX°T(l - p,C°+3) + u_kTX°T(l - p.C°_,), (2.148) where the superscript 0 indicates the infinite dilution limit, and therefore XT *S *NE isothermal compressibility of the pure solvent. This solute partial molar volume is a well-defined quantity, but as has been pointed out by Kusalik and Patey [34], the individual Chapter 2. Theory 57 ion partial molar volumes are somewhat ambiguously defined by Vi = v+V++v-V-, (2.149) which leaves room to arbitrarily split the experimentally available total into almost any proportions. A number of physical assumptions are used in order to get the individual ion quantities and some of these methods are given in reference [101,102]. As has been shown [34], within the framework of the Kirkwood-Buff theory there are at least two reasonable theoretical definitions of these quantities. The first is done by comparing equations 2.148 and 2.149 which suggests the definition i ? - = ferXr(l-P.<7?).' (2-150) It can be shown that this is the quantity that is given by the Kirkwood-Buff theory in the limit of infinite dilution, if the limit is taken while charge neutrality is enforced. It is most likely that this quantity is closely related to experimentally determined values because charge neutrality is rigorously enforced by nature. Another definition that is possible if charge neutrality is not enforced is to define V® as the volume change when a single ion is added to an infinite bath of pure solvent. While this is not possible experimentally, it is nevertheless this single ion picture which is used to interpret experimental results, and is therefore an appropriate definition. Using the Kirkwood-Buff expressions for a two component system of solvent and ions all of the same charge qi, and taking the pi —> 0 limit the immediate result is V? = kTx°T - G°u- (2-151) The relationship between F^0* and is shown [34] by writing equation 2.151 in the form V? = kTX°T(l - P.C? - p.T?), (2.152a) Chapter 2. Theory 58 where and where e is the dielectric constant of the pure solvent. The quantity T°s depends only on ion charge and pure solvent properties. For ions of equal but opposite charge, T°e = — T° a so for charge neutral systems these quantities will always cancel exactly. The well-defined quantity V2° is always recovered no matter which of these two definitions is used. This term is a model dependent term which arises from the coupling between ion-solvent and solvent-solvent orientational correlations which occur only when an individual ion is unscreened. As will be seen below, it is the definition which includes the coupling of correlations which is in agreement with physical notions of electrostriction, while the alternate definition, V®~ is not. Finally we note that the quantity T? must be non-zero for these two definitions to be distinct. For models which require that the solvent-solvent projection c^}at(r) be identically zero (such as purely dipolar solvents), we have t^°* = t^°. This projection is non-zero for N H 3 , however, and the distinction between definitions is important. Chapter 2. Theory 59 2.3 Ideas of Ion Pairing. An ion pair is conceptually easy to envision, but an unambiguous mathematical de-scription requires defining a new species in a somewhat arbitrary way. The question "What is an ion pair?" has no definitive answer, but it is considered to be a distinct species held together by electrostatic forces only. The mutual potential energy due to the electrostatic interaction of oppositely charged ions is proportional to q+q_/r and as the ions approach each other it becomes large with respect to the thermal energy of the molecules of the solvent. Thus it is supposed [11] that the closer a pair of oppositely charged ions approach each other, the more stable the species will be, and hence longer lived. If there are other short ranged interactions, however, the molecular picture could be quite different, and solvent separated pairs or other clusters of species could be more stable, and this has in fact been observed in such models [53]. The distance an anion and a cation must be within before they are considered an "ion pair" is clearly an arbitrary decision. Statistical mechanical theories of ionic distribution [10-12] generally define ion pairs in terms of a distribution function (r) and a distance R. Ions of opposite charge are considered associated if they are within this distance of each other, and the average number of ions within this distance can be calculated from the distribution function. Pair functions do not distinguish between ion pairs and other associated groups of ions such as ion triples so an assumption that ion pairs are the only significant contribution to g+-(r) for r < R is usually made in order to calculate the proportion of associated ions. 2.3.1 The Bjerrum Theory. The Bjerrum theory [10,11] was the first statistical mechanical theory in which inter-ionic forces in electrolyte solutions were used to calculate the proportion of "free" ions a Chapter 2. Theory 60 or the proportion of ion pairs 1 — a. Bjerrum began by assuming the number of ions of type i contained in a shell of width dr at a distance r from a central ion of type j (with charge opposite to type i) is 4 W 2 exp ) dr, (2.153) with u(ij), the mutual potential energy of the ions, given by (ij) = m - (2-154) u er By comparison with the Debye-Hiickel potential of mean force given in equation 2.136a we can see the potential energy is only that of a pair of charges in a medium of dielectric constant e, and the contribution to the potential energy from the ionic atmosphere has been neglected. As r increases, equation 2.153 has two competing terms. The exponential term decreases and the volume of the shell, r2dr, increases. This produces a minimum which can easily be found by differentiating and equating the result to zero. The value of r which corresponds to the minimum local charge density of ions of type i (opposite in charge to the central ion of type j) is ~ - Qiqj' = R. . (2.155) n u n 2ekT Bjerrum then defined ion pairs as those ions within a distance R of the central ion. This corresponds to defining ion pairs as those for which the electrostatic potential energy is greater than 2kT. Now by integrating the function in equation 2.153 from contact to the separation R the degree of association (1 — a) is found, 1 - a = ATPi f* r 2 exp (^ jg*) dr. (2.156) If one defines the association constant in the usual way, Chapter 2. Theory 61 a relationship between y± and KA consistent with Bjerrum's definition of ion pairs is possible, and hence an association constant for the low concentration Hmit. Expressions for activity coefficient immediately follow [3] from equations 2.147a and b by assuming that ions outside the distance R obey the Debye-Hiickel law for unassociated ions. We have, lny± = v (2.158a) using the DH limiting law, and + < 2 1 5 8 b > for the related d-dependent quantity, where B = d± (^JA-KV\q+q~ \/ekT ) . 2.3.2 Other Definitions Fuoss [11,12] objected to the Bjerrum theory on several grounds. He argued that the value R which Bjerrum took as the upper limit of integration (equation 2.156) had only a mathematical significance, not a physical one. He also noted that while this definition of R is not energetically unreasonable, for low dielectric constant solvents R corresponds to ions separated by several solvent molecules. His main objection, however, was that the larger the value of R one chooses, the larger the integral in equation 2.156 becomes until eventually all the ions are counted, giving unphysical (negative) values for the number of free ions. Fuoss [11,12] proposed an alternate definition that an anion would be considered paired to a reference cation if there was no other cation closer to it. If for notational simplicity we introduce the Bjerrum parameter b = lg+g-1 ekTd±' the Fuoss probability for finding an anion paired to a reference cation is bd+ tr , (bd (2.159) G(r)dr = 4:TTp2r exp — — — 47iyj2 / x2 exp [ —— ) dx r Jd± \ x I dr. (2.160) Chapter 2. Theory 62 It can be seen that this function is composed of two terms, one due to ions interacting directly, and the other a screening term. If the screening term is set to zero the Bjerrum result is recovered. It is still necessary to choose a distance R inside which ions will be considered paired, but for this function the integration G{r)dr = 1, (2.161) so Fuoss' main objection to the Bjerrum result is circumvented. From here the calculation of a and K~A then proceeds as for the Bjerrum theory. A second theory of Fuoss [11,13] finds an expression for the ratio of paired to free ions, 1 —a Aird+a (—u(ij)\ — = - f - ^ e x p H H ' <2'162) and uses the d-dependent equation 2.136b evaluated at r = d± to get, 1 - a Airdl ( b \ — = - j i p , « x p ^ T T - 5 - ] , (2.163) which gives a KA consistent to within a constant of that derived by Denison and Ramsey [11,103] using purely thermodynamic arguments. 2.3.3 Ion Pair Definitions for the Present Study. The calculations in this study all yield a probability r2g+_(r)dr as do the theories above. The necessity to define an ion pair remains, however, and the following definitions will be used here. For all models except Model V there are many local minima in g+~(r), so we choose R at the local minimum closest to ion-ion contact and write rR 1 - a = ATps / r2g+_(r)dr, (2.164) J d± which defines ion pairs as those found within the first coordination sphere of radius R. Unlike the Bjerrum theory however, the molecular models in this study give a function Chapter 2. Theory 63 <7+_(r) which has distinct local minima, whereas it is necessary to multiply the Bjerrum g+-(r) by r2 in order to construct a minimum. A possible definition of free (rather than paired) ions can be found by observing the long range behavior of g+-(r) which, if it followed the DH limiting law would be given by equation 2.139 with the exponential decay given by the DH screening parameter K. The results of the calculations here show that the long range behavior is in fact exponential, but with an effective screening parameter somewhat reduced from that of the DH theory. This suggests a new method of calculating the proportion of free ions, a. We write the concentration of free ions as ap2 and assume that the free ions are behaving as predicted by the Debye-Hiickel theory, and the remaining associated ions do not affect this behavior in any way. Within this assumption the effective DH screening parameter « e / / will be given by neff = 47T \q^q-\vap2 (2.165) UkT If we now use equation 2.135b for the DH K we find an alternate route to a, « = ^ T . (2-166) where neff is obtained from the slope of a (linear) plot of ln[r(flf+_(r) — 1)] versus r at long range. We note that there is an inherent advantage in obtaining the proportion of free ions directly, which is that no assumption about the numbers of other clusters of ions (triples, quadruples, etc.) is necessary. This method of finding a and that given by equation 2.164 are not formally equivalent and give different but similar values. Chapter 2. Theory 64 2.4 Computing Considerations. 2.4.1 Method of Numerical Solution. The solution for a given set of input parameters was found by an iterative process. Beginning with an initial estimate for the function c(12)^ the OZ equation (section 2.1.4) is solved for the function ^(12)^. The function 7/(12)^ is then used in the RHNC closure (section 2.1.5) and a new function c(12)^ e u J is found, c(12)L, °4 ,,(12)1, c(l2)Tr, (2.167) where i indicates the i-th iteration. The initial estimate for c(12)^ and the newly found c(12)l^ew are then mixed proportionately, c i + 1(12)^ = (1 - <TK(12)q/3 + <rcijnew {12)^, (2.168) where 0 < a < 1, and the result becomes the initial estimate for the (z -f l)-st iteration. The proportion a of the newly found c(12)^ e u ; which is taken is made larger as the difference between the functions for each successive iteration decreases. Convergence can usually be hastened if the value of cr used for an individual projection is varied separately, according to the amount of change per iteration in that projection. The initial starting point is usually taken as the converged solution for a system with slightly different input parameters. 2.4.2 Choice of Basis Set. The basis sets used for Models I, II and III throughout the calculations are determined by the properties of the models as defined in section 2.1.2, with the only choices being the maximum values of m, n and / which are used. Model III for pure ammonia has been studied extensively [34] with regard to basis set convergence, and it was found that Chapter 2. Theory 65 m,n < 4 and I < 8 gave results for dielectric constants, configurational energies and all projections within a few percent of larger basis sets, while smaller basis sets were clearly inadequate. These maximum values of m,n and I have also been used in similar models for water, and the basis set convergence was found to be adequate [47,59,63,64]. While it was possible to use larger basis sets for Model III, it is beyond available computing resources to make the same calculations for Models I and II. Table 2.3 shows the number of projections Unique Projections cpu time (s) Basis set Model Pure solvent Ionic solution Pure solvent Ionic solution m, n < 3 I < 6 1,11 38 51 10 14 III 20 31 4 7 m, n < 4 I < 8 I,H 90 107 21 26 III 35 48 9 12 m,n < 5 I < 10 I,H 158 179 70f 80T III 56 71 30 36 m, n < 6 l< 12 1,11 344 371 200t 250T III 84 101 60 68 Table 2.3: Computing resources needed for calculations using the indicated basis set. The number of projections is a measure of the amount of storage space needed. for given basis sets for pure solvent and ionic calculations, and the approximate amount of cpu time for one iteration of the OZ and RHNC equations. Typically this has to be done 200 to 400 times for a solution to converge for one set of input parameters. The number of projections are admeasure of the-amount of storage space needed. A rough estimate of the amount of space needed to perform a calculation is space •= 1.5+ (j^j + (^)* Megawords, (2.169) where p is the number of unique projections. The calculations not attempted are marked with a | , and the cpu times for these calculations are estimates assuming that the space required for this many projections is available in main memory. If this amount of space is Chapter 2. Theory 66 not available the disk caching that would be necessary would greatly increase the times given. All cpu times given apply to an Amdahl 5860 mainframe computer. 2.4.3 Handling Long-Ranged Potential Tails. An inherent problem in systems with long ranged potentials is numerically Fourier transforming these functions accurately. The necessity to use FFT techniques restricts the ranges and sampling grids of the discrete points in the functions. For functions which have a lot of structure at short range but also have benign long-ranged parts which are not known analytically, a method of Fourier transforming using a logarithmic grid was proposed by Talman [104]. Although this method seems ideal for the present models, and has been used successfully for weakly screened ionic systems [105], it did not prove useful here despite persistent attempts. As pointed out by Talman there are serious problems -with the use of logarithmic transforms for functions which have discontinuities in themselves or in any of their derivatives. In light of the recent appearances in the literature [105,106] of these logarithmic FFT methods, some explanation of the problems which make them unreliable are in order. In addition, the methods which are used instead for handling the long-ranged tails are given. We require functions of the form (see equation 2.31), where as before ji(kr) denotes a spherical bessel function, and m is an arbitrary integer in the range 0, 1, . . . ,/. If we make the variable transformations r = ep and k = eK we can rewrite equation 2.170 as (2.170) (2.171) Chapter 2. Theory 67 where and [104] 1 roo ^ ( i ) = 2 ^ / - o o e i f P e ( m + 3 / 2 ) P / ( e P ) ^ ' V 8 7 r j = 1 2 j = 1 2 (2.172) cos , ?fj » + sin ( j e (2.173) Here p — I — m, (2.174a) * i ( < ) = I m { l n r ( l / 2 - i i ) } , • (2.174b) and $ 2(i) = tan-1(tanh(7rf/2)). (2.174c) With the logarithmic transform method it is no longer possible to transform a pair of real functions in the same complex FFT process when / is even as it is for the linear transform case. The logarithmic method also requires the execution of two FFT's: one to i-space and then another to K-space. Also, because the number of floating point calculations is increased when other parts of a calculation are performed on a logarithmic grid, computing time will only be saved by using this system if the necessary ranges of r and k are sampled effectively with one-quarter or less points than are needed with the linear grid system. As noted by Talman, there are inaccuracy problems for transforms of high order in I, which become more troublesome as / increases. Simply stated, it is possible to get good accuracy for the transform for either negative or positive values of either logarithmic variable, but not both. A choice of m = / will give the most accurate negative K values and least accurate positive K values. A choice of m = 0 will improve Chapter 2. Theory 68 the high K range at the expense of the low K range. If the function is well behaved (not too steep anywhere) a reasonable compromise can be achieved with a value of m from the middle of the allowed range, or by performing the transform twice and "splicing" the result together at K. = 0. The functions that need to be transformed in this study are very badly behaved, however, and even if the splicing method is used, there is a middle range where the accuracy seems to be of the order of only 3 or 4 digits. Even if this were acceptable the cpu time for these high order transforms would be doubled. The only choice would be to use hat transforms (equations 2.32a,b) on a logarithmic scale and then use equation 2.171 for order I = 1 or 0. However, the most serious problems arise due to the discontinuities in the functions and their derivatives. When any projections from a hard-sphere type system are trans-formed, extremely long-ranged functions arise in the variable t. It can be easily shown that a discontinuity in the function gives rise to a 1/i-dependent tail, but if the hard-sphere discontinuity were the only problem a step function added to make the function continuous would remove the problem in f-space. Similar tails arise from discontinuous derivatives of all orders, and these can not easily be removed analytically. The range of £-space directly determines the size of the logarithmic grid according to Sp = iv/NSt, where Sp and St are the distances between successive points in the variables p and t, re-spectively, and TV is the number of points. Even if the discontinuity in the function itself is removed, the range of t needed to avoid truncation of the tail due to discontinuous derivatives would cause the step in p to be so small that the number of points necessary to sample r-space would be larger than the linear grid system requires. Even if the functions are continuous in all their derivatives as with soft-core type calculations, the oscillatory nature of the functions still seems to give rise to relatively long-ranged i-space tails. If these tails in i-space are simply truncated and the further transformation to fc-space is performed, errors from two sources arise. The truncation of the tail causes an oscillation Chapter 2. Theory 69 at low K with period 2K to be superimposed on the true function values. It has been suggested [106] that this error can be overcome by averaging consecutive points in the result where the oscillation is significant with respect to the function value. This does not remove the error caused by truncating the tail, however, as can be seen in Figure 2.2. The solid line is the absolute error (numerical result minus the analytical result) of the transform of the discontinuous function f(r) = 0 for 0 < r < 1 and f(r) = e _ r elsewhere, after this averaging process has been performed. Clearly the oscillation is superimposed on another error. This is an absolute error, and will be roughly proportional to the size of the discontinuity, which for ion-ion projections in very dilute systems is of the order of 105. The dotted line is the absolute error in the linear transform, and is due to the error in the trapezoidal rule. The other source of error occurs when a function is transformed into f-space and it has not reached its large t limit (zero) in half the sampled range of t. The result of transforming such a function over a finite "window" will no longer be the same as the result of transforming over an infinite range. The conditions of continuity, evenness or oddness, and particularly periodicity which accompany transformations over a finite range imply that the initial function which is transformed is only one period of an infinite periodic function and the result in a window in i-space is similarly just one period of a periodic function. If the resulting transform fits in one window then it accurately represents the infinite function we wished to transform. If not, the numerically enforced periodicity will cause the tail from one window to "leak" into the next and all subsequent (higher t) windows and be superimposed upon the functions there. Similarly, the tails from the previous (lower t) windows are superimposed on the window we have calculated. The real part of the i-space function must also be even, so the window must be symmetric about its center. The tail from this negative frequency section is also superimposed on the function, as well as negative frequency tails from previous and subsequent windows. Chapter 2. Theory 70 Figure 2.2: The absolute error in the transform of the function f(r) = 0 for r < 1 and f(r) = e~T elsewhere. The solid line indicates the error in the logarithmic transform when the oscillations are removed by averaging as suggested in [106], and the dotted line indicates the error in the linear transform. Chapter 2. Theory 72 Figure 2.3 shows one window of the real part of the i-space transform, Re[^(i)], of the simple step function f(r) = 1 for 0 < r < 1 and f(r) = 0 elsewhere. The dashed hne represents the correct analytical result, ^/ly+t2 • '^^ i e a D S C i s s a i s plotted on a logarithmic scale to show that the tail is still significant in the middle of the window. This causes the unwanted negative frequency "mirror image" tail (dotted hne) to be superimposed on the function at low t. The dot-dash hne is the sum of these two pieces. It is still significantly different from the result of the numerical transform of f(r) using Sp = 0.008, pmin — —10.0 and N = 2048, given by the solid hne. This additional error is due to the superimposition of tails from windows at lower and higher t, as it can be seen that the analytical result is not zero as it leaves the window at the right hand side. In order to get the value given by the numerical transform at the center of the window to 9 digits accuracy it is necessary to superimpose the tail contributions from 500 windows. In summary, the logarithmic FFT was examined carefully for many test functions and for many projections from actual model calculations. Much time and effort was spent in an attempt to overcome the unsuitability of this method for functions with discontinuous derivatives of all orders, but to no avail. There are a large class of functions which produce long range tails in i-space which cannot be removed analytically unless the discontinuity in all significant derivatives is known exactly. These problems have not been previously fully appreciated and the analysis of the method given in the literature is clearly over-optimistic. The method should only be used with caution and for very smoothly varying potentials. We now turn to the numerical methods which were actually employed in our calcula-tions. All integral transforms in this study are performed on equally spaced points. The functions which are forward transformed into /c-space are the projections of the function cap(12). As the interparticle distance r becomes large c™^ /3(r) —> — u^^^r) so a subset of the projections to be transformed have tails due to potential terms of the form l / r n Chapter 2. Theory 73 Figure 2.3: The result of the logarithmic transform of a step function into t-space (solid line) compared to the analytical result (dashed line). The dotted line represents the negative frequency "mirror image" tail, and the dash-dot line is the sum of the positive and negative frequency analytical results. Transform of a Step Function t Chapter 2. Theory 75 (see equation 2.13. For n > 2 these potential tails can be subtracted from individual pro-jections leaving a relatively short-ranged remainder. The transform of the long-ranged tail is calculated analytically and then recombined with the numerical transform of the short-ranged remainder. For the ion-ion projections however, the form of the ion-ion potential gives C ^ ) - — asr-^oo, (2.175) r where Bij is a constant and i and j represent (possibly the same) ionic species. These potential tails can not be transformed analytically. For projections with this long range component the function [38,54,55,107] eij(r) = ^ { l - e - a r ) (2.176a) r is subtracted from each ^^-(r), with the value of a sufficiently large to effectively remove the tail. This function is analytically transformed to e<^ = kW^Fy <2;176b» The remaining short range part of CQO°J(T") is transformed numerically, and then recom-bined with 6ij(k). We must also back transform the projections of the function T/(12)a/g, ClU') = h7i*<r) - (2-177) From equation 2.139 we can see that for low concentration the ion-ion projections b^.-(r) have a longest component of the form AB^e^-fr which is shorter ranged than the potential tail of -(r). Thus, if the ionic concentration is not extremely low 7 / ™ ° •(&) can be treated similarly to c™°j(T").and back transformed by subtracting 6{j(k), transforming the remainder numerically and again recombining with the analytical part &ij(r). This method is used for ionic calculations for all models down to a concentration of 2 x 10~2 molar. Chapter 2. Theory 76 At very low concentration (10 - 4 to 10~5 molar) the longest ranged contribution to /ioo°ij(r) (the exponential function given above) has become significantly longer than at higher concentrations. The properties of FFT's require that all projections for a given calculation have the same number of points, so if this component of ^o(Wj( r ) a n d there-fore also T/OS°J(7") were not removed and treated analytically, the demands on computing resources would become unmanageable. In order to obtain these extremely low concen-tration calculations for Models IV (McMillan-Mayer) and V (primitive), the following method was used. The ion-ion functions C o o ° i j ( r ) a r e transformed from r-space to fc-space using 9ij(r), but the inverse transform of the functions 7/So°»j(T') c a n n ° t be performed this way because the r-space result would still be too long. As the concentration dimin-ishes the longest contribution to (r), which has the exponentiated form as given above, becomes longer and longer ranged as the damping due to the exponential factor decreases. To remove both these long components from T]g°°a^(r) it is necessary to use the function •^•(r) = ^ ( 1 - e-~) + ( 1 ~ A ) B i j ( l - e-"), (2.178a) r r where B{j is as above, and which has the transform Ax2 (1 - A)ct2 (2.178b) k2 [k2 + x2 (k2+a2) where A = exdi> / ( l -+- xdij), and it is expected that x should approach the Debye K as the concentration is decreased. The x values used for calculations using these corrections were chosen to ensure self-consistency in the charge neutrality relation between the Kirkwood and (see equation 2.120a). In all cases the x value required was within 1% of the Debye K, but had to be specified to nine figures to ensure the Kirkwood G's were consistent to six figures. Reliable results could be obtained using this method at concentrations as low as p2 — 1 x 10 - 6 . Chapter 3 Results and Discussion 3.1 Results for Models of Pure Ammonia The results of the RHNC/SCMF calculations for pure ammonia Models I, II and III, as defined in section 2.1.2 are given here. These models are defined in terms of various input parameters. The values of the parameters used for all three models are given in Table 3.4, as well as their sources. The particle diameter has been chosen to parameter value source T -35°C chosen d 3.2 A chosen P 0.68400 cm,-3 [108] Ql •1.47D [109] Ql -3.35 [110,111] Ql -1.786 x 10~34esii cm3 [112] Ql -2.154 x 10-34esu cm3 [112] Ql 5.187 x lfr^esu cm4 [112] Ql -2.292 x 10~42esu cm4 [112] an 2.412 A3 [113] a j . 2.124 A3 [113} Table 3.4: Input parameters for polarizable pure ammonia models. 77 Chapter 3. Results and Discussion 78 be consistent with the neutron scattering data of Narten [40], from which it can be seen O 0 that the repulsive core in real ammonia takes effect over a range of 3.1 A to 3.3 A- A o few calculations with a diameter of 3.1 A were done and the changes which resulted in the thermodynamic quantities and dielectric constant were 3 — 5%. Where possible, experimental data were used for the input parameters, but the values of the octupole and hexadecapole come from quantum calculations. The effective dipoles and renormalized polarizabiHties from the SCMF theory (equa-tions 2.86 and 2.87) are given for each model in Table 3.5. It is clear that the octupole Model I Model II Model III 2.10 2.07 1.98 m' 2.03 2.00 1.90 a'\\ 3.324 3.282 3.124 2.801 2.771 2.657 Table 3.5: SCMF effective dipoles and renormalized polarizabiHties. The dipoles are in Debyes and the polarizabilities in cubic angstroms. and hexadecapole act to reinforce the effective dipole, and it is this increased dipole that precipitates most of the differences between the models. The few thermodynamic quantities which can be compared with experiments are given in Table 3.6. The total configurational energy from Model I compares quite favorably with the experimental re-sult, whereas the Model III results is not as good. The opposite is true of the isothermal compressibilities, with the Model III result within about 2% of the experimental value. The results for all models for isothermal compressibility factor, Z, are clearly unphysical. This quantity is extremely sensitive to the model in this type of calculation [59], with very small changes in the model causing large differences in the calculted value. The Chapter 3. Results and Discussion 79 Model. I Model II Model III Experiment e 22.7 21.7 18.4 21.6 [41] -UT 10.66 10.06 7.01 10.94 [114] - u D D 3.14 3.02 2.67 -— UDQ 3.41 3.38 3.11 --UQQ 1.18 1.19 1.23 -z -1.32 -1.01 1.50 1.21xl0- 3 [108] XT 0.134 0.128 0.088 0.090 [114] Coordination number 10.25 10.28 11.05 12f [40] Table 3.6: Thermodynamic quantities, dielectric constants and coordination numbers for pure ammonia at equilibrium vapor pressure. The energies are in units of kT per particle, the coordination numbers in numbers of particles and the isothermal compressibilities in cm 3 mol - 1 . All quantities were found for ammonia at — 35°C except f which was found at 4°C coordination numbers compare reasonably well, but it should be noted that the exper-imental result is an estimate from neutron scattering data for a different temperature (4°C). The radial distribution functions for pure ammonia for Models I, II and III are given in Figure 3.4. The most significant difference between the models with an octupo-lar interaction (I and II) and Model III which has no octupole is due to the difference in SCMF effective dipole, (me in Table 3.5), which results in larger contact values for the octupolar models. This in turn leads to sharper second (particle separated) peaks, and the shifting of all maxima and minima to slightly lower radial positions. This is consistent with the increased interaction of larger dipoles, shown by the dipole-dipole energies in Table 3.6. The further addition of the hexadecapole (Model I) gives a radial distribution function which is essentially unchanged from Model II. The slightly higher me value results in the same general trends as when the octupole was added. The small Chapter 3. Results and Discussion 80 Figure 3.4: goo°{r) for pure ammonia. The solid, dashed and dotted lines represent Models I, II and III, respectively. Chapter 3. Results and Discussion 82 effect of the hexadecapole is partly due to its relatively small value, and partly due to the short range of the additional interactions. Figure 3.5 shows the average of the cosine of the angle between the dipoles of particles 1 and 2, as given by equation 2.113. It is clear from the contact values that there is a predominance of structures with acute angles between the dipoles. The effect of the negative quadrupole in ammonia is to reinforce the negative end of the dipole and almost completely neutralize the positive end, giving the effect of redistributing the positive charge towards the "equator" of the molecule. The tendency of dipoles to align is offset by the tendency for quadrupoles to form orthogo-nal T-shaped structures, resulting in a compromise of acute angles. If we approximate (8), the average angle between dipoles, as arccos(cos 9), the average contact angles for Models I, II and III are about 80°, 82° and 84°, respectively. As expected, a slightly larger effective dipole induces the tendency toward dipole alignment. A likely candidate to contribute to this average structure is the lone pair of one molecule getting as close as possible to a proton in another molecule, with its dipole almost lining up with the NH-bond. This would yield the slightly flattened tetrahedral structure indicated by the contact angles. Closely related to the average dipole-dipole angle is the function hgl°(r), given in Figure 3.6, which can be seen as the average angle between dipoles at a given sep-aration r, multiplied by the probability, relative to the bulk density, of finding a molecule there. The relationship between this function and the dielectric constant (see Table 3.6) is given by equations 2.66a-c. The dielectric constant is one of the few quantities that can be compared to experimental values, and it compares very favorably. A recent study shows that Model III gives good agreement with experiments over a range of tempera-tures [64,92], but consistently underestimates the dielectric constant by a small amount. The inclusion of the higher multipoles increases the calculated results to make up this slight underestimation. Computer simulations have shown [39] that dielectric constants from similar RHNC/SCMF calculations for ammonia are quite accurate. It should be Chapter 3. Results and Discussion 83 Figure 3.5: -A2) v s - r f°r P u r e ammonia. The solid, dashed and dotted lines are as in Figure 3.4. Chapter 3. Results and Discussion 85 Figure 3.6: hl0°(r) for pure ammonia. The solid, dashed and dotted lines are as in Figure 3.4. Chapter 3. Results and Discussion 87 noted that the dielectric constant calculated for a non-polarizable model with the same parameters as Model III is 9.02 [64,92], which seriously underestimates the experimental values. The dipole-dipole energy is related by equation 2.76 to hH2(r) which is shown in Figure 3.7. The higher effective dipoles in models containing the higher multipoles are consistent with the higher contact values and higher dipole-dipole energies, and again there is only a small difference between Models I and II. The peak at unit radial sep-aration for Models I and II is a consequence of the increased alignment in molecule separated dipoles as well as the size of the effective dipoles themselves. The dipole-quadrupole and quadrupole-quadrupole energies (Table 3.6) are calculated from the func-tions hH3(r) and hH4(r), which are given in Figures 3.8 and 3.9. The same features due to the increased dipole size and alignment are present in both these graphs. Since the quadrupole-quadrupole energy does not depend directly on the magnitude of the dipole, the differences between the models must be due to correlation effects of dipoles only, so the change in energy is much smaller than for the dipole-dipole term. The ab-solute value of the quadrupole-quadrupole energy decreases (note the abscissa of Figure 3.9 is inverted) because the aligning of dipoles is unfavorable for quadrupole-quadrupole interactions. Figure 3.10 shows hH6(r), from which is calculated the linear part of the octupole-octupole energy for Models I and II. This function is present in Model III but does not contribute to the total energy. The function displays the range where the octupole-octupole interaction is significant, showing that the effect of the octupole has disappeared outside a separation of about one particle diameter. This function, as well as those in Figures 3.5 to 3.9 all show that more detailed structure exists in Models I and II than in Model III, consistent with the change to lower (i.e. C3, ,) symmetry. Chapter 3. Results and Discussion 88 Figure 3.7: hH2(r) for pure ammonia. The solid, dashed and dotted lines are as in Figure 3.4. 0.15-1 0 .12H 0 .09H 0.06H 0 .03H 0 .00-hll2(r) 8 . 0 - 1 4 . 0 0 . 0 0 . 0 0 . 8 i — i — i — n ~ 1.6 2 . 4 (r-d.)/i. 3.2 0 .15 0.3 4 . 0 4 . 8 GO CO Chapter 3. Results and Discussion 90 Figure 3.8: hH3(r) for pure ammonia. The sohd, dashed and dotted-lines .are-as in Figure 3.4. 0 . 0 0 . 8 1.6 2 . 4 3 . 2 4 . 0 4 . 8 (r - d.)/d. Chapter 3. Results and Discussion 92 Figure 3.9: h™(r) for pure ammonia. The solid, dashed and dotted lines are as in Figure 3.4. - 0 . 0 8 - 1 i I0xhll4{r) _ 6 < 0 - 0 . 0 6 H - 0 . 0 4 H - 0 . 0 2 H - 0 . 0 0 -0 . 0 2 0 . 0 0 . 8 0 . 1 5 0 . 3 1 . 6 2 . 4 (r-dt)/dt i—i—r 3 . 2 "1 1 1 4 . 0 4 . 8 CD CO Chapter 3. Results and Discussion 94 Figure 3.10: hH6(r) for pure ammonia. The solid, dashed and dotted lines are as in Figure 3.4. Chapter 3. Results and Discussion 96 3.2 Results for Macrospheres at Infinite Dilution. 3.2.1 Reference Hard-Sphere Bridge Functions for Large Particles. The reference bridge function &R(r) (see equation 2.25) used for the calculations of pure ammonia in the previous section was calculated from the exact hard-sphere radial distribution function which in turn was calculated from the Verlet-Weis [72] fit to data from Monte Carlo calculations. Obtaining bn(r) in this manner requires Fourier inver-sion and solution of the OZ equation for hard spheres. Unfortunately, extremely subtle differences in <7i?(r) lead to significant differences in 6R(7-), SO the radial distribution func-tions must be known with enough accuracy to determine acceptable bridge functions. The Lee-Levesque [73] radial distribution functions for mixtures of hard spheres differing in diameter by a factor of 2 or less are sufficiently accurate (3 significant figures every-where) to give bridge functions which should differ from exact results by not more than a few percent. Subsequent sections will examine results for ammonia near the surface of macrospheres of up to 20^, where ds is the solvent diameter, so it is first necessary to determine whether the Lee-Levesque reference functions are reliable for mixtures of particles which differ in diameter by this amount, and to compare results with those using other choices of reference functions. Figure 3.11 shows the radial distribution functions for several hard sphere systems, all with a density of p = 0.793d73- The solid line is an RHNC calculation for hard spheres around a single macrosphere with diameter equal to 20 solvent diameters. The fci?(r) for this calculation was obtained directly [116] from a diagrammatic expansion of the bridge function in powers of p, keeping terms up to order p3. The integrals in this expansion were calculated using Monte Carlo techniques where analytical results were not possible. This result will be referred to below as the "direct method". The dashed line is the radial distribution function obtained directly from the Lee-Levesque [73] fit to Monte Chapter 3. Results and Discussion 97 Figure 3.11: gii(r) for hard spheres. The solid hne represents the gR(r) corresponding to the direct method of calculating bR(r) for hard spheres around a single 20*^ microsphere. The dashed hne is the Lee-Levesque [73] fit to Monte Carlo data for the same system. The dotted hne represents the gR(r) fit to Monte Carlo calculations by Henderson [115] for hard spheres by a flat wall, and the dash-dot hne is an HNC calculation (bR(r) = 0). All calculations are for a density of p = 0.793d~3. Chapter 3. Results and Discussion 99 Carlo data for the same system of hard spheres around a single macrosphere. This result will be referred to below as the "Lee-Levesque" result. The dotted hne is taken from the Henderson [115] fit to Monte Carlo calculations for hard spheres of the same density next to a fiat wall, and is referred to below as the "fiat wall" result. In order to compare the flat wall and macrosphere data the effects of curvature of the surface must be small compared to other differences. The contact values for the Lee-Levesque gn^rfs for 10ds, 20ds and Z0de diameter macrospheres are 6.45, 6.74 and 6.83, respectively, showing the curvature effects are minimal in macrospheres of this size. However, these contact values do not appear to be converging to the flat wall contact value of 7.58 which is known to a very good approximation, indicating the difference between these two results is not primarily an effect of curvature. The dash-dot hne is an HNC (bR.(r) = 0) calculation. The negative of the bridge functions ( — bR(r)) can be considered as additional "poten-tials" added to the approximate system to improve the approximation. Comparison of the three reference system radial distribution functions calculated with non-zero bridge functions to the HNC result shows the effect these corrective potentials have on the HNC result. In general for hard spheres the HNC result gives contact values which are too large, and all three bridge functions considered correct for this. The flat wall result is probably close to exact at contact, and the direct method of calculating bn(r), although for a 20d macrosphere, very closely approximates this value. The Lee-Levesque result seems to over-correct the HNC solution. The HNC result also in general gives a first minimum which is too shallow, for which the Lee-Levesque and flat wall functions com-pensate, but the direct method does not. The bridge functions also have an effect on the period of oscillation of the gii(r) functions, the flat wall extending the period, and the Lee-Levesque reducing it. It is not clear which of these adjustments is more reliable, al-though the fiat wall fit imposes a constant period on gii(r) which is probably not correct, and previously published Monte Carlo data [117] seem to indicate the second maximum Chapter 3. Results and Discussion 100 to be at a slightly lower r value than is given by the HNC gR.(r). In general, however, all the radial distribution functions are remarkably similar, especially when one considers the bridge functions from which they are calculated. Figure 3.12 shows the bridge functions corresponding to the gR(r) functions in Figure 3.11. The calculation of the bridge function for the flat wall gB.(r) was done by assuming the curvature effects near a 20d„ macrosphere were minimal and solving the OZ equation as if this gn(r) was that of a macrosphere. The bridge function which results is not exactly the same as that which would be calculated if the OZ equation were solved for a flat wall system, but this is a very good approximation because the bridge functions calculated this way for 20<fs and 30d„ macrospheres are essentially identical everywhere. Therefore, in order to distinguish between the bridge functions and identify them with their sources, this particular bridge function will still be referred to as the flat wall result. It should be noted that the scale of this plot is similar to that used for the radial distribution functions, so the differences in the bridge functions are very significant. The flat wall and direct method, which give the same contact value for gii(r) have contact values in t>R{r) which are not only quite different in magnitude, but also are of opposite sign. If we regard —fcfl(r) as a corrective potential, the direct method adds a repulsion at contact, while the flat wall adds an attraction, yet both result in the same contact value in gji(r). It is clear that the high gii(r) contact peak in the HNC calculation can be corrected in several ways. The direct method simply adds a repulsion at contact, forcing the hard spheres away from the surface. The flat wall and Lee-Levesque results add an attractive potential well between contact and the first minimum which draws spheres into this area from neighboring regions, decreasing the minimum in gR.(r) as well as reducing the peak height at contact. It is also clear that, depending on viewpoint, either gii(r) is extremely insensitive to changes in bfi(r), or bji(r) is extremely sensitive to small changes in gR.(r). The difference in maxima in the Lee-Levesque and flat wall bridge functions is greater Chapter 3. Results and Discussion 101 Figure 3.12: bn(r) for hard spheres. The hne types refer to the same approximations as in figure 3.11. 4.5-. ' 3 . 0 H 1.5H 0 . 0 -- 3 . 0 bR(r) / \ / \ \ \ \ V \ \ \ \ ~ ~ 1 — I — I — I — I — I — I — I — I — I — I — I 0 . 0 0 . 4 0 . 8 1 . 2 1 . 6 2 . 0 2 . (r - dms)/da Chapter 3. Results and Discussion 103 than a factor of 2, whereas the corresponding differences at similar values of r in the radial distribution functions are almost negligible. The effect of the different reference approximations on the radial distribution function from a full Model I calculation is given in Figure 3.13. Unfortunately, unlike the rela-tively small changes in which arise from the different choices of the reference bR(r) function, the changes are more pronounced when Model I ammonia is used as the solvent. While all choices for bR(r) reduce the contact values in 5foo0Tn«(r), the difference between the calculation using the direct method reference bR(r) and those using the other two non-zero bridge functions is quite significant. It should be noted that while the flat wall hard-sphere gR(r) had the same contact value as the direct method hard-sphere gR(r), when the the corresponding bR(r) functions are used as reference functions the resulting contact values in goo^fr) are not the same. In fact, the result using the flat wall bR(r) is now paired with the result from using the Lee-Levesque bR(r), although their behavior just outside the wall is quite different, with the latter showing a density maximum at a separation of about 0.1ds away from the surface of the macrosphere. Again it should be emphasized that this is not a curvature effect. It is not clear which reference bridge function is the best, but they all seem to agree that the HNC calculation gives a func-tion which is too steep at contact. The sensitivity to reference bridge function indicates that even if the exact reference bridge functions bR(r) were known, it is probable that the corresponding distribution functions would differ significantly from those calculated using the correct angular dependent bridge function 6(12). Some structural effects of the model are shown in Figures 3.14 to 3.16. These functions .are important components in the angular probability distributions described in section 2.1.8, and their relationship to the structure of ammonia at a surface will be discussed below, so the dependence of these quantities on the choice of reference bridge function is of importance here. The differences between the functions in each figure is solely due to Chapter 3. Results and Discussion 104 Figure 3.13: Model I goo°ms(r) functions using different bR(r) approximations. The line types refer to the same approximations as in Figure 3.11. Chapter 3. Results and Discussion 106 Figure 3.14: Model I h^^^r) functions using different bR(r) approximations. The hne types refer to the same approximations as in Figure 3.11. 0.35-n 0 . 2 5 H 0.15H 0.05-^ -0 .05H -0.15--3.0-n -1.5H > i H i • i \ / / i W 0 . 0 • r/Ai 0 . 0 0 . 8 1.6 2 . 4 ~ dm8)/da 3 . 2 4 . 0 Chapter 3. Results and Discussion 108 Figure 3.15: Model I fcgj|?m,(r) functions using different bji(r) approximations. The line types refer to the same approximations as in Figure 3.11. 0 . 1 5 - 1 0.10 0.05H h022 (r) - 0 . 0 0 - 0 . 0 5 H -0 - 1 . 0 -(r - dma)/da Chapter 3. Results and Discussion 110 Figure 3.16: Model I h°^ms(r) functions using different fcfl(r) approximations. The line types refer to the same approxmatipns systems as in Figure 3.11. 0.2<H 0 . 1 6 H 0 . 0 8 H h 0 3 3 (r) 0 .50-L - 0 . 0 0 H -0.08H - 0 . 1 & 1.0 1.5 2 . 0 Chapter 3. Results and Discussion 112 the different choices of reference function. The angular dependent potential of the solvent seems to emphasize the differences between these choices significantly, even though the reference functions themselves have no angular dependence. Small peaks and shoulders which appear in these figures indicate that small differences in the curvature of g2o,°m«(T') couple strongly with the projections reflecting angular correlations by emphasizing and de-emphasizing regions where the different multipoles have different amounts of influence. Some general trends do appear in all these figures, however. The contact peak for the HNC calculations are all concave upward, whereas all the calculations with non-zero reference functions are concave downward. This strongly suggest that the HNC results predict too great a density directly in contact with the surface, and shows the need for the correction that the bridge functions afford. The relative size of the peaks, or the amount of structure present, reflects the relative size and structure of the bridge reference functions themselves. The direct method bridge function is identically zero beyond a separation of 1 solvent diameter and is extremely benign. This gives rise to a similar lack of structure in the full solvent model, both in absolute size and in range. The Lee-Levesque reference system gives the largest and most structured bR(r), which is reflected throughout these figures. The difference between these extremes is significant. The result of the full Model I calculation using the flat wall reference bridge function seems to agree with the equivalent result using the Lee-Levesque bridge function, having similar but less pronounced structure everywhere, which suggests that the direct method does not adequately describe the true hard-sphere bridge diagrams. The final comparison between Model I results using different reference functions is given by.Figures 3.17 to 3.19. These figures display the probability density of finding a solvent molecule with a given angle between its dipole vector and a vector orthogonal to the surface, p^ (#) (see equation 2.105). The principal differences between the values of Pn(6) for each of the calculations is one of degree rather than kind. There are no peaks Chapter 3. Results and Discussion 113 Figure 3.17: Probability density for the angle between dipole and surface normal at contact using different 6^(r) approximations. The line types refer to the same approxi-mations as in Figure 3.11. Chapter 3. Results and Discussion 115 Figure 3.18: Probability density for the angle between dipole and surface normal at a separation of 0.5 solvent diameters using different bn(r) approximations. The line types refer to the same approximations as in Figure 3.11. Chapter 3. Results and Discussion 117 Figure 3.19: Probability density for the angle between dipole and surface normal at a separation of 1.0 solvent diameters using different bji(r) aproximations. The hne types refer to the same approximations as in Figure 3.11. Chapter 3. Results and Discussion 119 distinct to a single reference system, and all maxima and minima are at the same value of 8. This is due to the fact that each function has been normalized by its own fl,So°(7')) removing most of the individual distinctions. The results which had the most structure in the previous figures have the greatest amplitudes here. The differences are most obvious in Figure 3.19 where the direct method bn(r) has become identically zero. This bji(r) has not only not added any structure, but it has actually reduced the structure found in the H N C calculation. Qualitatively, the most probable angular distribution of solvent next to a wall is unambigupus and independent of the reference function is chosen. It seems clear that guaranteeing a reasonably accurate radial distribution function in the l imit as the potential approaches that for hard spheres does not necessarily guaran-tee accurate pair distribution functions for the present model. Fortunately, despite this problem unambiguous angular probabilities can be calculated using any of the reference bridge functions examined here. We wil l use Lee-Levesque [73] reference functions for all subsequent surface calculations, simply for consistency with all the other calculations in this study. This is a reasonable choice because calculations using the flat wall reference function do not differ enormously from those using the Lee-Levesque reference function, and the direct result is identically zero for r > 2ds, which is probably incorrect. Fur-thermore, the range of structure i n results using the Lee-Levesque or flat wall reference function have a similar range of structure as results from Monte Carlo calculations [117]. Chapter 3. Results and Discussion 120 3.2.2 Model Dependence in Large Particle calculations. One of the quantities of interest when examining the structure of ammonia near the surface of a macrosphere is the probability density of the NH-bond angle with respect to a surface normal vector. Since the arrangement of NH-bonds in ammonia has C3,, symme-try, and Model III has linear symmetry, it is expected that there should be significant dif-ferences between the structure predicted by Models I and II and that predicted by Model III. Figure 3.20 gives g^ma(r), the radial distribution function for solvent molecules near the surface of a 20da diameter macrosphere. Models I and II are essentially identical everywhere, and Model III differs only at contact, showing a higher contact value than the other two models. The contact values for the three models are qualitatively inversely proportional to the size of their SCMF effective dipoles, as we would expect. This is the opposite of the result given by the solvent-solvent radial distributions (see Figure 3.4). The stronger dipole-dipole interactions induce stronger solvent-solvent interactions, effec-tively pulling the ammonia molecules away from the macrosphere surface. The function ^oo-msC71), shown in Figure 3.21, also has very little variation between Models I and II, and the region very close to the surface is again the only place where Model III differs sig-nificantly. This effect is also induced by the weaker solvent-solvent interaction in Model III, as the effects in 5oo°m«(T') a r e coupled to all the projections. The projection ^oo2m«(r) given in Figure 3.22 displays the same trend. The projection hQ0i^ms(r) (Figure 3.23) does not even have a significantly different contact value. These similarities in the projections of gm8(12), as expected, give a similar probability distributions for the angle between the solvent dipole and a surface normal vector, as given in Figure 3.24. Even though Model III doesn't have the same sjTnmetry, the angle of the molecular dipole is unaffected, with both the most probable orientation and the fluctuations about this most probable angle being equivalent for all three models. Model III can be interpreted as a full C3,, model Chapter 3. Results and Discussion 121 Figure 3.20: Model dependence of goc?me(r) f° r ammonia near the surface of 20cZ, macro-spheres. The solid line represents Model I, the dashed line represents Model II and the dotted line represents Model III. 122 Chapter 3. Results and Discussion 123 Figure 3.21: Model dependence of /ioo?m«(r) f° r ammonia near the surface of 20d„ macro-spheres. The lines are the same as in Figure 3.20. 0 . 5 0 -0 . 3 5 H 0 . 2 0 -0 . 0 5 H - 0 . 1 0 H 0 . 0 i -1 .5H i — i — r 0 . 0 0 . 1 5 0 . 3 - 0 . 2 5 -0 . 0 0 . 8 1.6 2 . 4 3 . 2 {r - dms)/dt 4 . 0 4 . 8 t o Chapter 3. Results and Discussion 125 Figure 3.22: Model dependence of ^oofm*(r) for ammonia near the surface of 20ds macro-spheres. The lines are the same as in Figure 3.20. 0 .30 -1 0 . 2 2 -0.14H 0.06H -0.02H -0.10-- 0 . 5 - 1 0 . 0 0 . 8 -1 . O H I 1 0 . 0 0 . 1 5 ~i—i—i—i—i—r 1.6 2 . 4 3 . 2 {r - dma)/da 0 . 3 4 . 0 4 . 8 to 0 5 Chapter 3. Results and Discussion 127 Figure 3.23: Model dependence of /ioofm.C7') f° r ammonia near the surface of 20^ macro-spheres. The lines are the same as in Figure 3.20. 0 . 1 6 -0.08H -o.ooH -0.08H - 0 . 1 6 -0.0 0.8 1 I 1 1 1 1 1.6 2 . 4 3 . 2 {r - dma)/da 4 . 0 4 . 8 to OO Chapter 3. Results and Discussion 129 Figure 3.24: Model dependence ofpM(0) for ammonia at the surface of 20ds macrospheres. The lines are the same as in Figure 3.20. Chapter 3. Results and Discussion 131 with the potential averaged over the Euler angle representing the rotation of a molecule about its dipole axis. Clearly this averaging has little effect on the distribution of dipole angles throughout the solvent, nor does it have much effect on the radial distribution of particles near a surface. As would be expected, the distribution of the angle between the NH-bond of an ammonia molecule and the surface normal vector is quite different in Model III than that obtained for Models I and II (Figure 3.25). It should be emphasized that Model III has linear symmetry, and as such has no preferred angle about its dipole axis for the location of an NH-bond. But, for purposes of comparison and for the calculation of PNH{@) o nty> an NH-bond is defined to be at rotations of 0, 27r/3 and 47r/3 radians about the dipole axis in the molecular coordinate frame. The angle between an NH-bond and the dipole is chosen to be consistent with the bond angle given by the quantum calculation from which the octupole and hexadecapole for Model I were taken [112]. The fact that .there is no preferred angle for an NH-bond about the dipole in Model III does not mean there is no preferred angle between the NH-bond and a surface normal, because the rigid structure of the molecule and the dipole distribution enforce an angular distribution on the N H -bond direction. This is extremely useful, because the Model III line of Figure 3.25 can be interpreted as the angular distribution of NH-bonds given the imposed dipole distribution and no hindrance to free rotation about the molecular dipole axis. This makes it possible to see the effect of the hindrance to free rotation invoked by the addition of the octupole to the model. The intensification and narrowing of peaks indicates that the rotation about the dipole is quite severly hindered, and that the structure for Models I and II has a tendency towards a structure near a surface that is more defined than the structure observed in the bulk solvent. Chapter 3. Results and Discussion 132 Figure 3.25: Model dependence of PNH(O) for ammonia at the surface of 20ds macro-spheres. The lines are the same as in Figure 3.20. Chapter 3. Results and Discussion 134 3.2.3 Dependence of Radial Distribution on Macrosphere Size. The dependence of the radial distribution function goo?ms(r) on the size of the macro-sphere is compared in Figure 3.26. In particular, the changes in local density near the surface of the sphere, and the intensity and range of structure, as affected by the size of the macrosphere, are of most interest. As has been seen above, the intensity and range of structure in the projections of total distribution functions are significantly affected by the choice of reference bridge functions. However, we find that the type of changes which occur as the macrosphere size is increased are common to all choices of bridge function. When the Lee-Levesque reference bridge function is used, the result for the smallest diameter sphere (lde) is concave up in the region close to contact with its maximum directly at contact. When the sphere has reached a diameter of 5d„ the first maximum has already moved away from contact and the function has become concave downward, with only small changes in curvature and peak position for larger spheres. Thus the greatest changes in the general shape occur as the spheres are increased from lds to bdt in diameter. The intensity of the structure increases beyond this size sphere, however, with roughly equal increases in the local density near contact for each 5ds increase in macrosphere size. While the absolute shape of the curves from calculations using other reference bridge functions would differ somewhat, the types of changes that occur with sphere size are quite similar. The size of the oscillations in <7o(Hn« ir) fr°m the second minimum and beyond to larger values of r also follow the trend of increased intensity following the size of the macrospheres: In the intermediate range the 20^, functions no longer show any increased amplitude in structure, with the first minimum and second maximum having a lower amplitude than those for the 15^ spheres. The position of each relative minimum and maximum also moves closer to contact as the macrosphere is enlarged, but once again the 20de macrospheres no longer follow this trend. While the Chapter 3. Results and Discussion 135 Figure 3.26: Dependence of goo°ms(r.) on hard sphere size. The solid, short-dashed, dot-ted, dot-dashed and long-dashed lines represent macrospheres in Model I ammonia with diameters equal to lds, 5da, 10ds, 15da, and 20da, respectively. Chapter 3. Results and Discussion 137 exact position and height of these middle ranged peaks is highly dependent on bridge function, the breaking of these trends for 20tis macrospheres indicates that the changes in local density as a result of the changes in surface curvature are becoming small, and will not alter significantly when the macrosphere is increased in size beyond 20d„. This is consistent with curvature effects no longer being significant for any of the reference radial distribution functions examined above. In Figure 3.27, the function h^}me(r) for the same calculations displays similar tendencies of position and intensity of structure. The first peak beyond contact increases with sphere size, and continues to do so for the largest spheres considered. This indicates the angular distribution near the surface is slightly more sensitive to surface curvature than is the local density. 3.2.4 Dependence of Angular Structure on Macrosphere Size. The angular structure of the solvent near the surface of a macrosphere is also enor-mously affected by the size of the sphere. It is most easily represented by plotting the probability distribution pM(#,r") or PNH{Q>T) a s a function of both the angle 6 between the vector in the molecular coordinate frame (the dipole vector p or the NH-bond) and a surface normal vector, and also of the separation distance r between a solvent molecule and the macrosphere. In all such figures the angular distribution at contact is that fur-thest into the page. All angular figures are symmetric about 180°, which arises from the definition of an angle between a vector in the molecular coordinate frame and a surface normal vector rather than molecular symmetry. Figures 3.28, 3.29, 3.30 and 3.31 rep-resent Pfi(@,r). for hard spheres of diameters lds, 3da, 5de and 20de, respectively. The dipole distribution is relatively unstructured for lds spheres (Figure 3.28), with about equal tendencies for the dipoles to be pointing at the surface and approximately parallel to the surface. As a solvent molecule comes away from the surface the dominant angles are relatively unchanged, with the molecular dipoles which were parallel to the surface Chapter 3. Results and Discussion 138 Figure 3.27: Dependence of h^mt(r) on hard sphere size. The lines represent the same sphere sizes as in Figure 3.26. 0 .35-n 0.25H 0 . 1 5 - J 0 . 0 5 H - 0 . 0 5 H - 3 . 0 0 - 1 A / \ / \ w ft // \l I // \v I; //, -1.50H 0 . 0 0 w /// - ° - 1 5 n — i — i — i — i — i — i — i — i — i — i — r 0 . 0 0 . 8 1.6 2 . 4 3 . 2 4 . 0 ( r - d m s ) / d a Chapter 3. Results and Discussion 140 Figure 3.28: p^(6,r) for lds diameter hard spheres in Model I N H 3 . 141 Chapter 3. Results and Discussion. 142 turning slightly away from it. The structure is only significant in the region within ld„ of the surface. For 3<£„ diameter hard spheres (Figure 3.29), the same structure as in the lds case is apparent, except a new feature has appeared close to the surface. The solvent molecules in contact with the surface now have a tendency to point at the surface with an angle away from the normal vector of about 120°. Slight undulations in p^(0,r) are beginning to appear outside the r = 1.0da range and these new peaks occur at the same angles as the new structure at contact. When the diameter of the hard sphere is increased to 5da (Figure 3.30), the effects which appeared at lower diameters are more pronounced with the structure which filled the region from contact to lds for the smallest spheres now restricted to the range 0.2ds <r< 0.7d„. Figure 3.31 shows the same result for a 20da diameter macrosphere. All the peaks have intensified and narrowed, but their positions are essentially the same as those in Figure 3.30 for 5da spheres. The dipole orientational structure has a marked dependence on the size of the macrosphere, but be-yond a size of about 3da the changes are in degree rather than kind. No additional peaks appear or disappear beyond this size, but when the sphere has become large compared to the solvent the structure has become quite definite and three distinct layers in the angular structure are apparent. Between the surface and a separation of about 0.2da there is a predominance of dipoles pointing at an angle of 120° to the surface normal. Between r = 0.2da and about 0.7da there is a predominance of dipoles pointing at the surface and in a direction about 60° away from the surface normal. Beyond r = 0.7da to about r = 1.3d,'the dominant orientations are again similar to those at the surface with the additional tendency for dipoles to point directly away from the surface. The dominant angles in each layer are approximately at a tetrahedral angle to those in the previous layer, and the amplitude of the structure in each layer decays with distance from contact. This is entirely consistent with the multipole structure of Model I. The charge distribution around an N H 3 molecule has a strongly negative region around the Chapter 3. Results and Discussion Figure 3.29: p^(6,r) for 3d, diameter hard spheres in Model I N H 3 . 144 Chapter 3. Results and Discussion 145 Figure 3.30: p^(0,r) for 5ds diameter hard spheres in Model I NH : 146 Chapter 3. Results and Discussion 147 Figure 3.31: p^(8,r) for 20ds diameter macrospheres in Model I NH: 148 Chapter 3. Results and Discussion 149 electron lone pair, and the corresponding positively charged region is located around the "equator", concentrated around the three hydrogen atoms at slightly flattened tetrahe-dral positions. In order for the negatively charged region of one ammonia molecule to be in close proximity to the positively charged region of another, the successive layers should have dipole vectors at angles roughly tetrahedral (100°-120°) to each other. The angular distribution P N H ( @ , r ) of NH-bonds is given for lds hard spheres and 20ds macrospheres in Figures 3.32 and 3.33, respectively. The onset of angular structure is very similar to that of dipoles. The same layering of dominant angles and the surface structure appears when the spheres are 3ds or less in diameter, and intensifies as the spheres grow to 20ds in diameter. The same approximately tetrahedral difference between the angles in the successive layers is observed, except that the layers are oppositely placed in the r-direction with an NH-bond pointing directly at the surface at contact, and directly away from it at separations of between 0.2d„ and 0.7ds. As was noted above, it is expected that in rigid models, a dipole angular distribution will have associated with it a natural distribution in the NH-bonds, due to the fixed angle between an NH-bond and a dipole. As was seen above, this NH-bond distribution is intensified when the free rotation about the dipole is hindered in non-linear models such as Model I. The intensification of structure is due to the NH-bonds in one layer favoring a direction towards the negative end of a dipole in surrounding layers. This "hydrogen bonded" structure is similar to the lattice found in frozen ammonia. The most probable orientation for an ammonia molecule in contact with the surface is given in Figure 3.34. This orientation is consistent with Figures 3.31 and 3.33. 3.2.5 Combined Positional and Orientational Distributions. The positional and orientational probability distributions, as they depend on macro-sphere size, have been considered in the two previous sections. The range of the structure Chapter 3. Results and Discussion Figure 3.32: PNH{G,T) for lds diameter hard spheres in Model I N H 3 . Chapter 3. Results and Discussion Figure 3.33: PNH{8,T-) for 20da diameter macrospheres in Model I NH3. 153 Chapter 3. Results and Discussion 154 Figure 3.34: The dominant orientation of ammonia at a macrosphere surface. One of the three NH-bonds is pointing directly at the surface, at an angle of 180° to the surface normal vector. . Chapter 3. Results and Discussion 156 in the angular distributions was found to be from contact out to about r = 1.2ds. While the structure in <7oo°ms(r) w a s s e e n *-o be several solvent diameters in range, the contact peak greatly dominated the structure for all diameter spheres, ranging in magnitude from 2.1 for Id, spheres to 3.4 for 20dB spheres. The combined probability distribution 9oo°ms(r)Pfi{e,r), can be considered as the relative probability distribution for finding a solvent molecule located a given distance from the macrosphere surface with a given dipole orientation. Figure 3.35 shows this quantity for 20d„ macrospheres. The radial distribution emphasizes the angular structure in a region within a distance of 0Ada from contact which is the structure in the innermost layer of Figure 3.31, and it de-emphasizes the structure beyond this region. This is also true of the lds case, which is plotted in Figure 3.36 using the same scale. There is very little angular structure near the surface so even when it is amplified by multiplying by <7oo™«(r) the result is extremely unstruc-tured. A direct comparison of Figures 3.35 and 3.36 shows that a highly structured layer of approximately 0.4 solvent diameters forms at the surface of an uncharged macrosphere, and there is no similar structure around a solvent sized sphere. The average local den-sity within this layer is approximately 1.5 times the bulk density. The most probable direction of the dipole in this layer is at approximately a tetrahedral angle away from a surface normal vector. Since Figure 3.31 shows that the orientational structure carries on into the bulk and is significant up to a separation of about 2d„. The angular distribution beyond 0Ade is probably induced by the strongly defined relatively dense surface layer. Figures 3.37 and 3.38 show the corresponding distributions for the NH-bond vector. The dominant direction here is consistent with the dominant dipole direction, pointing the NH-bond.directly at the surface. Chapter 3. Results and Discussion 157 Figure 3.35: g%o°me{r)P»(0>r) f° r 20ds diameter macrospheres in Model I NH 3 . Chapter 3. Results and Discussion F i g u r e 3.36: 5™° m , ( r ' ) / v (^ T ' ) f ° r diameter h a r d spheres i n M o d e l I N H 3 . Chapter 3. Results and Discussion Figure 3.37: g™°ms(r)PNH(8,r) for 20de diameter macrospheres in Model I NH3. 162 Chapter 3. Results and Discussion Figure 3.38: go°-ms(r)PNH{0,r) for ld„ diameter hard spheres in Model I NH 3 . Chapter 3. Results and Discussion 165 3.3 Results for Ions at Infinite Dilution. The results given in this section are all for singly charged hard-sphere ions at infinite dilution in Model I ammonia. Some sample calculations for Models II and III have been performed, but the relative differences and the conclusions which can be drawn from them are similar to those which have already been discussed in the previous sections on model dependence. 3.3.1 Dependence of Solvation on Ion Size. The contact values, goo°is{die), for positive and negative ions as a function of ion diameter are shown in Figure 3.39. The positive ions of diameter Id, have a contact peak double that for negative ions, indicating a much stronger interaction between positive ions and solvent molecules for ions of this size. When the ions are increased in size to a diameter of 2d3 the value for both ions drops, but at a faster rate for positive ions, so the difference in contact values drops quickly to a few percent. At an ion diameter of slightly greater than 3da the difference in the rate that the contact values are dropping causes them to reverse order, and the value representing the negative ion remains slightly higher for all larger sized ions calculated. This asymmetry even for large ions is induced by solvent-solvent interactions and is discussed below. The contact values for ions of both charges go through a minimum at an ion diameter of about 5d8 after which there is a steady increase, and the positive and negative values converge. As the ion size increases, the electric ion field which the neighboring solvent molecules experience becomes smaller, to the point where the solvent-solvent interactions dominate the structure. Thus, the bigger the (monovalent) ion the more the solvent behaves as if it were just a macrosphere, and the structures which were seen to build at a macrosphere surface in the previous section are also seen to build here. This building of surface Chapter 3. Results and Discussion 166' Figure 3.39: goo°is{die) for hard-sphere ions in Model I NH3 as a function of ion diame-ter. The solid hne represents values for cations, and the dotted hne for anions, a) Ion diameters in the range ld„ < di < 2ds. b) Ion diameters in the range Id, < di < 20ds. Chapter 3. Results and Discussion 168 structure results in an increase in radial distribution function contact value for ions larger than about bd8. From these results it can be seen that solvent molecules are more closely held around small positive ions than they are around small negative ions, and as the ion size increases the differences diminish. The approximate range of the effect of ion field on solvent structure (about 5d,) can also be deduced. The interaction energy between an ion and the solvent Uj!?T, as given by equation 2.75, is another measure of solvation strength, and is shown in Figure 3.40. Once again the lower energy for positive ions of diameter lds indicates a greater interaction between the positive ions and the solvent, and hence a stronger solvation. For ions larger than 2d, in diameter, however, the trend is reversed and large negative ions are significantly more strongly solvated than large positive ions. The difference between the interaction energy for cations and anions appears to be significant over the entire range of ion size shown here, in contradiction to the results seen for radial distribution contact values. As will be examined below, the most energetically stable angular distribution of solvent molecules around oppositely charged ions are radically different, which is the underlying cause for the unequal interaction energies. This gives rise to different values for the multipole components of the energies as well. Figure 3.41 shows the charge-dipole contribution to the total ion-solvent energy. This is clearly the largest component of the total energy and for large ions it is the only significant component. The orientation of dipoles about a small positive ion is more energetically favorable than that for a negative ion. This can be understood in terms of the charge distribution around a solvent molecule. The charge density near the negative end of the dipole is strong and localized because it is reinforced by the linear quadrupole, but the positive charge density is less strong and more spread around the molecule because the positive end of the dipole has its charge partially cancelled by the quadrupole. Thus solvent molecules can align the negative end of their dipoles towards a positive ion while directing most of their positive charge Chapter 3. Results and Discussion 169 Figure 3.40: UjgT /NkT for hard-sphere ions in Model I N H 3 as a function of ion diam-eter. The lines represent the same ions as in Figure 3.39. a) Ion diameters in the range ld„ < d{ < Zda. b) Ion diameters in the range 2d„ < di < 20ds. Chapter 3. Results and Discussion 171 Figure 3.41: Um/NkT for hard-sphere ions in Model I N H 3 as a function of ion diameter. The lines represent the same ions as in Figure 3.39. a) Ion diameters in the range lds < di < 3de. b) Ion diameters in the range 2de < di < 20ds. Chapter 3. Results and Discussion 173 density away from the ion. The most energetically stable configuration around a negative ion, however, is a compromise between aligning the more spread out positively charged region towards the ion and keeping the strong negative end of the molecule away from the ion. The reason that the large difference is reduced as the ion size increases can be seen by examining the charge-quadrupole energies given in Figure 3.42. The contribution of UJQ is much less significant for ions above about 2.5ds in diameter. This effective range of the charge-quadrupole interaction reflects the range over which the quadrupole influences the orientation of ammonia molecules near an ion. To ions beyond this range the solvent molecules appear as dipoles, because the effect of the quadrupolar interaction, which decays as 1/r3 is no longer significant. As will be seen below, the total energies for big negative ions are larger (more negative) than those for big positive ions because the structure built by the solvent-solvent interaction, which is very similar around both big cations and big anions, is compatible with that found around a small negative ion. Thus the bulk solvent builds around big ions of both charges, the exact structure which is favorable to negative ions. The charge-octupole and charge-hexadecapole interactions are even shorter ranged, so the corresponding energy contributions to the total energy (Figure 3.43) are only significant for small ions. These energy contributions are both positive for lds diameter cations, indicating that the strong ion-dipole interaction comes from solvent orientations which are unfavorable for both the octupole and hexadecapole interactions. The orientation taken by the solvent around negative ions, however, is favorable for these higher multipoles, as indicated by negative energies. The coordination numbers, as defined by equation 2.78, are shown in Figure 3.44 and these too show larger numbers of solvent molecules in the coordination spheres around small positive ions than around small negative ions. As the ions size increases, the difference between the coordination numbers also decreases, and as the macrosphere structure builds, the coordination numbers are essentially unaffected by the single charge Chapter 3. Results and Discussion 174 Figure 3.42: Uio/NkT for hard-sphere ions in Model I N H 3 as a function of ion diameter. The lines represent the same ions as in Figure 3.39. a) Ion diameters in the range Id, < di < 5d,. b) Ion diameters in the range 5d, < di < 20d,. 175 Chapter 3. Results and Discussion 176 F i g u r e 3.43: a) Uio/NkT for hard-sphere ions i n M o d e l I N H 3 as a function of i o n d i a m -eter. T h e lines represent the same ions as i n F i g u r e 3.39. b) Um/NkT for hard-sphere ions i n M o d e l I N H 3 as a function of i o n diameter. T h e lines represent the same ions as i n F i g u r e 3.39. UIO/NkT Uw/NkT Chapter 3. Results and Discussion 178 Figure 3.44: Coordination numbers for Model I N H 3 molecules around hard-sphere ions as a function of ion diameter. The lines represent the same ions as in Figure 3.39. a) Ion diameters in the range Id, < di < 5de. b) Ion diameters in the range 5ds < di < 20ds. Chapter 3. Results and Discussion 180 on the ions. It is somewhat surprising that the coordination numbers for positive and negative ions converge as the ion diameter is decreased, and even cross for ions of diameter ldt. If, however, we observe the coordination sphere density, defined by equation 2.79 (Figure 3.45), it is clear that the solvent molecules are more closely packed around a small positive ion and the densities for anions and cations are in fact growing rapidly as the cation size decreases. The difference between the coordination numbers and the coordination sphere density is simply due to the volume of the coordination sphere, which is significantly less for cations than for anions of about ld„ diameter. This is itself an indication of how tightly solvent molecules are packed around ions. The way that solvent molecules are packed around oppositely charged ions is related to single ion partial molar volumes. Two possible definitions of single ion partial molar volumes for ions at infinite dilution were given in section 2.2.4. Each of these definitions for positive and negative ions, as well as the unambiguously defined solute partial molar volume, t^0, are shown in Figure 3.46. It is clear that only one of these definitions, that of V°, behaves in accord with the other quantities observed here. Using the definition V?" of single ion partial molar volumes one would predict from the results for this model of ammonia that the local density of solvent around positive ions would be lower than that for negative ions at all ion sizes, but that is the opposite of what is actually found. This definition is not in agreement with the total ion-solvent interaction energies for small ions either. The cation-solvent energy is more negative than the anion-solvent energy which should predict a stronger interaction and hence a larger change in volume when added to solution. Trying to deduce single ion behavior using the definition t^°* to divide solute partial molar volumes, V°, could lead to erroneous conclusions, whereas using the definition Vi 0 would not lead to such errors. Both definitions of the single ion values give the correct, unambiguous definition of the solute partial molar volume. As the ion size increases the relative difference between the single ion quantities for both Chapter 3. Results and Discussion 181 Figure 3.45: Coordination sphere density for Model I N H 3 molecules around hard-sphere ions as a function of ion diameter. The lines represent the same ions as in Figure 3.39. a) Ion diameters in the range lds < di < 2dt. b) Ion diameters in the range 2d, < di < 20d,. 182 Chapter 3. Results and Discussion 183 Figure 3.46: Partial molar volumes at infinite dilution for hard-sphere ions in Model I N H 3 as a function of ion diameter, a) Ion diameters in the range ld„ < di < 2d„. The solid, dashed, dotted and dot-dashed lines represent V°, V2, V+* and repectively. The long-dashed Hne represents V2 . b) Ion diameters in the range lds < di < 20dB. The long-dashed hne represents V2°, and on the scale used in this figure, the four single ion quantities in a) are superimposed forming the solid line. Chapter 3. Results and Discussion 185 definitions and both ion charges diminishes, and approaches V2°/2. This is consistent with the diminishing ion field. 3.3.2 Solvent Distribution About Ions. The ion-solvent radial distribution functions for positive and negative ions (Figure 3.47) are remarkably similar. Each shows greatly increased structure over the uncharged hard spheres, and a tendency for all maxima and minima to occur at lower values of r. The electric field due to each ion draws the bulk solvent in around itself, and the large contact peaks result. This is perhaps the most natural illustration on a molecular level of the macroscopic electrostriction effect, and the relative size of the contact structure is consistent with the conclusions above. The contact peaks are the only significant differences between the results for oppositely charged ions, and they give rise directly to the differences in coordination sphere density. The manifestation of the asymmetric solvation is found in the angular contributions to the total correlation function. Figure 3.48 shows the function h°l\3(r) for the same result as in Figure 3.47, with the positive ion function multiplied by —1 to facilitate comparison. This function is directly related by equation 2.76 to the charge-dipole energy, and the large differences in the size of the contact peak, as well as more overall structure, leads to the greater ion-dipole energy for cations than for anions. By comparison the result for a hard-sphere solute is almost featureless. Even greater differences in the function h°l2ie(r), related to charge-quadrupole energy, can be seen in Figure 3.49, with the anion-quadrupole function closer to the hard-sphere solute case than to the cation result. The function /&o(Hs(r) i s given in Figure 3.50. The positive ion projection has been multiplied by -1 in this figure also, so the difference in the sign of the charge-octupole energy contribution from oppositely charged ions is expected. The most dominant overall feature of Figures 3.47 to 3.50 is the very large contact peak. This dwarfs all the other structure, and dominates the energy contributions Chapter 3. Results and Discussion 186 Figure 3.47: g%™e(r) for ions in Model I NH 3 . The sohd, dotted and dashed lines are results for positive ions, negative ions and hard spheres, respectively. Chapter 3. Results and Discussion 188 Figure 3.48: h^^r) f° r i ° n s in Model I N H 3 . The solid, dotted and dashed lines are results for positive ions, negative ions and hard spheres, respectively. The result for positive ions (solid Hne) has been multiplied by -1 to facihtate comparison. 189 Chapter 3. Results and Discussion 190 Figure 3.49: ^oofistj) ^ o r i ° n s m Model I N H 3 . The sohd, dotted and dashed lines are results for positive ions, negative ions and hard spheres, respectively. The result for positive ions (sohd hne) has been multiplied by -1 to facilitate comparison. 0.0 0.8 1.6 2 . 4 3 .2 4.0 4 .8 (r-du)/d, Chapter 3. Results and Discussion 192 Figure 3.50: h^ia(r) for ions in Model I N H 3 . The sohd, dotted and dashed lines are results for positive ions, negative ions and hard spheres, respectively. The result for positive ions (solid line) has been multiplied by -1 to facilitate comparison. Chapter 3. Results and Discussion 194 dependent on these projections. The angular distribution about a Id, diameter negative ion is given by the functions Pn(9,r) and PNH(9,r) is shown in Figures 3.51 and 3.52, respectively. These plots are remarkably similar to the graphs in the previous section for the structure of ammonia near the surface of a macrosphere (Figures 3.31 and 3.33), except the peaks are higher and more defined. The negative ions seem to enhance the tendency for ammonia to form this ordered structure. It is consistent with the charge distribution of an ammonia molecule to form this structure near a negative ion, because one of the three equivalent points of highest positive charge, an NH-bond, is pointing directly at the surface of the negative ion. This necessitates that the dipole not point directly at the ion which is the most favorable alignment for the dipole alone. This is the nature of the energetic compromise discussed above, and gives rise to the smaller charge-dipole interaction energy term for small negative ions, compared to the equivalent quantity for positive ions. The corresponding plots for the angular distribution about a positive ion are given in Figures 3.53 and 3.54. The structure here is distinctly different from that for ammonia near a negative ion or a macrosphere. The only feature in the dipole distribution (Figure 3.53) is the massive peak at contact (note the scale on the vertical axis), at an angle of 0° from the surface normal vector, indicating that the solvent dipoles are lined up directly away from the positive ion, with the negative end of the dipole pointed directly at the ion. This is extremely favorable energetically and gives rise to the large charge-dipole energy which in turn dominates the total energy of interaction. The NH-bond angle at contact is consistent with the dipole structure, with the only peak at about 75° to the surface normal, or at about a tetrahedral NH-bond/lone-pair angle. The structural effects observed near ions here bring new light to the previous section on ion size. The ordered structure of ammonia, driven by the strong solvent-solvent in-teractions, builds up around macrospheres and macroions of both charges alike, because Chapter 3. Results and Discussion 195 Figure 3.51: p^{9,r) for Model I N H 3 molecules near an anion of diameter ldg. 19G Chapter 3. Results and Discussion 197 Figure 3.52: PNH{8,T) for Model I NH3 molecules near an anion of diameter ld„. 198 Chapter 3. Results and Discussion Figure 3.53:.p (^ff,r) for Model I N H 3 molecules near a cation. 200 Chapter 3. Results and Discussion Figure 3.54: pNii(Q,r) for Model I N H 3 molecules near a cation. 202 Chapter 3. Results and Discussion 203 for large particles the ion field is too small to compete with the solvent's internal struc-turing. This structure produced by the solvent-solvent interactions is more favorable energetically for anions. As a negatively charged ion gets smaller the curvature of the ion surface disrupts the solvent structure, and decreases in contact value and coordina-tion sphere density are observed. As this happens, however, the ion field becomes more significant, and it serves to rebuild the same structure, and so further decreases in ion size reverse the trends in contact value and coordination sphere density after they go through a minimum for medium sized negative ions. The interaction energy is increased as the ion field increases. The story for positively charged ions is similar at first, with the disruption of the macrosphere structure as the cation decreases in diameter. But, as the cation field increases it induces strong orientational changes in the solvent around it, and the dipoles which were pointing down toward the surface at an angle of about 110° from a surface normal vector, swing around and point directly away from the surface of the ion. This is is extremely favorable energetically for cations, as the increase in ion-solvent energy which surpasses that for negative ions attests. This favorable cation structure is consistent with the big increases in contact value and coordination sphere density which occur as the ions are reduced in size. Chapter 3. Results and Discussion 204 3.4 Dependence of Ionic Behavior on Concentration. This section will focus on the marked deviations between experimentally determined mean molar activity coefficients at low concentrations and the Debye-Huckel hmiting law. These deviations have been attributed in the literature to ion pairing, or ion associ-ation, so some discussion of these ideas will be necessary. The quantities and equations that will be used are primarily those from the Kirkwood-Buff theory, which are given in section 2.2.1. The reliability of the HNC approximation when applied to the primi-tive model (Model V) has been questioned in the literature [18,57] for certain ranges of temperature and density. Model V can be defined in terms of the reduced temperature, T" = ekTd±/\qiqj\, and the density, p2 (which will be given throughout this section in dj3 units), and using these quantities the ranges within which HNC results should be treated carefully are approximately 0.015 > p2 > 0.00015 and T* < 0.14. We note that the numerical values of these reduced parameters are slightly different from those usually given because we have chosen the diameter of the ions to be di = 1.52dt rather than the usual di = 1 so Models III, IV and V can be readily compared. As has been well documented [17-22, 24, 50], the primitive model seems to have a phase change at relatively low density, but the exact position of the critical point and the ionic gas/ionic liquid coexistence curve is not yet known. Using the stability criteria introduced in section 2.2.2, the HNC approximation can also be used to predict this phase behavior, and as will be described below a spinodal decomposition point for a wide range of p2 values can be found. The spinodal decomposition line found in this manner is given by the dotted line in Figure 3.55. It can be seen from Figure 3.55 that almost the entire range of T* and p2 where the HNC results are considered unreliable is either subcritical or close to the spinodal decomposition hne. For purposes of comparison the critical points from other literature sources [17,19,20] are also given in Figure 3.55. As has been Chapter 3. Results and Discussion 205 Figure 3.55: The Primitive Model Spinodal Decomposition Line (dotted line) as Predicted by the HNC Theory. The numbered points all represent critical points from the following sources. 1) HNC theory; this work, 2) The Mean Spherical Approximation [19], 3) Monte Carlo data [24,25]; this value should be regarded as a lower bound for both T* and p2, 4) Reference [20], 5)-8) The theories denoted L, PUB 2 , riUB2, and r 2 U B 2 from reference [17], respectively. The solid hne segment represents the temperature and the range of concentrations for the experimental data for AgN0 3 taken from reference [8]. 0. u 0 . 1 2 H T' 0 . 1 0 H 0.08 H • 5 6* 7#* •4 •2 0.06 | i i i 11 m i i i i 11 m i i i i 11 • i• | — 0.000001 0.00001 0.0001 0.001 ?2 I I I I I I I I I I 11 I I I 0.01 0.1 to o G i Chapter 3. Results and Discussion 207 pointed out by Gillan [19,32] the MSA result is not reliable and the purpose of the work by Stell, Wu and Larsen [17] was to demonstrate the existence of a critical point in theories which differed substantially from each other, rather than give definitive numerical values. These alternative theories, therefore, can not be considered more accurate than the HNC theory. The internal energies given by the HNC theory are in good agreement with recent Monte Carlo results [24] in all regions except near the critical point predicted by the HNC theory, where they differ from the Monte Carlo values by as much as 10%. Figure 3.55 also contains an estimate of the critical point from these Monte Carlo calculations, although this point should be regarded as a lower bound for both T" and p2 [24,25]. It is probable that the energies calculated by the HNC theory deviate from the Monte Carlo results in this region because the Monte Carlo simulations predict the critical point at higher density. It is also probable, therefore, that the HNC theory predicts the critical point at a temperature T" which is too high and a density p2 which is too low. It has been pointed out [57] that the addition of a bridge diagram term to the HNC theory changes the distribution functions significantly, so it is clear that care is needed in interpreting the results of the HNC theory for Model V. 3.4.1 Examination of Experimental Results. There is not a large amount of experimental data on electrolytes in liquid ammonia, compared to the vast wealth of experimental data for electrolytes in water. Mean molal activity coefficients, f±, have been determined in a series of papers by Baldwin, Gill and others [5-9], and these have been converted to the concentration units used in this study and plotted in Figures 3.56 and 3.57. These experiments use conductance measurements to determine the ratio of the activity coefficient for the electrolyte in a cell where the concentration of electrolyte is varied, to the activity coefficient for the electrolyte in a reference cell in which the concentration is held fixed. The determination of the activity Chapter 3. Results and Discussion 208 Figure 3.56: Experimental mean molar activity coefficients for 1:1 Electrolytes. The solid hne represents the Debye-Hiickel limiting law for 1:1 electrolytes in a solvent with the experimental dielectric constant of N H 3 (22.1). The triple-dot dashed Hne repre-sents experimental results for AgN0 3 [6], the dot-dashed Hne and dashed Hnes represent N H 4 N 0 3 and NH 4I results from reference [6], respectively. The dotted and long dashed lines represent the re-analysis of the N H 4 N 0 3 and NH 4I results given in reference [7], respectively. AU results are for — 40°C. to o CO Chapter 3. Results and Discussion 210 Figure 3.57: Experimental mean molar activity coefficients for 1:2 Electrolytes. The sohd hne represents the Debye-Hiickel hmiting law for 1:2 electrolytes in a solvent with the experimental dielectric constant of N H 3 (22.1). The dashed and triple-dot dashed lines represent the Cu(N0 3 ) 2 and Pb(N03)2 results given in reference [8], respectively. All results are for —40°C. Chapter 3. Results and Discussion 212 coefficients from the ratio is done in a number of ways. If the experimental ratios were for a system which followed the Debye-Hiickel hmiting law at low concentration, a plot of the natural logarithm of the ratio versus the square root of the concentration would be linear and the extrapolation to zero concentration would unambiguously determine the activity coefficients. The experimental activity coefficient ratios do not follow such a linear relationship over any concentration range, but plots versus the 4-th root of the concentration are remarkably linear. This 4-th root relationship is the analysis which was used to determine the individual results in Figure 3.57 and some results in Figure 3.56 [6-8]. The use of n-th root graphical methods may at first seem arbitrary, but has been used enough that there have been attempts to justify the use of cube root relationships [118] instead of the usual Debye-Hiickel square root behavior. Using different ra-th root relationships makes a significant difference to the final activity coefficients obtained, corresponding to different vertical displacements of the entire curve for each n. The activity coefficients obtained for the electrolyte-ammonia systems clearly do not exhibit Debye-Hiickel behavior, particularly for AgN03 and the 1:2 electrolytes. Had the usual square root law been used the results would have deviated even further from the limiting law. It should be noted that none of the salts used were near their solubility limits over the experimental concentration ranges. The behavior found in the experiments was attributed to ion pairing [7], a contention which is not confirmed by Raman spectroscopy [9], and a different analysis of the same raw data for N H 4 N O 3 and N H 4 I , in terms of the relationship between the degree of association, 1 — a, ion diameter and activity coefficient (equation 2.158b) was performed in reference [7]. Relationships between ion association constants, K^, the degree of association and activity coefficients were found using the analysis of Bjerrum (see section 2.3.1). The association constants obtained in this manner were extremely dependent on concentration. When ion triples were included in the analysis some of the concentration Chapter 3. Results and Discussion 213 dependence was removed [7]. The resulting curves appeared as if the original curves had simply been shifted up towards the limiting law (see Figure 3.56). The non-linear shape of the curves did not appear to change. It should be noted that these results depended on the minimization of the standard deviation from the functional form given in equation 2.158b, achieved by varying the association constant. The value used for the distance of closest approach of positive and negative ions was taken to be the distance inside which ions are considered paired by Bjerrum's theory (Rin equation 2.156). This was evaluated to 16.22A which corresponds to 5-07dt. The effective crystal radii [119] of N H 4 , N O 3 and I - are 0.47^, Q.52dt and 0.64d!a respectively, so the distance of closest approach used corresponded to ions separated by about 4 solvent molecules. This is not consistent with the ion pairing theories of either Bjerrum [10,11] or Fuoss [11-13], nor with the basic concept of ion pairing. These experimental results simply display the problems inherent in fitting results which display behavior distinctly different from that of the Debye-Hiickel limiting law to rather arbitrary functional forms. When analyses are performed in such a way that parameters representing physical quantities are determined, these quantities are often decidedly unphysical in size or even sign. It seems clear, however, that while it is difficult to exactly determine experimental values for ln y±, regardless of the method of analysis used, the experimental results lie below the limiting law and follow a curve which is concave upwards. The idea that all ion pairing theories share is that some ions form a new species and no longer behave as individual ions in solution, therefore the experimental quantities should be looked upon,as functions of the concentration of "free" ions, ap, rather than of the stoichiometric concentration p. It is clear from equation 2.158a, however, that if fewer free ions are present, lny± will deviate above the limiting law, not below it. The inclusion of ion diameter dependence also causes a deviation above the limiting law. While deviations above the limiting law are often [3] (and perhaps correctly) attributed Chapter 3. Results and Discussion 214 to ion pairing, deviations below the hmiting law seen in these experimental activity coefficients for electrolytes in ammonia are also attributed to this phenomenon. There seems to be no simple theory which supports this contention. It can also be seen from the functional form of equation 2.158b, that deviations due to ion pairing and ion size should result in a plot of \ny± versus square root concentration which is concave upwards. This is seen in the experimental data given, and could be due at least in part to ion pairing, but the position of the curve, particularly that of AgN0 3 , must also be affected by other phenomena which result in much more drastic changes in In y± than those caused by ion pairing. The experimental activity coefficients could possibly be explained by some kind of chemical interaction between ions and the solvent. The formation of ion-solvent com-plexes is probable with some metal cations, but any simple analysis of this situation [3] simply gives the cation a larger effective diameter, for which simple theories (see equation 2.158b) predict a deviation in the wrong direction. In summary, experimental activity coefficients deviate from the Debye-Hiickel hmiting law in two ways. The graphs of lny± versus ^/p~l are concave upwards which is probably due to the finite size of ions, but may also be due in part to ion association. These graphs also lie below the hmiting law, which can not be explained by simple ideas of ion pairing. Clearly some explanation for activity coefficients which deviate below the limiting law is needed, particularly for such large deviations as are seen for AgN0 3 in Figure 3.56 and the 1:2 electrolytes in Figure 3.57. 3.4.2 Results for Electrolytes in Liquid Ammonia. Calculations for Model III (molecular solvent), IV (McMillan-Mayer) and V (primitive model) were carried out for equal size ions of diameter 1.52d„ at — 35° C, over a range of concentrations between infinite dilution and 0.25M. For Model III, calculations were Chapter 3. Results and Discussion 215 carried out up to a concentration of just under 3.0M in an attempt to find a solubility limit, but the charge neutrality condition calculated above this high concentration was such that the results were considered to be unreliable. The calculated values of the Kirkwood G + _ , as given by equation 2.119, are plotted versus ion number density in Figure 3.58. The dependence on ion diameter is clearly an important effect as the Model III result can be seen to closely follow the d-dependent function given by equation 2.140b. As a test of the importance of ion pairing, two definitions of ion pairs were used and the results rescaled accordingly. In the first, the coordination number (C n) route as given in equation 2.164 was used to find the proportion of free ions. The same G + _ values were then plotted versus ap2. In the second, the long range behavior of w +_(r) was assumed to have the same functional form as the low concentration limit (equation 2.136a) and the difference between the calculated neff and the debye K was used to find the proportion of free ions (see equation 2.166), and the same rescaling technique was used. For the Model III calculations at concentrations above p2 — 2.5 X 10~"3 the exponential form was sufficiently damped that short range contributions to w+_(r) made it impossible to numerically determine K e / / . The proportion of free ions a found by each of these methods is given in Table 3.7. The two routes agree more closely as concentration is decreased, but the C n route predicts more free ions as the concentration is decreased whereas the the neff route predicts fewer free ions. The value of a given by the neff route is probably a more realistic measure because it finds the number of free ions directly. The increase with concentration probably arises from a reduction ih the range of ion-ion correlations due to Debye screening. The C n route depends on the assumption that all negative ions found around a central positive ion inside the first local minimum in g+-(r) are paired (rather than in ion triples), which becomes less valid as the concentration is increased. As can be seen in Figure 3.58, for both definitions of a the estimates of ion pairing effects are small compared to the effect of ion diameter, and also Chapter 3. Results and Discussion 216 Figure 3.58: G + _ vs. p2 for 1.52<ia diameter ions in NH 3 . The sohd hne represents the Debye-Hiickel hmiting law, the dashed hne is the related d-dependent quantity, the dotted hne is the result from Model III. The long-dashed hne represents the Model III result where the result has been rescaled to the number of free ions calculated using coordination numbers. The dot-dashed hne is similarly rescaled with the number of free ions calculated from the long-range behavior of w +_(r). 0 1 2 3 4 5 6 7 Chapter 3. Results and Discussion 218 P2 a (C n route) a (neff route) 0.0004853 .6637 .7179 0.0004972 .6642 .7227 0.0005550 .6562 .7345 0.0005950 .6492 .7364 0.0006960 .6325 .7437 0.0007930 .6165 .7493 0.0009940 .5868 .7633 0.0011890 .5608 .7769 Table 3.7: Proportion of Free Ions a Calculated From Model III. in both cases these adjustments move the Model III results away from the Debye-Hiickel limiting law rather than towards it. In summary, the dependence on particle diameter is clearly more important than ion pairing concepts in explaining the deviation of G+_ from the Debye-Hiickel limiting law. The Model III result in Figure 3.58 does not go below p2 = 0.0004853, because this is the lowest concentration where numerical solution could be obtained with these ion diameters. Much time and effort was expended in trying to get solutions below this con-centration by varying concentration, ion diameter, ionic charge and dielectric constant (indirectly by varying the magnitude of the solvent dipole). Models IV and V were also calculated over this same concentration range, and each of these models had a minimum solvable concentration limit, p2 = 0.0003659 and p2 — 0:0003878 respectively. At ex-tremely low concentrations it was again possible to get Model IV and V solutions using the methods described in section 2.4.3, but in this case there was a similar concentration maximum for each model beyond which no solution could be achieved. This behavior in the primitive model (Model V) has been attributed to a phase change, as mentioned above. It is somewhat surprising for similar behavior to be found at such a similar Chapter 3. Results and Discussion 219 concentration in a full molecular model. The interpretation of this phenomenon will be examined more closely below. The numerical procedures used to get these low concentra-tion solutions reached their reliable lower limit at a concentration of p2 = 1 x 10 ~ 6, where, as can be seen in Figure 3.59, there are still significant deviations from the Debye-Hiickel hmiting law, but the calculations now deviate in the opposite direction to that for higher concentrations. The reason for these deviations from the hmiting law is not clear, but is perhaps due to the proximity to a spinodal decomposition point as predicted by the HNC theory. Equation 2.129 gives the derivative of the logarithm of the activity coefficient with respect to ion density. The equivalent expressions in the low concentration limit are given by equations 2.146b and c from which it can be seen that ^fp^ multiplied by the Debye-Hiickel limit gives a constant. This product is plotted in Figure 3.60 and again it can be seen that the calculations from Models IV and V are closer to the DH hmiting law than the related d-dependent result, but they can not really be considered as being close to hmiting behaviour. In fact, in the range where results were feasible, they appear to be approaching the wrong zero concentration limit. As a measure of the approximations made in equation 2.146b for the DH hmiting law, the more fundamental DH hmiting law expression for G+- (equation 2.140a) was inserted into equation 2.129 and is also given in Figure 3.60. The concentration range shown in Figure 3.60 is still not low enough for the limit taken in equation 2.146b to be seen. The same quantity was calculated for Models III, IV and V in the higher concentration range, and is given in Figure 3.61. Here none of the models approach the DH hmiting law, but they do follow the related d-dependent quantity closely, and the results from Models III and V are almost identical. Again the similarity between the primitive model and the full molecular model is surprising, but the proximity of both results to the d-dependent equation indicates that errors in the HNC theory are probably not large here. Because calculations were not possible over the entire range from zero concentration Chapter 3. Results and Discussion 220 Figure 3.59: vp2G+- vs. p2 for 1.52d, diameter ions in NH 3 . The solid line represents the Debye-Hiickel limiting law, the dashed hne is the related d-dependent quantity, the long-dashed Hne is the result from Model IV and the dot-dashed Hne is the result from Model V. Chapter 3. Results and Discussion 222 Figure 3.60: y/p~2(dlny±/dp2)T,P v s - \/p~2 for l-52d„ diameter ions in N H 3 . The sohd Hne represents the Debye-Hiickel Hmiting law, the dashed Hne is the related d-dependent quantity, the long-dashed line is the result from Model IV and the dot-dashed Hne is the result from Model V. The dotted line is obtained by inserting the Debye-Hiickel Hmiting law for G+_ from equation 2.140a into equation 2.129. The similarity between the sohd and the dotted Hne is simply the coincidence that \qiqj\/ekTd± is close to 8 for these calculations. Chapter 3. Results and Discussion 224 Figure 3.61: ^/P2(dlny±/dp2)T,p v s - yfPi f° r l-52<iB diameter ions in NH 3 . The sohd hne represents the Debye-Hiickel limiting law, the dashed hne is the related d-dependent quantity, the dotted line is the Model III result, the long-dashed hne is the result from Model IV and the dot-dashed hne (which almost coincides with the Model III result) is the result from Model V. p 2 \ eP2 ) T ) P 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 • 0.08 t o t o Chapter 3. Results and Discussion 226 up, it is not possible to calculate \ny± for comparison with experimental results because the Kirkwood-Buff theory only gives the derivative of this quantity. In order to make such a comparison some assumption about the low concentration behavior must be made. Using the assumption that on a plot of lny± versus yfpi, the missing range of concentra-tion can be replaced by a straight Hne of the same slope as that given by the derivative in equation 2.129 for the lowest available concentration, the values plotted in Figure 3.62 can be calculated. It should be noted that any assumption made for the values of lny± in the missing concentration range corresponds to choosing a starting point for the curves. The only thing that can be asserted from this graph is that the curves are aU concave upward as are the experimental results. This particular choice of starting point is the highest (least negative) starting point possible without imphcitly assuming that either the curve is discontinuous, or it must have an inflection point and be concave downward over some range of concentration. It is as close to the .Debye-Hiickel limiting law as one can place the curves without incurring these side-effects. If these calculated curves were placed on Figure 3.56 with the experimental results (assuming the 5° difference in temperature is negHgible) they would faU between the two pairs of curves representing the two analyses of the data for N H 4 N O 3 and N H 4 I . 3.4.3 Indication of a Molecular Model Phase Change. The previous section leaves two important questions unanswered. These questions are: 1. Is the range of concentrations where calculations were not possible for the molecular solvent model a numerical artifact, or is it related to the phase change that is known to exist in the primitive model? Chapter 3. Results and Discussion 227 Figure 3.62: lny± vs. ^/pl for l.b2d„ diameter ions in NH 3 . The starting points of the curves are explained in the text. The sohd hne represents the Debye-Hiickel hmiting law, the dashed hne is the related d-dependent quantity, the dotted hne is the Model III result, the long-dashed line is the result from Model IV and the dot-dashed hne is the result from Model V. - 2 . 5 H 1 1 1 1 1 1 1 1 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 to to oo Chapter 3. Results and Discussion 229 2. What is the explanation for the large deviations in experimental activity coefficients below the DH hmiting law, particularly those of AgNC-3 and the 1:2 electrolytes? There is considerable evidence that the existence of a concentration range over which no solutions were possible for Model III is not a numerical artifact. Instabilities in calculations which have a numerical source usually have a gradual onset. As an input parameter (in this case density of one of the components) is adjusted, various quantities become less and less precise, and results appear less and less rehable. This was not the case with any of the calculations considered here. As discussed above, many input parameters were varied near the end of the solvable range in concentration, but no matter which input parameter was being considered, the onset of the divergence of the numerical solution was quite sudden. The last obtainable point in all cases showed no indication of numerical instability. The similarity in the concentrations at which Models III and IV no longer have a numerical solution indicates some connection to the phase change in the primitive model. If such a phenomenon were occuring in any of the models, certain thermodynamic quantities should be seen to behave in predictable ways. Ursenbach and Patey [100] have derived four stability criteria in terms of the quantities given by Kirkwood-Buff theory, and these are given in equations 2.131a to 2.134c. For Models IV and V the compositional stability criteria are constant, equal and positive, so only mechanical stabihty need be considered for these models. For Models IV and V the mechanical stabihty criteria both reduce to the quantity l / / 9 2 G + _ , which is plotted in Figure 3.63, over a very narrow concentration range near the limit where the theory could be solved. At a spinodal point this quantity should approach zero. This quantity has also been calculated for the full molecular solvent model (Model III) and is given in Figure 3.63 but it should be noted that this quantity is not a stabihty criterion for Model III. The results for the full Chapter 3. Results and Discussion 230 Figure 3.63: • 1/>2G+- vs. concentration for Models III, IV and V. The dotted hne represents Model III, the long-dashed hne represents Model IV, and the dot-dashed hne is the Mode] V result. P2<2 + _ 1 . 2 8 1 . 2 6 -1 . 2 4 -1 . 2 2 -1 . 2 -1 . 1 8 - • i 1 . 1 6 -] 1- 1 1 1 1 r 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 : ( M ) to CO Chapter 3. Results and Discussion 232 molecular model (Model III) and and the primitive model (Model V) are both increasing (as concentration is decreased towards the instability point) over a long range of higher concentrations. As all models approach the concentration at which they become unstable, the quantity l/p2G+- decreases rapidly towards zero, not only for the models for which this indicates an instability, but also for Model III. The same quantities for the other side of the concentration break are plotted in Figure 3.64 for the two models (IV and V) for which calculations were feasible. Again this stabihty criterion behaves as one would expect near a spinodal point, but the curvature at each concentration limit is not nearly as distinct here because of the range of concentrations needed to display the results from both models, and because the onset of the instabihty here is extremely sudden. However, the derivative quantity for these models which is plotted in Figure 3.60 by the nature cf equation 2.129 amplifies the behavior of the instabihty criterion l/p2G^ , and the downward turn can be seen more clearly. When a full molecular model is considered, the ions are now only one component in the system. The instabihty criterion examined above is essentially the reciprocal of the isothermal compressibility which is a measure of mechanical instabihty. If a similar phenomenon of ionic condensation within the molecular solvent medium is occuring, the corresponding instabihty will be compositional rather than mechanical. Both mechanical and compositional stabihty criteria from equations 2.131a to 2.132b and are plotted in Figure 3.65 and it can be seen that the mechanical instability that was present in Models IV and V is now replaced by a compositional instability in Model III, due to essentially the same phenomenon. There exists in Model III a solute concentration at which the system becomes compositionally unstable. In answer to the first question above, therefore, we find there is no evidence of numer-ical problems in any of the calculations below an ionic concentration of 3.0M. Since the stabihty criteria for all three models behave, at the limits of the available concentration Chapter 3. Results and Discussion 233 Figure 3.64: l/p.2G+- vs. concentration for Models IV and V. The long-dashed hne represents Model IV and the dot-dashed Hne represents Model V. Chapter 3. Results and Discussion 235 Figure 3.65: a) Mechanical Stability Criteria for Model III. The dashed hne represents the quantity SG/{DG -f VG) and the solid hne represents the quantity SG/DG- b) Com-positional Stability Criteria for Model III. The dot-dashed Hne represents the quantity 1/SG and the dotted line represents the quantity (DG + VG)/DQSG-Mechanical Stability 14.5-13.5-.020 T 1 1 1 1 1 r 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 :{M) 0.66 0.64-0 .62-0.60 0.58H 0.56 0.54 Compositional Stability * - \ — i — i — i 1 1 — i — i — 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 c(M) Chapter 3. Results and Discussion 237 range, in the way that would be expected for systems approaching a spinodal point, we conclude that a combination of the models and the theory used indicate the existence of an instability in the full molecular solvent model which we interpret as a separation into two ionic solutions of differing concentrations. While the indication of a phase separation in the theories used here does not in turn indicate that such a phenomenon occurs in real electrolyte-ammonia solutions, it is possible to examine the behavior of the activity coefficient, y±, assuming such phase behavior does exist. If we now re-calculate the quantities given in Figure 3.62 we can use the assumption that as the concentration of electrolyte p2 is increased the solution will reach a point where it begins to separate into two hquid phases. There will be a range of concentration where the two hquid phases exist in contact, and as more electrolyte is added to the solution, only the relative volumes of the "dilute" (d) and "concentrated" (c) phases will change, while the concentrations within each phase, p2d^ and p2°\ will remain constant. When enough electrolyte has been added for the whole system to have a stoichiometric concentration equal to p2c\ the dilute phase will disappear leaving only the concentrated phase, and further addition of electrolyte will result in the whole solution increasing in concentration while remaining homogeneous. Since the activity of each ionic species must be equal within the two phases in equilibrium, and since the activity in each phase will remain constant while both phases exist in contact, a relationship between the mean molar activity coefficient of the salt in each phase can easily be found. We simply write (<0 _ ( 3 . 1 7 9 ) Chapter 3. Results and Discussion 238 where the subscript i denotes either ionic species and we have used p + = p_ = p2. From this equation we have The ratio in the final term must always be greater than one by definition and hence the logarithm must always be positive. The quantity lny±^ is well defined by the Kirkwood-Buff theory and can be calculated by integrating the derivative quantity in equation 2.129 from p2 = 0 to p2 — p2^- Since for our calculations p2d^ is very small we assume y±^ % 1 so we set lny±^ to zero. If we make the final assumption that the concentration of phase separation is close to the spinodal decomposition point for both dilute and concentrated phases, we can now completely determine lny± for p2 > p2c^ for Models IV and V. It was not feasible to calculate the low concentration spinodal point for Model III, but since the results for Models III and V were very similar in the concentration range p2 > p2°^ we assume for Model III is similar to that of Model V in order to get a qualitative idea of the behavior of lny ± . Equation 3.180 thus gives \ny±^ and equation 2.129 gives \ny± for p2 > p{2]. The result is given in Figure 3.66, and an enormous drop in lny± can be seen. If such a phase change were to exist, it is not clear what experiments would measure in the stoichiometric concentration range between the stability limits. Nor is it clear that the phases would completely separate and form an observable interface. It is probable, however, that if such phase phenomena existed the measured activity coefficients would be considerably below the Debye-Hiickel hmiting law as is the case for AgNC^. It is also probable that attempting to interpret the data by fitting it to functional forms such as equation 2.147b will result in decidedly unphysical parameters unless the phase change is accounted for explicitly. This has also been observed for the experimental data examined here. Although this suggestion is obviously highly speculative, there is some (3.180) Chapter 3. Results and Discussion 239 Figure 3.66: lny± vs. yfp], for 1.52d„ assuming a phase change. The sohd hne represents the Debye-Hiickel Hmiting law, the dashed Hne is the related d-dependent quantity, the dotted Hne is the Model III result, the long-dashed Hne is the result from Model IV and the dot-dashed Hne is the result from Model V. to O Chapter 3- Results and Discussion 241 experimental evidence that it might be valid. Such phase separations in low concentration solutions of alkah metals in hquid ammonia have been observed since the beginning of this century [26,27]. The further addition of electrolytes has been seen to greatly affect the temperature at which such systems separate into two hquid phases [30-33]. The sohd hne segment in Figure 3.55 indicates the temperature and range of densities at which the experimental AgNOs activity coefficients [8] were measured. For the HNC theory for Models III, IV and V this hne segment falls in a subcritical region of the phase diagram. It is not known where this line segment would he with respect to the coexistence curve for the exact solution to the primitive model, nor how closely Model III places the spinodal hne to that of real ammonia solutions, nor even if such phase separations occur in this range of temperatures and concentrations. It does seem clear, however, that if activity coefficients were measured below the critical temperature the results in the concentrated phase should deviate below the Debye-Hiickel hmiting law. In summary, we find phase instabilities are predicted by the theories used here when solved using three widely differing models which in turn are apphed to low concentration electrolyte solutions in hquid ammonia. If we then assume that such a phase separation occurs in real electrolyte-ammonia solutions, we find that one consequences of this as-sumption is that activity coefficients must deviate from the Debye-Hiickel hmiting law in the same direction as the experimental results. Further, the usual ideas of ion pairing do not adequately explain the experimental data because they predict deviations from the limiting law in the opposite direction to that which is actually found. We suggest, therefore, that at low concentrations some electrolyte-ammonia solutions may undergo phase separation and that this may be an explanation for the unusual concentration dependence of lny± observed for some systems. Chapter 4 Summary and Conclusions This thesis has been concerned with the structural, thermodynamic and dielectric properties of pure ammonia, ammonia near charged and uncharged curved surfaces, and the very large deviations from the Debye-Hiickel hmiting law in the behavior of some electrolyte-ammonia solutions at low concentrations. Five polar/polarizable models for ammonia, ranging from a full molecular model with C3,, symmetry and electric multipoles up to hexadecapole order (Model I), gradually simplified through to the primitive model of ions in a dielectric continuum (Model V) were considered in order to discover the effects of features of individual models. Ions and macrospheres when used were modelled as hard spheres with or without charge. The RHNC/SCMF integral equation theory was solved for Models I, II and III, and the HNC theory for Models IV and V. The formalism of Kirkwood and Buff [58] was used to calculate some thermodynamic properties of ionic solutions. Models I, II and III were compared for pure ammonia and the effects of truncating the multipole series in the model definition for an ammonia molecule were examined. The differences in structure between the models were interpreted in terms of the structure of charge density of an ammonia molecule, and the series of multipoles which approximates that charge density. The axially symmetric quadrupole of ammonia has the effect of reinforcing the negative end of the dipole while partially neutrahzing the positive end, giving the effect of redistributing the positive charge about the "equator" of a molecule and setting up a compromise between the tendency of dipoles to ahgn and quadrupoles to 242 Chapter 4. Summary and Conclusions 243 form T-shaped structures. The effect of adding an octupole to the system was to reinforce the dipolar correlations, resulting in an increased SCMF effective dipole moment, m e, and the increased dipole-dipole energies. The forming of structures more favorable for dipoles was correspondingly less favorable for quadrupoles and was marked by lower quadrupole-quadrupole energies. The larger effective dipole in models containing an octupole resulted in a value for the dielectric constant which was very close to the experimental result, and which had been slightly underestimated by the quadrupole level model. The effect of the change of symmetry from linear to C3,, which occured with the addition of an octupole was found to have only small effects on the structural and thermodynamic properties of pure ammonia. The further addition of a hexadecapole had a negligible effect, partly due to its small absolute magnitude, and partly due to the short-range nature of the added interactions. The reference hard sphere bridge functions needed for use in the RHNC theory for mixtures of large and small particles are not accurately known, so three approximations to the reference system bridge function bn(r) were chosen and compared for systems of macrospheres at infinite dilution in a system of hard spheres. Although all three bridge functions examined were surprisingly different, they all gave very similar hard sphere radial distribution functions gji(r). When these bridge functions were used as reference functions for macrospheres at infinite dilution in Model I ammonia, however, projections of the resulting angular dependent pair distribution functions differed quite considerably when bji(r) was varied. Thus, guaranteeing an accurate hard-sphere limit does not nec-essarily lead to accurate results in these systems. All the bridge functions considered reduced the contact peak in <7oo°m<.(r) calculated by the HNC theory, but different func-tions achieved this result in different ways. The "direct method" bridge function [116], derived directly from a truncated diagrammatic expansion of pair distribution functions in powers of density, was a rnonotonically increasing function. If the bridge function is Chapter 4. Summary and Conclusions 244 considered as a corrective additional (negative) potential, the direct method bR(r) func-tion acts as a short-ranged repulsion which simply decreased the contact peaks. The bridge functions taken from fitting flat-wall [115] and macrosphere [73] Monte Carlo data had more structure and not only added a repulsion at contact, but also an attraction at shght separation, thus both repelling solvent from the wall and pulling it out into the bulk. The range of the structure in the projections of the distribution functions was found to be similar to the range of the structure of the reference bridge functions themselves. The angular probabihty densities were found to be much less sensitive to the reference hard sphere bridge functions, due to their being normalized by the correspond-ing radial distribution function. Thus, the relatively unambiguous angular probabihty distributions can be expected to give a qualitatively accurate picture of the orientational structure of the solvent near a macrosphere surface. The absolute range of the ordered region near the macrosphere surface, however, was somewhat dependent on the choice of reference bridge function. In the absence of any absolute criteria for making a choice, the Lee-Levesque [73] macrosphere reference function was used in most calculations for consistency with the pure solvent and small ion results. The model dependence of results for calculations of ammonia near a large uncharged sphere was examined and the contact values for ammonia molecules near the wall were found to vary inversely with the small changes in the size of the SCMF effective solvent dipole for each model. This reflected the tendency of strong solvent-solvent interactions to pull the solvent away from the surface. As with the pure ammonia calculations, the differences between Models I and II were found to be negligible and Model III only showed significant differences very close to the macrosphere surface where the effects of the solvent-solvent interactions were the most pronounced. The orientational structure of ammonia near the surface of a macrosphere was also examined and the difference in symmetry of the solvent models was found to have very httle effect on the probabihty Chapter 4. Summary and Conclusions 245 distribution for the angle between a molecular dipole vector and a surface normal vector. But, as expected, the probabihty distribution for the angle between an NH-bond vector and a surface normal vector was found to be considerably more pronounced for the models with C3, , symmetry (I and II) than for the model with axial symmetry (III). The dependence on model symmetry was interpreted as a distinct hindrance to rotation about the dipole vector for models with C 3 v symmetry. The effect of macrosphere size on the angular structure was examined, and the most significant changes were found to occur as an uncharged hard sphere was enlarged from lde to 5d„. As the sphere was enlarged further to 20d„ in diameter, the intensity of the peaks in the angular probabihty densities and the local particle density continued to increase forming a highly structured layer of solvent within about 2 solvent diameters of the macrosphere surface. The angular distribution was also seen ,to form several distinct layers parallel to the surface in this region. When orientational and positional probability densities were combined, a layer . *j of about 0Ad„ which had an average density of approximately 1.5 times that of the bulk solvent dominated the resulting functions. The orientational structure in this ordered layer was similar to that of the hydrogen-bonded structure found in frozen ammonia. The next system examined was that of hard sphere ions of diameters ranging between ldB and 20d3 at infinite dilution in Model I ammonia. The unsymmetric charge distribu-tion about an ammonia molecule was found to directly give rise to asymmetric solvation of ions of opposite charges. In order to discover consequences of this asymmetric solva-tion, singly, charged ions at infinite dilution in Model I ammonia were investigated for ion diameters ranging from lds to 20d„ from a number of viewpoints. The interaction energies as well as their contributions from interactions between an ion and specific mul-tipoles were considered as functions of ion diameter. Small cations were characterized by larger interaction energies than were small anions, but as the diameters were increased the anions became preferred over cations. The closely related quantities of contact value Chapter 4. Summary and Conclusions 246 and coordination sphere density were also investigated as a function of ion size. Both of these quantities revealed a denser packing of solvent molecules around a small cation than a small anion, in agreement with the relative size of the interaction energies. As the ion size was increased the density around both cations and anions decreased, but when they grew larger than about 5ds in diameter the density at the ion surface was seen to rebuild, having passed through a minimum value. The reason for these observa-tions was made clear by examining the angular probability distributions of dipole and NH-bond vectors, and was again interpreted in terms of the solvent charge distribution. Around small positively charged ions, the negative end of the dipole, reinforced by the quadrupole, was found to orient itself directly at the ion, which allowed a lower configu-rational energy than for small negatively charged ions, which had to compromise between the orientation preferred by the dipole and that preferred by the quadrupole. The same angular .probability distributions for ammonia around large ions showed results almost identical to those for uncharged macrospheres, indicating that the solvent-solvent inter-action was dominating over the ion-solvent interaction when the electric ion field at the surface was reduced by the macroion size. The dissipating and then rebuilding in coordi-nation sphere density as the ion size was increased was attributed to the initial reduction of ion field effects, and then the solvent-solvent interactions becoming dominant. The macroion-solvent interaction energies were larger for anions than for cations because the surface structure that the solvent preferred was essentially the same as that which small negative ions preferred, and was hence less stable for cations. Two legitimate definitions for single ion partial molar volumes were considered within the framework of the Kirkwood-Buff theory. The definition denoted V° [34], which corresponds to the change in volume when a single ion is added to an infinite bath of solvent, was found to be consistent with the energies, coordination sphere densities and radial distribution function contact values as the ion diameter was varied. The Chapter 4. Summary and Conclusions 247 other definition, denoted V^ °* [34], corresponding to the infinite dilution limit of the volume change when a small amount of solute is added to a constant quantity of solvent (maintaining charge neutrality) was not consistent with V®. We have also considered the behavior of low concentration solutions of electrolytes in liquid ammonia. The concentration dependence of experimentally determined mean molal activity coefficients deviated enormously from the Debye-Hiickel Hmiting law, par-ticularly for AgN03 [6] and 1:2 electrolytes [8]. This behavior was examined closely in terms of ion pairing concepts, and while the shape of the experimental curves can be interpreted to some extent using these ideas, there seems to be no simple explanation for the absolute size of activity coefficients in terms of ion association theories. Furthermore, the fitting of experimental data to equations derived in terms of ion pairs can lead to quite unphysical values for parameters with direct physical interpretation. In particular the parameter for the distance of closest approach of an ion pair was taken to be over 5 solvent diameters [7] which does not agree with the usual ideas of ion pairing. Two definitions of ion pairs were considered in this study and their effects while not neghgible were found to be smaU in comparison to the effects due to ion size. Furthermore, when calculations were corrected for presumably associated ion pairs the result was found to move further away from the Debye-Hiickel Hmiting law. It was concluded that ion pairs in the usual sense can not be the primary cause of the deviations in the activity coefficients from the Hmiting law. It was found that over the range of concentrations given approximately by 2 x 10~4M to 2 x 1 0 _ 2 M no theoretical solutions were possible for any of the models. This phe-nomenon was examined in terms of stability criteria, and these criteria behaved in such a way as would be expected if a phase change were occuring. The phase change in the primitive model (Model V) can be interpreted as the condensation of an ionic gas into an ionic Hquid [17]. A corresponding instabiHty found to exist in a molecular solvent model Chapter 4. Summary and Conclusions 248 (Model III) was interpreted as a phase separation into a relatively dense and relatively dilute phase. While a phase change has been seen in theoretical treatments of the primi-tive model before [20,21,23], here it was demonstrated that (at least at the RHNC level) this phenomenon occurs in molecular solvent models as well. The effect that such a phase change would have on activity coefficients was considered, and it was found that devia-tions in ln y± below the Hmiting law are quahtatively consistent with Hquid-liquid phase separations. The suggestion was made that a phase change could explain the behavior of the experimental activity coefficients. While unproven in real electrolyte-ammonia so-lutions, the existence of phase separation in alkaH metal-ammonia solutions with added salts has been observed [31—33] and lends some support to this suggestion. Bibliography [1] Debye, P., and Hiickel, E . , Phyzik. Z., 2 4 , 185, 305 (1923). [2] Onsager, L., Physik. Z., 2 8 , 277 (1927). [3] Robinson, R. A., and Stokes, R. H., Electrolyte Solutions, 2nd. ed., Academic Press, New York, 1959. [4] Davies, C. W., Ion Association. Butterworth and Co., Washington, 1962. [5] Baldwin, J., and Gill, J. B., Chem. Soc. A, 2040 (1971). [6] Baldwin, J., Evans, J., and Gill, J. B., J. Chem. Soc. A, 3389 (1971). [7] Gill, J. B., and Lowe, B. M., J. C. S. Dalton, 1959 (1972). [8] Evans, J., GiU, J. B., Prescott, A., J. C. S. Faraday I, 1023 (1978). [9] Gans, P., and Gill, J. B., Disc. Faraday Soc, 6 4 , 150 (1977). [10] Bjerrum, N., Kgl Danske Vidensk. Selskab, 9 , 9 (1926). [11] Petrucci, S., Ionic Interactions. Academic Press, New York, 1971. [12] Fuoss, R. M., Trans. Faraday Soc, 3 0 , 967 (1934). [13] Fuoss, R. M., J. Am. Chem. Soc, 8 0 , 5059 (1958). [14] Hansen, J. P., and McDonald, I. R., Theory of Simple Liquids, Academic Press, London, 1976. [15] Friedman, H. L., Ann. Rev. Phys. Chem., 3 2 , 179 (1981). [16] Ramanathan, P. S., and Friedman, H. L., J. Chem. Phys., 5 4 , 1086 (1971). [17] Stell, G., Wu, K. C , and Larsen, B., Phys. Rev. Lett, 3 7 , 1369 (1976). [18] Friedman, H. L., and Larsen, B., J. Chem. Phys., 7 0 , 92 (1979). [19] Gillan, M. J., Molec Phys., 4 9 , 421 (1983). [20] Pitzer, K. S., Chem. Phys. Lett, 1 0 5 , 484 (1984). 249 Bibliography 250 [21] Belloni, L., Phys. Rev. Lett., 57, 2026 (1986). [22] Jedrzejek, C , Konior, J., Streszewski, M., Phys. Rev. A, 35, 1226 (1987). [23] Caccamo, C , Malescio, G., J. Chem. Phys., 90, 1091 (1989). [24] Graham, I. S., and Valleau, J. P., In press. [25] Valleau, J. P., personal communication. [26] Ruff, 0., and Zedner, J., Ber., 41, 1948 (1908). [27] Kraus, C. A., and Lucasse, W. W., J. Am. Chem. Soc, 44, 1949 (1922). [28] Cubicciotti, D. D., J. Chem. Phys., 53, 1302 (1949). [29] Sienko, M. J., J. Am. Chem. Soc, 71, 2707 (1949). [30] Schettler, P. D. Jr., and Patterson, A. Jr., J. Phys. Chem., 68, 2865 (1964). [31] Schettler, P. D. Jr., and Patterson, A. Jr., J. Phys. Chem., 68, 2870 (1964). [32] Doumaux, P. W., and Patterson, A. Jr., J. Phys. Chem., 71, 3535 (1964). [33] Doumaux, P. W., and Patterson, A. Jr., J. Phys. Chem., 71, 3540 (1964). [34] Kusalik, P. G., and Patey, G. N., J. Chem. Phys., 89, 5843 (1988). [35] Carnie, S. L., and Patey, G. N. , Molec Phys., 47, 1129 (1982). [36] Kusalik, P. G., and Patey, G. N., J. Chem. Phys., 79, 4468 (1983). [37] Torrie, G. M., Kusalik, P. G., and Patey, G. N., J. Chem. Phys., 88, 7826 (1988). [38] Kusalik, P. G., and Patey, G. N., J. Chem. Phys., 88, 7715 (1988). [39] Caillol, J. M. , Levesque, D., Weis, J. J., Perkyns, J. S., Patey, G. N., Molec Phys., 62, 1225 (1987). [40] Narten, A. H., J. Chem. Phys., 66, 3117 (1977). [41] Baldwin, J., and Gill, J. B., Physics and Chemistry of Liquids, 2, 25 (1970). [42] McMillan, W. G., and Mayer, J. E., J. Chem. Phys., 13, 276 (1945). [43] Kusalik, P. G., Patey, G. N., J. Chem. Phys., 89, 7478 (1988). Bibliography 251 [44] Friedman, H. L., and Dale, W. D. T., Modern Theoretical Chemistry, Volume 5, edited by B. J. Berne, Plenum, New York, 1977. [45] Krishnan, V., Friedman, H. L., J. Solution Chem., 3, 727 (1974). [46] Lee, L. Y., Fries, P. H., and Patey, G. N., Molec. Phys., 55, 751 (1985). [47] Perkyns, J. S., Fries, P. H., and Patey, G. N., Molec. Phys., 57, 529 (1986). [48] Fries, P. H., and Patey, G. N., J. Chem. Phys., 82, 429 (1985). [49] Ornstein, L. S., and Zernike, F., Proc. K. Ned. Akad. Wet. Amsterdam, 17, 793 (1914). [50] Lebowitz, J. L., and Percus, J. K., Phys. Rev., 144, 251 (1966). [51] Percus, J. K., and Yevick, G. J., Phys. Rev., 110, 1 (1958). [52] Patey, G. N., Levesque, D., Weis, J. J., Molec. Phys., 38, 1635 (1979). [53] Levesque, D., Weis, J. J., Patey, G. N., Molec. Phys., 51, 333 (1984). [54] Lado, F., Phys. Rev. A., 135, 1013 (1964). [55] Lado, F., Molec. Phys., 31, 1117 (1976). [56] Perera, A., Kusahk, P. G., Patey, G. N„ Molec. Phys., 60, 77 (1987). [57] Rossky, P. J., Dudowicz, J. B., Tembe, B. L., Friedman, H. L., J. Chem. Phys., 73, 3372 (1980). [58] Kirkwood, J. G., and Buff, F. P., / . Chem. Phys., 19, 774 (1951). [59] Kusahk, P. G., "The Structural, Thermodynamic and Dielectric Properties of Electrolyte Solutions: A Theoretical Study", Ph. D. Thesis, University of British Columbia, 1987. [60] Friedman, H. L., A Course in Statistical Mechanics, Prentice-Hall, Englewood Cliffs, N. J., 1985. [61] Messiah, A., Quantum Mechanics, Wiley, New York, 1958. [62] Rotenberg, M. , Bevins, R., Metropolis, N., and Wooten, N., The 3-j and 6-j Symbols, Massachusetts Institute of Technology, Cambridge, 1959. [63] Kusahk, P. G., "A Theoretical Study of Dilute Aqueous Electrolyte Solutions", M. Sc. Thesis, University of British Columbia, 1984. Bibliography 252 [64] Perkyns, J. S., "The Solution to the Reference Hypernetted-chain Approximation for Fluids of Hard Spheres with Dipoles and Quadrupoles with Application to Liquid Ammonia", M. Sc. Thesis, University of British Columbia, 1985. [65 [66 [67; [68 [69 [70 [71 [72 [7.3 [74 [75 [76 [77 [78 [79 [80 [81 [82 [83 Carnie, S. L., Chan, D. Y. C , and Walker, G. R., Molec. Phys., 43, 1115 (1981). Hirschfelder, J. 0., Curtiss, C. F., and Bird, R. B., Molecular Theory of Gases and Liquids, Wiley, New York, 1954. Stone, A. J., Molec. Phys., 29, 1461 (1975). Price, S. L., Stone, A. J., and Alderton, M., Molec. Phys., 52 , 987 (1984). Rushbrooke, G. S., and Scoins, H. I., Proc. Roy. Soc. A., 216, 203 (1953). Croxton, C. A., Introduction to Liquid State Physics, John Wiley and Sons, Lon-don, 1975. Croxton, G. A., Liquid State Physics, Cambridge University Press, London, 1974. Verlet, L., and Weis, J. J., Phys. Rev. A., 5, 939 (1972). Lee, L. L., and Levesque, D., Molec. Phys., 26 , 1351 (1973). Blum, L., and Torruella, A. J., J. Chem. Phys., 56, 303 (1972). Blum, L., J. Chem. Phys., 57, 1862 (1972). Blum, L., J. Chem. Phys., 58, 3295 (1973). Abramowitz, M. , and Stegun, I. A., eds., Handbook of Mathematical Functions, Dover, New York, 1970. Edmonds, A. R., Angular Momenta in Quantum Mechanics, Princeton, New Jer-sey, 1960. Frolich, H., Theory of Dielectrics, 2nd. ed., Oxford University Press, Oxford, 1958. Kirkwood, J. G., J. Chem. Phys., 4, 592 (1936). Stell, G., Patey, G. N., and H0ye, J. S., Adv. Chem. Phys., 38 , 183, (1981). Patey, G. N. , Molec. Phys., 34, 427 (1978). H0ye, J. S., and Stell, G., J. Chem. Phys., 64 , 1952 (1976). Bibliography 253 [84] Levesque, D., Weis, J. J., and Patey, G. N., J. Chem. Phys., 72, 1887 (1980). [85] Stillinger, F. A., and Lovett, R., / . Chem. Phys., 49, 1991 (1968). [86] Chan, D. Y. C , Mitchell, D. J., Ninham, B. W., and Pailthorpe, B. A., J. Chem. Phys., 69, 691 (1987). [87] Wolynes, P. G., Ann. Rev. Phys. Chem., 31, 345 (1980). [88] Falkenhagen, H., Electrolytes, trans: R. P. Bell, Clarendon Press, Oxford, 1934. [89] Hubbard J. B., Colonomos, P., and Wolnes, P. G., J. Chem. Phys., 71, 2652 (1979). [90] Barker, J. A., and Henderson, D., Rev. Mod. Phys., 48, 587 (1976). [91] Steele, W. A., J. Chem. Phys., 39, 3197 (1963). [92] Perkyns, J. S., Kusahk, P. G., and Patey, G. N., Chem. Phys. Lett, 129, 258 (1986). [93] Wertheim, M. S., Ann. Rev. Phys. Chem., 30, 471 (1979). [94] Gray, C. G., and Gubbins, K. E., Theory of Molecular Fluids, Oxford University Press, Oxford, 1984. [95] Kusahk, P. G., and Patey, G. N., J. Chem. Phys., 86, 5110 (1987). [96] Friedman, H. L., and Ramanathan, P. S., J. Chem. Phys., 74, 3756 (1970). [97] Harned, H. S. and Owen, B. B., Physical Chemistry of Electrolyte Solutions, 3rd. ed., Reinhold Publishing Corp., New York, 1959. [98] Prigogine, I., and Defay, R., Chemical Thermodynamics, translated by D. H. Ev-erett, Longmans Green and Co., London, 1954. [99] Tisza, L., Generalized Thermodynamics, M. I. T. Press, Cambridge, 1966. [100] Ursenbach, C. P., and Patey, G. N., to be published. [101] Millero, F. J., Chem. Rev., 71, 147 (1971). [102] Zana, R. F., and Yeager, E., J. Phys. Chem., 71, 521 (1967). [103] Denison, J. T., and Ramsey, J. B., J. Am. Chem. Soc, 77, 2615 (1955). [104] Talman, J. D., J. Comput. Phys., 29, 35 (1978). Bibliography 254 [105] Rossky, P. J., and Friedman, H. L., J. Chem. Phys., 72, 5694 (1980). [106] Caillol, J. M. , Weis, J. J., J. Chem. Phys., 90, 7403 (1989). [107] Springer, J. F., Pokrant, M. A., and Stevens, F. A., J. Chem. Phys., 58, 4863 (1973). [108] Vancini, C. A., Synthesis of Ammonia, Macmillan, New York, 1971. [109] McClellan, A. L., Tables of Experimental Dipole Moments, W. H. Freeman, San Francisco, 1963. [110] Kukolich, S. G., Chem. Phys. Lett, 5, 401 (1970). [Ill] Kukolich, S. G., Chem. Phys. Lett., 12, 216 (1971). [112] Amos, R. D., Handy, N. C., Knowles, P. J., Rice, J. E., and Stone, A. J., J. Phys. Chem., 89, 2186 (1985). [113] Bridge, N. J., and Buckingham, A. D., Proc. Roy. Soc. A., 205, 135 (1966). [114] Haar, L., and Gallagher, J. S., J. Phys. Chem. Ref. Data, 7, 635 (1978). [115] Henderson, D., and Plischke, M., Proc. Roy. Soc. A., 400, 163 (1985), [116] Attard, P. and Patey, G. N., to be published. [117] Snook, I. K., and Henderson, D., J. Chem. Phys., 68, 2134 (1978). [118] Frank, H. S., "Solvent Models and the Interpretation of Ionization and Solvation Phenomena", Chemical Physics of Ionic Solutions, Conway, B. E., and Barradas, R. G., eds., John Wiley and Sons, New York, 1966. [119] Huheey, J. E., Inorganic Chemistry: Principles of Structure and Reactivity, 3rd. ed., Harper and Row, New York, 1983. [120] Kirkwood, J. G., Theory of Solutions, Salzburg, Z. W., ed., Gordon and Breach, New York, 1968.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A theoretical study of ammonia-salt mixtures in bulk...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A theoretical study of ammonia-salt mixtures in bulk solutions and interfacial regions Perkyns, John Stephen 1990
pdf
Page Metadata
Item Metadata
Title | A theoretical study of ammonia-salt mixtures in bulk solutions and interfacial regions |
Creator |
Perkyns, John Stephen |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | Five models, ranging from a full molecular polar/polarizable model with C₃v symmetry, to the primitive model of ions in a dielectric continuum, have been used to study the properties of ammonia both as a pure liquid and as a solvent. Ammonia is modelled as a multipolar polaxizable hard sphere and ions as charged hard spheres. Using these models, in conduction with the Reference Hypernetted-chain integral equation theory, ammonia has been studied as a pure liquid, as a solvent near charged and uncharged surfaces, and in electrolyte solutions of finite concentration. The formalism of Kirkwood and Buff was used to obtain thermodynamic quantities of ionic solutions from calculated distribution functions. Structural, thermodynamic and. dielectric properties were calculated for pure ammonia and were compared with experiments where possible. Values for the dielectric constant and the configurational energy were found to compare favorably with experiment. Ammonia formed a relatively dense, highly structured layer within two solvent diameters of an uncharged surface. This structure was analyzed in terms of angular probability distributions of the molecular dipole vector and the NH-bond direction, and was found to be similar to that of frozen ammonia. The extreme asymmetry of solvation of unlike charges in ammonia was also investigated. Small cations were found to be more favorably solvated than small anions, but as the ion size was increased, the situation reversed. Estimates of the number of ion pairs in liquid ammonia and their effects on the behavior of mean molar activity coefficients were examined. Large differences between experimental activity coefficients and the Debye-Hückel hmiting law could not be explained by the usual ideas of ion pairing. It was found that the integral equation theories appear to have no solution between ionic concentrations of about 2 x 10⁻⁴M and 2 x 10⁻²M for either molecular or continuum models. Using rigorous stability criteria, this behavior was shown to be consistent with the onset of a phase change. It is proposed that such phase separation phenomena might explain the unusual behavior of the experimental activity coefficients. |
Subject |
Ammonia -- Analysis Ammonium salts |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-01-12 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0059785 |
URI | http://hdl.handle.net/2429/30588 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1990_A1 P39.pdf [ 8.71MB ]
- Metadata
- JSON: 831-1.0059785.json
- JSON-LD: 831-1.0059785-ld.json
- RDF/XML (Pretty): 831-1.0059785-rdf.xml
- RDF/JSON: 831-1.0059785-rdf.json
- Turtle: 831-1.0059785-turtle.txt
- N-Triples: 831-1.0059785-rdf-ntriples.txt
- Original Record: 831-1.0059785-source.json
- Full Text
- 831-1.0059785-fulltext.txt
- Citation
- 831-1.0059785.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0059785/manifest