S T A T I S T I C A L M E C H A N I C S O F S O L V A T I O N O F M A C R O P A R T I C L E S B y Reka Z. Vasarhelyi B . Sc. (Chemistry) University of Bri t ish Columbia , 1990 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S C H E M I S T R Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 1992 © Reka Z. Vasarhelyi, 1992 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 Chemistry The University of British Columbia Vancouver, Canada Date Oct. 13 1992 DE-6 (2/88) Abstract Integral equation theories are employed to study the solvation of large particles. The solutions studied consist of spherically symmetric solvent particles which interact through hard sphere and Lennard-Jones potentials, and hard sphere solutes. In particular, the Ornstein-Zernike (O-Z) equation is solved with the hypernetted-chain ( H N C ) closure, to obtain the pair correlation functions of mixtures at infinite di lut ion. The pair correlation functions in the O-Z equation and the H N C closure are expressed as power series in solute density to yield a pair of coupled equations which determine the derivatives with respect to solute density of the solvent-solvent pair correlation func-tions. The latter describe the perturbation of the solvent upon addition of a single solute particle. The derivatives are analysed to yield components that scale as the volume and surface area of the macroparticle, and which are then identified as changes in solvent structure due to the presence of a finite size particle and a flat surface respectively. From the pair correlation functions and their derivatives the excess internal energy, Helmholtz free energy, and entropy of solvation are calculated. The excess quantities are also separated into contributions from finite size and surface effects. Both components of the excess internal energy are negative at low densities, and become less negative for high density liquids. The magnitude and sign of the two contributions to the energy depend on physical conditions such as temperature and pressure. The excess entropy of solvation is negative for all systems studied, indicating that introduction of a macroparticle or a flat surface increases the spatial ordering of a bulk l iquid. n Table of Contents Abstract ii List of Tables v List of Figures vi Acknowledgement viii 1 Introduction 1 2 Theory 8 2.1 The Radial Distribution Function 8 2.2 Choice of Model 9 2.3 Thermodynamics 11 2.4 Integral equations 15 2.4.1 The Ornstein-Zernike equation 15 2.4.2 Choice of Closure 18 2.4.3 Numerical Details 22 3 Results and Discussion 23 3.1 Input Parameters 23 3.2 Radial Distribution Functions 24 3.3 Thermodynamics 42 4 Conclusions 53 i i i Appendices 56 A Helmholtz Free Energy for a Mixture 56 Bibliography ^QQ iv List of Tables 3.1 Input Parameters; H S / H S refers to a hard sphere solvent with a hard sphere solute; I I S / L J refers to a Lennard-.Jones solvent with a hard sphere solute 23 3.2 Excess Energy of Solvation from extrapolation of < Uex > /N^d2 vs. d plots 47 3.3 Excess Energy of Solvation from integration of <*)</!(r)and 8g2(r) 47 3.4 Excess Energy of Solvation at constant total density from extrapolation of < Uex > /N^P vs. d plots 49 3.5 Excess Energy of Solvation at constant total density, from integration of (fyl(r)and 8g2{r) 49 3.6 Excess Free Energy of Solvation, from Eq . 2.26 50 3.7 Excess Free Energy of Solvation, from Eq . 2.61 50 3.8 Excess Entropy of Solvation 52 List of Figures 3.1 Solute-solvent radial distribution functions for H S / H S solutions; p* = 0.8. The solid, dotted and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively 25 3.2 Solute-solvent radial distribution functions for H S / L J solutions; p* = 0.6, T* = 1.35. The solid, dotted, and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively 26 3.3 Solvent-solvent radial distribution function for H S / L J solution; p* = 0.6, T* = 1.35 28 3.4 Derivatives of solvent-solvent radial distribution functions for H S / H S so-lutions; p* = 0.8. The solid, dotted, and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively 29 3.5 Derivatives of solvent-solvent radial distribution functions for H S / L J so-lutions; p* = 0.6, T* = 1.35. The solid, dotted, and dashed curves corre-spond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively. 30 3.6 8gvv(r) at r = Id for the H S / L J solution; p* = 0.6, T* = 1.35 31 3.7 Sg^(r) at r = 5rf for the H S / L J solution; p" = 0.6, T* = 1.35 32 3.8 Sg„v(r) at r = lOd for the H S / L J solution; p* = 0.6, T* = 1.35 33 3.9 6gvv{r) at r = 30rf for the H S / L J solution; p" = 0.6, T* = 1.35 34 3.10 Sgl(r) for H S / L J solution; p* = 0.6, T* = 1.35. The solid, dotted, and dashed lines correspond to solutions of Eqs. 3.61 using different solute diameter pairs 36 vi 3.11 8g2(r) for H S / L J solution; p' = 0.6, T" = 1.35. The solid, dotted, and dashed lines correspond to solutions of Eqs. 3.61 using different solute diameter pairs 37 3.12 6gl(r) for H S / L J solution on an expanded scale; p* = 0.6, T* = 1.35. . . 38 3.13 8g2(r) for H S / L J solution on an expanded scale; p* = 0.6, T* = 1.35 . . . 39 3.14 8gvv(r) at constant total density for H S / L J solution; p = 0.6, T* = 1.35. The solid, dashed, and dotted curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively 40 3.15 Excess Energy of Solvation as a Function of Solute Diameter 43 3.16 Excess Free Energy of Solvation as a Function of Solute Diameter 44 3.17 Excess Entropy of Solvation as a Function of Solute Diameter 45 VI1 Acknowledgement I would like to acknowledge the help and support I have received from several people during the course of this work. M y supervisor Dr. G . N . Patey suggested an interesting topic, and provided expert guidance throughout. My fellow workers contributed to an excellent work environment under the palm tree. Discussions with Dr . Phi l ip At ta rd and Daniel Berard were especially helpful and enjoyable. Finally, I would like to thank my husband Paul , who was, as always, understanding and supportive of my efforts. v 111 Chapter 1 Introduction A great deal of recent work in l iquid state theory has focussed on inhomogeneous liquids. Such systems include dense fluids in contact with a wall or macroparticle, and liquids in slits. The elucidation of the microscopic structure of a l iquid in contact with a macro-scopic, body is important for the understanding of electrode processes, the phenomenon of wetting, and interactions between large particles. The general purpose of any statistical mechanical theory is the development of micro-scopic models which reproduce macroscopic properties such as thermodynamic quantities. A theoretical treatment of the l iquid state is unique in the sense that a simple reference system, such as the harmonic solid and ideal gets, for solids and gases respectively, does not exist. Unlike a solid or gas, a real l iquid cannot be thought of as a departure from an ideal state. In this thesis, only classical fluids are considered. This means that the vari-ables of kinetic energy, velocity, and momentum, are decoupled from those of potential energy, position, and orientation. Since the kinetic contribution to thermodynamic quan-tities is known, only the contribution of the potential energy will be studied here. Also, we restrict ourselves to simple liquids. Specilically, we mean dense fluids of monatomic, spherical, argon-like particles. Central to the microscopic description of liquids of simple spherical particles is the radial distribution function, g(r). This function, which represents the probability of finding a second particle at a distance r from a particle at the origin, is particularly useful since it can be related to the structure factor. The latter is a physical observable 1 Chapter 1. Introduction 2 of X-ray and neutron scattering experiments. Furthermore, if we make the simplifying assumption that the total potential energy of the system is a sum of pairwise additive potentials, then all thermodynamic quantities can be expressed as integrals over the radial distribution function. Examples of two such relations are the compressibility equation, which is also known as the vir ia l equation of state [1]. If an exact g(r) were available, both equations would give the same value for the pressure. For approximate radial distribution functions, the discrepancy between the pressures calculated from the two equations is an indication of the thermodynamic consistency of the approximate theory [1]. The function g(r) can be calculated using two different, but complementary, tech-niques, computer simulations and analytical theories [2]. Two types of computer ex-periments are generally performed for liquids. A Molecular Dynamics ( M D ) simulation involves the solution of the classical equations of motion for a set of particles, assigned ini t ial velocities in accordance with the Maxwell -Bol tzmann distribution for a given tem-perature. This technique is appropriate for the study of time-dependent processes. The Monte Carlo ( M C ) method takes a probabilistic, rather than a deterministic approach. From an init ial set of randomly assigned coordinates, further configurations are generated by random displacements. The configuration is accepted or rejected depending on the Bol tzmann factor, and therefore the probability of such an arrangement of the particles of a system at equilibrium. A M C simulation yields equil ibrium ensemble averages. Since particles in a dense fluid are influenced by many other particles, the interaction between two molecules is a sum of pair and many-body forces. For argon-like liquids, potentials (1.1) and the pressure equation, (1.2) Chapter 1. Introduction 3 which accurately fit experimental data have been developed by Barker [3], using quantum mechanical perturbation methods. In most cases however, simple analytic expressions for the effective interaction between two particles are employed. A pair potential is assumed which retains the main features of the interparticle interaction. For simple liquids, the dominant features are a short-range strong repulsion and a long-range weak attraction. The potential energy of the model fluid is then taken to be a sum of pair interactions over al l pairs of particles. In a computer simulation, the form of the pair-potential is part of the input, and is therefore defined exactly. That is, the results of a computer simu-lation are essentially exact for the model l iquid. Analyt ical theories, on the other hand, involve further approximations. When comparing the results of analytical calculations to experimental data, the source of deviations is ambiguous. In contrast, comparison with computer simulation tests the mathematical approximations of the theory. Current analytical theories can be placed in two broad categories, perturbation and integral equation theories. The former, due originally to the work of Zwanzig [4], treats a fluid whose molecules interact v ia a potential containing both a repulsion and an attraction as a perturbation of a fluid with a purely repulsive pair potential. Whi le a reference system consisting of hard spheres is not an exact ideal state, it has proven to be a practical reference state. This is due to the fact that the structure of a fluid is determined almost entirely by the repulsive part of the potential. Also, hard sphere liquids have been completely characterised by computer simulations. Perturbation theories are easy to work out, but their success has been l imited to simple systems [5]. Integral equations are somewhat more demanding, but have been more successful for complex systems. Integral equation theories are based on approximate equations for the pair distribution function. One of the first such relations was developed by Kirkwood in the 1930's [6]. The differentiation of the equilibrium n-particle density with respect to a coupling parameter gives the Ki rkwood hierarchical equations. If the differentiation is Chapter 1. Introduction 4 carried out with respect to the coordinates of a given particle, another set of equations, called the Born-Green-Yvon hierarchy is obtained. As the names suggest, both theories consist of sets of coupled equations which give the n-particle distribution functions in terms of the (n+l)-particle distribution functions. Another set of integral equtions, which was derived in the 1950's, is based on a simple physical idea, proposed by Ornstein and Zernike several decades earlier [7]. The mathematical statement of their idea is, This equation, and its relation to g(r), will be discussed in greater detail in the next section. For now, it is sufficient to say that the total influence of particles 1 and 2 on each other, given by the total correlation function / / ( r i , ! ^ ) , is a sum of a direct part described by the direct correlation function c ( r i , r2 ) , and an indirect part described by the integral term in E q . 1.3. The quantity p represents the bulk density. For a homogeneous isotropic fluid, the interparticle forces are central; therefore, the correlation functions depend only on the magnitude of the separation between particles [7]. In such a case, the Ornstein-Zernike (O-Z) equation contains only two unknowns, namely h(r) and c(r). Hence, given another relation between the total and direct corre-lation functions, a complete solution can be obtained. Approximate equations relating h(r) and c(r), such as the Perçus-Yevick ( P Y ) , mean spherical ( M S A ) , hypernetted-chain ( H N C ) , and reference hypernetted-chain ( R H N C ) approximations, combined with the O-Z equations form several integral equation theories for g(r)-The homogeneous integral equation theories have been extended to the case of a macroparticle in contact with a bulk fluid, by considering a l imi t ing case of a homoge-neous two component mixture [8]. In the l imit of large solute radius and infinite solute dilution, the mixture becomes a single component fluid in contact with a flat surface [9]. (1.3) Chapter 1. Introduction 5 Perram and Whi te [8], have shown that it is physically reasonable to take this l imi t by demonstrating that for a growing hard sphere solute, g(r) becomes vir tually independent of the size of the macroparticle for solute sizes as small as fifteen times the solvent diam-eter. In a concurrent but independent study Henderson et cd. [9] evaluated analytically the infinite radius l imit of the 0 - Z equation. They combined this result with the P - Y approximation to obtain a one dimensional wall-solvent distribution function. The wall-solvent 0 - Z equation, in combination with various closure relations, has been used to calculate the density profile of systems such as hard sphere and Lennard-Jones fluids in contact with hard or attractive walls [10]. The technique has also been successfully applied to the problem of developing a molecular description of the electrical double layer [11, 12]. For liquids near a phase transition, density profiles calculated from the homogeneous 0 - Z equation do not agree well with those obtained from simulations [3], The inho-mogeneous O-Z equation has been applied to such cases [13]. The inhomogeneous O-Z equation is obtained if the bulk density in Eq . 1.3 is replaced with the single particle den-sity pfa), and is moved inside the integral. This means that />(r3) is another unknown, and an additional relation between the single particle density and the pair correlation function is required to determine the variables completely. Although for liquids near a phase transition improved results are obtained from the inhomogenous equations, the implementation of the technique is significantly more difficult. Both versions of the O-Z equations are definitions for the direct correlation function, and are therefore formally exact. W i t h improved approximate closures, the homogeneous O-Z equation is expected to be useful lor systems near phase transitions, while maintaining computational simplic-ity [14]. In this thesis, the radial distribution functions of homogeneous liquids in contact with a hard macroparticle are calculated using the homogeneous O-Z integral equation with the H N C closure. Chapter 1. Introduction 6 As indicated previously, the thermodynamic quantities of the system can be evaluated from integrals over the radial distribution function. In particular, we calculate the excess energy, Helmholtz free energy, and entropy of solvation of a single macroparticle. Here, the term "excess" refers to the difference between the thermodynamic quantity of a mixture and that of the pure solvent. In addition, we separate the excess quantities into two components; one due to the insertion of a Unite size particle into a l iquid , and the other due to the rearrangement of a solvent near a hard surface. The latter contribution is then associated with surface excess thermodynamic quantities. To obtain the excess internal energy at infinite dilution, we cannot simply evaluate the sum of pair energies over all solute-solvent interactions. The presence of a single solute particle perturbs the solvent. The structural changes in the solvent give rise to an additional contribution to the internal energy, which can be expressed as an integral over the derivative with respect to solute density of the solvent-solvent radial distribu-tion function. Following Garisto ct.al. [15], the O-Z equation for the derivatives of the distribution functions are derived, and solved with the H N C closure. Thereby, the excess internal energy at infinite dilution is obtained. The excess Helmholtz free energy associated with adding a single particle to the bulk fluid is evaluated using a coupling integration technique [7]. A coupling parameter is defined in such a way that by varying it the solute is effectively added to or removed from solution. For a hard solute, the infinite dilution l imit of the free energy is evaluated directly by applying some of the ideas of scaled particle theory [17]. The coupling pa-rameter which multiplies the diameter of the hard particle, physically scales in a single solute. The change in free energy due to this scaling process is the quantity sought. A n expression for the excess Helmholtz free energy of a mixture is obtained, and subse-quently, following Karplus et.al. [18], its infinite dilution l imit is evaluated. Final ly, the excess entropy of solvation is obtained from the thermodynamic expression relating the Chapter 1. Introduction 7 internal energy, Helmholtz free energy, and entropy. The remainder of this thesis is organized as follows. In Chapter 2, a detailed derivation of the equations required to calculate the density profiles and thermodynamic quantities are presented. Chapter 3, contains results obtained for a hard sphere solute in a hard sphere or Lennard-.Iones solvent at a number of different solvent densities and tempera-tures. In Chapter 4, results and conclusions are summarized. Chapter 2 Theory 2.1 The Radial Distribution Function We consider a system consisting of particles ( 1 , N ) with spatial coordinates ( r i , ...r/v)-A distribution function expresses the probability of various configurations of a subset (1, . . . ,n) of these particles. The thermodynamic relationships in this chapter are derived in the canonical, ( N V T ) ensemble. The choice of ensemble is somewhat arbitrary, as for the calculation of thermodynamic quantities the various ensembles are indistinguish-able [6]. Hence, the ensemble is chosen on the basis of computational convenience of thermodynamic quantities. We consider fluids of particles which interact through pairwise additive potentials. Given an N-particle system, pairwise additivity is equivalent to the condition where UN(TI , r ^ ) is the total potential energy, u(vij) is the pair potential, and the summation is over all pairs of particles. For such a system the probability that particle 1 is in volume element r / r i , particle 2 in dr? , and particle n in volume element drn, irrespective of the positions of the remaining particles is given by [2] UN(ru...,rN) = 5Z«(rtj) , (2.4) N\ f... fv exp[-[1UN]drn+l...drN (N-n)l ZN (2.5) where the configurational integral is defined as (2.6) Chapter 2. Theory 9 The factor /? equals j^f, where is the Boltzniann constant. It can be shown, that in a bulk homogeneous fluid, is equivalent to the bulk density, p. For an isotropic fluid, p^ depends only on the separation between particles 1 and 2, and wil l be denoted ^ 2 ' ( r r 2 ) , or when no confusion is possible, simply p(2\r). If the probability of every particle i being in dvi at r,, were independent of the probability of every other particle j being in dtj at r-j, then p^ would be the product of single particle distributions, p^K When this independence condition does not hold, a correlation function can be defined as / > ( , l ) ( n , - , r w ) = / ) ( 1 ) ( r 1 ) . . y i ' ( r n ) / , ) ( r i , . . . , r n ) . (2.7) The function c / n ' ( r i , r , J is a factor that measures how much p ^ ( r - i , r n ) deviates from the product of independent probabilities. Taking n = 2 in E q . 2.7, and noting that in a fluid, g^(r) is again a function only of the interparticle separations 7-12, we obtain the radial distribution function </ 2 ' ( r i 2 ) , which lias a simple physical meaning; it is the factor which multiplies the bulk density, to give the local density. That is, given a particle at the origin the number of particles in a spherical shell of thickness dr and radius ?• centered at the origin is pg 2.2 Choice of Mode l The dominant features of pair potentials of simple liquids are a short-range strong repul-sion, and a long-range weak attraction. The repulsion which is an expression of the fact that two particles cannot completely overlap, is a dominant contribution to the short-range order exhibited by dense fluids. The attractive forces increase the stability of the fluid. From quantum mehanical perturbation theory, we know that at large separations the rare gas pair potential has the form [3] . u(r) = -cG/r6 - c s / r 8 + ... . (2.8) Chapter 2. Theory 10 The exact form of the repulsive forces is not known. Although, in principle, an exact pair potential can be derived using quantum mechanics, the problem is too difficult numerically to be practical. Instead, simple analytical expressions with parameters which can be adjusted to fit experimental data are employed. The simplest of these is the hard sphere potential, where d is the hard sphere diameter. The hard sphere potential does not include attractive forces. Since the Boltzmann factor e.x\)(—fiUN{TN)) is either zero or infinity for the potential defined by Eq . 2.9, temperature has no effect on the average configurational properties of a hard sphere fluid. Despite its simplicity, a hard sphere l iquid exhibits many of the structural properties of real liquids. The Lennard-Jones(L-J) potential, incorporates the asymptotic behaviour of the true potential, and replaces the repulsion with an inverse power law, to obtain the mathemat-ically simple expression, with cr and e as the adjustable parameters. This potential has been used extensively for the study of argon-like liquids. The underlying assumption that the total potential energy can be expressed as a sum of pairwise additive potentials means that all many-body interactions are ignored. For rare-gas atoms, the net effect of the dominant many-body forces is repulsive [2]. We also know from experiment that the true potential has a deeper well and a weaker tail than the L-. l function [2]. Nevertheless, probably as a result of appropriate cancellation of errors, the L-.J potential has been reasonably successful in accounting for man}' of the properties of argon-like liquids [2]. (2.9) (2.10) Chapter 2. Theory 11 2.3 Thermodynamics In this section, thermodynamic quantities in terms of radial distribution functions for a two component system wi l l be presented. For classical liquids, all thermodynamic quantities can be separated into ideal and non-ideal parts. The non-ideal term contains the contribution of interparticle interactions. The terms energy, Helmholtz free energy, and entropy wi l l refer only to this second term, as it is the only contribution that wi l l be considered here. The term "excess" wil l be used in the context of differences from the pure solvent. In particular, we calculate the change in energy, Helmholtz free energy, and entropy associated with taking a single solute particle from the ideal gas phase into solution in a pure solvent [15]. For the derivation of the excess potential energy, we appeal to the physical interpre-tation of g(r). Consider a two component mixture consisting of a solvent v and a solute fi. Around a central molecule of type //, the number of particles of the same component in the volume element A-Kr2dr is given by where pu = p \ y i p is the bulk density, and \ „ is the mole fraction, jf. Therefore, the potenial energy due to all solvent-solvent interactions in the mixture is The factor N„ accounts for all solvent particles and the | ensures that each interaction is counted only once. Similarly, the contributions to the total potential energy from solvent-solute and solute-solute interactions are Pv<jvv(r)Anr2dr, (2.11) (2.12) (2.13) Chapter 2. Theory respectively. The sum of these contributions divided by N is the potential energy per particle for the mixture The first term in E q . 2.18 represents the sum of all solute-solvent interactions. The second term is the contribution to the excess energy from the change in solvent structure induced by the presence of the solute. The next quantity sought is the change in Helmholtz free energy associated with adding a single solute to a pure solvent. An expression for the Helmholtz free energy in terms of the radial distribution function can be formulated via a coupling parameter integration technique [7]. A particle can be formally added to the system if the energy is expressed in terms of a coupling parameter £, as in Chapter 2. Theory 13 The parameter £ is defined so that u(rij;0) = 0 and u(rij;a) = u(vij), where a is a number signifying complete coupling. The value of a depends on the way we choose to introduce the effect of the particle. A solute can be effectively added to or removed from solution by physically scaling in its diameter, or by increasing the magnitude of its interaction potential from zero. The configurâtional integral, in terms of the coupling parameter is ZN{0 = J exp[-(WN(Ç)]dr1...drN . (2.20) From the fundamental relation connecting thermodynamics to canonical ensemble aver-ages, we have the configurational Helmholtz free energy as a function of £ [7], 0<A(Q>=-\n(?0) . (2.21) In order to calculate the free energy in excess of the pure solvent, the derivative of A with respect to £ is integrated, according to 0<A">= [ l d t d l i < m > . (2.22) Jo o( A t this point, it is convenient to introduce the function yu^{i"i 0 = exp[/3u(r; £)]flW( r ' which is continuous for all potentials. If the configurational integral is differentiated ex-pl ic i t ly with respect to £, an expression for the excess free energy in terms of £ and g(r;£) is obtained, f j < ^ > = ~ J*W*P J~ ^M-0u*Ar\mv**(rUydr. (2.23) The factor \/Nlt is impl ic i t ly included, since by coupling a single particle, the excess free energy per solute particle is computed. For hard solutes it is useful to define a coupling parameter which physically scales in a particle by increasing its diameter. This allows us to easily add particles of various sizes, and examine the free energy as a function of increasing solute diameter. Following the Chapter 2. Theory 14 approach of Reiss et.al. [17], the coupling parameter is defined in terms of the interaction of the hard sphere solute with the solvent, Here, dv is the solvent diameter and d is the variable representing the increasing diameter of the solute. The particle is completely coupled when d = d^ and £ - - (du + d^)/2 . The definition above implies that when the solute is completely uncoupled, its diameter is —dv. Whi le seemingly unphysical, this in fact is necessary, because the effect of the particle does not disappear when the diameter is zero. A n impenetrable point particle would sti l l be present when d = 0, which would exclude solvent particles from a volume element equal to the volume of a solvent particle. When d = —d„, the effect of the solute vanishes, and the particle is completely removed. Once the potential is defined explicit ly as a function of the coupling parameter, the derivative in the integrand of Eq . 2.23 can be evaluated. For the potential and coupling parameter defined in Eq . 2.24, the expression to be differentiated is The function in Eq . 2.25 is equivalent to 1 — H(Ç — r ) , where / / ( £ — is the Heav-iside step function, and its derivative with respect to £ is the Dirac delta function, —8(£ — r). Therefore the integration with respect to 7- becomes t r iv ia l , and the one-dimensional integral is obtained. To evaluate the free energy for a hard sphere solute, then, requires an integration over the contact values of </(/•), for solute particles of increasing size. (2.24) (2.25) (2.26) Chapter 2. Theory 2.4 Integral equations The final piece of information required for the calculation of thermodynamic quantities is the radial distribution function. In order to find an approximate g(r) we solve the Ornstein-Zernike (O-Z) integral equation with the hypernetted-chain ( H N C ) closure. 2.4.1 The Ornstein-Zernike equation In 1914, Ornstein and Zernike [2] introduced the direct correlation function, c(r), v ia the relation, The function /),(ri ,r 2 ) represents the total pair correlation between particles 1 and 2 , and is known as the total pair correlation function. We can see that / i ( r i , r 2 ) does indeed contain the total pair correlation between two particles, since when g(ri,r-i) = 1, and particles 1 and 2 are uncorrelated, /*(r i , r-2) is zero. The physical meaning of E q . 2.28 becomes clear if it is expanded recursively to give Chapter 2. Theory 16 The total correlation function consists of direct and indirect parts. The direct correlation function represents the direct influence of particles 1 and 2 on each other. The integral terms in Eq . 2.30 represent the influence particle 1 on 2 as propagated by an increasing number of intermediate particles. The O-Z equation can be viewed as the definition of c ( r i , r 2 ) , which is exact for a fluid with pairwise additive potentials. For an isotropic homogeneous fluid, the O-Z equation has a simpler form which, when generalized for a mixture becomes for the present calculations of thermodynamic quantities. In order to evaluate the internal energy at infinite dilution according to E q . 2.18. we need to find the derivative with respect to solute density of the solvent-solvent radial distribution function. This derivative can be determined, if following Garisto et.al [15], the pair correlation functions in the O-Z equation are expanded in powers of solute density as, Chapter 2. Theory Substituting Eqs. 2.34-2.37 for h(r) and c{r) into Eq . 2.28, and collecting terms of equal power in solute density gives the O-Z equations for the pair correlation functions at infinite dilution, and their derivatives with respect to p^. The equations are Upon Fourier transforming Eqs. 2.38-2.40, the integral equations in r-space become the following linear equations in k-space [2], Chapter 2. Theory 18 anc <"/^ = P» : , [ 0] ocuv{k) + ^ j — , (2.46) 1 PI/CM 1 (Vi/^A,'J respectively. In the formulation of the above equations the fact that /i„M = h^u is used. Since Eqs. 2.44 - 2.46 tire derived in the grand canonical ensemble [2], they are valid under constant volume conditions. Experiments however, are likely to be performed under constant pressure. In order to approximate a system under constant pressure, Eqs. 2.44 - 2.46 are rederived with the added restriction that upon addition of solute, an amount of solvent is removed to keep the total density constant. That is, pu = p — p^ where p is the constant total reduced density. The O-Z equations for a system at constant total density are, 1 - pcii(k) and 6 " M ~ T^ÏÏW) T T - m + i - p c K ( t ) ' ( 2 ' 4 8 ) The O-Z equation for the correlation functions at infinite dilution are the same for the constant density and constant volume formulations. The equation determining the deriva-tives of the solvent-solvent correlation functions however, contains an extra term. There-fore, as expected, thermodynamic quantities for the two cases wil l be different. 2.4.2 Choice of Closure Since, the O-Z equations are defining identities and both the total and the direct pair correlation functions are unknown, closures relating li(r) to C ( Ï ' ) , and Shfr) to 6c(r) are required to determine these functions completely. Several such closures have been proposed. The Perçus-Yevick ( P - Y ) approximation, //.(/•) - c ( r ) = e x p [ - i « ( r ) ] ( / ( r ) - l (2.49) Chapter 2. Theory 19 was derived by Perçus and Yeviek in 1958 [7], using effective field techniques. The same equation was derived by Stell in 1963 [7], through graphical analysis of g(r) and c(r). The P - Y approximation is important because it has an analytical solution for hard sphere fluids. The mean spherical approximation ( M S A ) , g(r) = 0, r < d , c(r) = -pu(r), r > d , (2.50) is superior to the P - Y equation for square well fluids. Furthermore, it yields analytical solutions for some systems such as simple electrolyte and polar l iquid models [2]. For hard sphere fluids the P - Y and M S A approximations are equivalent. A l l distribution functions defined here can be expressed as many-dimensional inte-grals over the particle coordinates. Eq . 2.5 is an example of such an expression for the n-particle density. These integrals can be represented by graphs, and manipulated accord-ing to the topological properties of the diagrams [2]. Through such graphical analysis, diagrammatic expansions of c(r), and h(r) have been derived. From these expansions, the exact equation, hij(r) = djir) - Haij(r)] - flu^r) + BtJ(r) (2.51) resulted [2]. Together with the O-Z equation, Eq . 2.51 yields two equations in the three unknowns, c(r), h(r) and B(r). The last term represents a class of diagrams known as elementary or bridge diagrams. Through further graphical analysis a third relation between g{r) and B(r) can be derived. This however, is of no practical use since it introduces an infinite series of complex diagrams. Therefore, in practice B(r) is approx-imated. B y setting B(r) = 0, the hypernetted chain ( H N C ) approximation is obtained. A generalization of the H N C equation, the reference hypernetted chain ( R H N C ) equa-tion, approximates B(r) by the bridge functions of a reference system with a short range Chapter 2. Theory 20 potential. For the present calculations, the H N C approximation is used. In order to obtain a closure relation for Eqs. 2.38-2.40, the distribution functions in the H N C equation are expanded in powers of density according to Eqs. 2.34-2.37. Collecting terms of equal power in ptl gives the closure relations One of the advantages of choosing the H N C approximation is that an analytical expression can be obtained for the Helmholtz free energy at infinite dilution [18]. The formulation begins with the expression for the Helmholtz free energy for a mixture at finite concentration, Chapter 2. Theory 21 The derivation of Eq . 2.54 is somewhat lengthy, and is therefore given in the Appendix. To evaluate the excess Helmholtz free energy of solvation per solute particle, we first write the difference in free energy per solute particle of the mixture and the pure solvent, To evaluate the infinite dilution l imit , h.(r) and c(r) are again expanded in powers of solute density. In the l imit of p^ = 0, this gives Chapter 2. Theory 22 solute particle at infinite dilution is /j < Aex > = Pu (2.61) The significance of this expression comes from the fact that in its derivation no assump-tions about the form of the solute-solvent pair potential was required. Hence, unlike Eq . 2.26, which is applicable only to hard solutes, Eq . 2.61 is valid for a general potential. 2.4.3 Numerical Details The H N C integral equations are solved by a standard iterative technique known as the Picard method [2]. A n ini t ia l guess is made for the values of c;j(r). For our calculations, this ini t ial guess for the homogeneous fluid is the analytical P - Y solution of a hard sphere fluid at the appropriate density. Calculations for mixtures of increasing solute size have the direct correlation functions of the mixture with the nearest smaller solute as input. These correlation functions are fast Fourier transformed, and used to calculate fj^\k) from equations (2.38) and (2.39).The function i]ff(r) is recovered by calculating the inverse Fourier transform of iiff{k). The ?/{^(r) are substituted in the closure relation to obtain a new value for cfj. This procedure is repeated until the values of cfj converge. In order to force convergence, each guess is a combination of the results of the two previous iterations, computed according to A value of 0.90 is used for a. A similar iteration procedure is used for the solution of Eqs. 2.46 and 2.53. The input consists of the infinite dilution results and the derivatives 6cjj(r) for the solute with the nearest smaller diameter. During each iteration SJ9J points are calculated, with grid size output + ac\ J-i 'output ' (2.62) of 0.01. Chapter 3 Results and Discussion 3.1 Input Parameters The two systems studied were hard sphere and Lennard-Jones solvents, each with a single hard sphere solute particle. The parameters for which calculations were performed are given in Table 3.1 , where p* = pa* is the reduced density and T* = ksT/e is the reduced temperature. System P* H S / H S 0.80 H S / H S 0.60 H S / L J 0.80 1.35 H S / L J 0.60 1.35 H S / L J 0.50 1.35 H S / L J 0.80 1.10 Table 3.1: Input Parameters; H S / H S refers to a hard sphere solvent with a hard sphere solute; H S / L J refers to a Lennard-Jones solvent with a hard sphere solute. 23 Chapter 3. Results and Discussion 24 3.2 R a d i a l D i s t r i b u t i o n F u n c t i o n s The functions relevant to the calculation of thermodynamic quantities are (/[^(r), and <5<7„„(r). To illustrate the major trends observed, a representative sample of the calculated distribution functions are plotted. The variable ?• in the figures is in units of solvent diameter, d. In F ig . 3.1 the solute-solvent radial distribution functions for solutions consisting of a single hard sphere solute in a hard sphere solvent are plotted. The three curves represent calculated for three different solute sizes. The solid, dotted, and dashed curves correspond to solutes with diameters one, ten, and one-hundred times the diameter of the solvent respectively. The clotted and dashed curves have very similar shapes, il lustrating the fact that <y|,0J(r) becomes independent of the diameter of the macroparticle for sufficiently large solutes. The discontinuity in each distribution function occurs at solute-solvent contact, given by (c/„ + d^)/2\ however, for comparison purposes, the plots oi gl°](r) are shifted toward the origin along the r-axis by a distance equal to difference between solute-solvent and solvent-solvent contact. F ig . 3.2 shows the solvent-solvent radial distribution functions for solutions consisting of a single hard sphere solute in a L-.J l iquid. The meaning of the curves is analogous to those in Fig . 3.1. Again , it is clear that as the diameter of the macroparticle increases <7^(r) becomes independent of solute size. Once the diameter of the solute is fifteen times that of the solvent, the shape of ^{^(r) remains vir tually unchanged as the diameter of the macroparticle is increased. The same trend was observed by Perram and Whi te [8] in their calculations for hard sphere fluids, using the P - Y approximation. The fact that becomes independent of solute size indicates that a large particle is a good representation of a flat wal l . We emphasize that this is true only for separations which are small compared to the solute diameter. Chapter 3. Results and Discussion 25 Figure 3.1: Solute-solvent radial distribution functions for H S / H S solutions; p* = 0.8. The solid, dotted and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively. Chapter 3. Results and Discussion 26 Figure 3.2: Solute-solvent radial distribution functions for H S / L J solutions; p* = 0.6, T* = 1.35. The solid, dotted, and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively. Chapter 3. Results and Discussion 27 If the distance between the solute and the solvent particles is of the order of magnitude of the solute diameter, finite size effects become significant. For completeness, the solvent-solvent distribution function for the pure L - J solvent is given in F i g . 3.3. The function </f°J(r) for the pure hard sphere solvent is, of course, identical to ^[^(r) for the hard sphere mixture with equal solute and solvent diameters, given in F ig . 3.1 by the solid curve. In F i g . 3.4, derivatives of the solvent-solvent distribution functions are plotted for a hard sphere solvent containing hard sphere solutes of increasing size. The same deriva-tives for a L - J solvent are given in F ig . 3.5. Since 8guu(r) scales as the volume of the solute particle, the function is divided by the cube of the diameter. As for g]fy(r), the shape of 8gvl/{r)l d* becomes independent of the magnitude of solute diameter as the macroparticle is increased. In order to extract the surface term which is expected to scale as the surface area d2, the general form of the derivative is assumed to be 8g„„(r, d) = Sgi{r)dJ + 8g2{r)d2 . (3.63) That is, we assume that the dependence of 8gui/(r,d) on solute diameter d, can be sep-arated from that on particle separation r. In order to test the validity of E q . 3.63, 8g1/u(r)/d2 for a fixed r as a function of solute diameter is plotted in Figs. 3.6-3.9. In the range where the assumption is valid, these plots should be linear. It is seen from Figs. 3.6-3.9 that at short distances from the macroparticle, the graph becomes linear for solutes as small as two to three times the solvent. As expected, for larger r , linearity is observed only for larger solutes. Given the values of 8gul/(r,d) for two different solute diameters, two equations in two unknowns are obtained, which can be solved to yield 8g\{r) and Sg2(r), 8gw{r,di) = 8g\(r)d{ + 8g2{r)d\ 8gvvii\d2) = 8g\{r)d\ + 8g2{r)d\. (3.64) Figure 3.3: Solvent-solvent radial distribution function for H S / L J solution; p* = 0.6, T* = 1.35. Figure 3.4: Derivatives of solvent-solvent radial distribution functions for H S / H S so-lutions; p* = 0.8. The solid, dotted, and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively. Chapter 3. Results and Discussion 30 Figure 3.5: Derivatives of solvent-solvent radial distribution functions for H S / L J so-lutions; p* = 0.6, T* — 1.35. The solid, dotted, and dashed curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively. Chapter 3. Results and Discussion 31 Figure 3.6: 6gui/{r) at r = id for the H S / L J solution; p* = 0.6, T* = 1.35. Chapter 3. Results and Discussion Chapter 3. Results and Discussion Chapter 3. Results and Discussion Figure 3.9: 8gvv(r) at r = 307/ for the H S / L J solution; p* = 0.6, T* = 1.35. Chapter 3. Results and Discussion 35 The functions obtained by solving Eqs. 3.64 using three different sets of solute diam-eters are plotted in Figs. 3.10-3.11. Specifically, Eqs. 3.64 were solved for the following macroparticle sizes d\ = I40dsoivent d-2 = I20dsoivent dx = I40daoivent d-2 = I00dsoivent d\ = V20daoivent d2 = 90dsû[vent . Figs. 3.10 and 3.11 show 8g\(r) and Sg2(r) for the full range of ?• at which the functions were calculated. The three curves obtained by solving Eqs. 3.64 using three different sets of solute diameters are identical up to a solute-solvent separation of approximately 20 solvent diameters. For r greater than 20d, the curves are considerably different. Therefore, given the macroparticle diameters used to calculate Sgl(r) and Sg2(r), the separation of surface and volume terms according to E q . 3.63 is valid up to solute-solvent distances of approximately 20d. The range of validity can be increased by performing calculations for larger solutes, but for the present calculations 20ri is more than sufficient. In Figs. 3.12 and 3.13, Sgi(r) and Sg2(r) are plotted on an expanded scale to show the short-range oscillatory structure of the functions. Returning now to Eq . 2.40, the mathematical origin of the volume term can be iden-tified. The term involving the solute-solvent correlation functions necessarily introduces a volume dependence due to the reduction of configuration space. This dependence on the volume of the macroparticle is the result of inserting a particle of finite volume into solution. B y identifying the Sgl(r), the finite size effect is accounted for and the surface term 8g2(r), which contains the surface effect is extracted. Fina l ly in F ig 3.14 the derivatives of the solvent-solvent distribution functions calcu-lated with the added restriction that the total density is constant are given. Chapter 3. Results and Discussion 36 Figure 3.10: Sgl(r) for H S / L J solution; p' = 0.6, T* = 1.35. The solid, dotted, and dashed lines correspond to solutions of Eqs. 3.61 using different solute diameter pairs. Chapter 3. Results and Discussion 37 Figure 3.11: 8g2{r) for H S / L J solution; p" = 0.6, T* = 1.35. The solid, dotted, and dashed lines correspond to solutions of Eqs. 3.61 using different solute diameter pairs. Chapter 3. Results and Discussion Figure 3.12: 8g\[r) for H S / L J solution on an expanded scale; p* = 0.6,T* = 1. Chapter 3. Results and Discussion Figure 3.13: 6g2(r) for H S / L J solution on an expanded scale; p* = 0.6, T* = 1. Chapter 3. Results and Discussion 40 Figure 3.14: 8g„„(r) at constant total density for H S / L J solution; p = 0.6, T* = 1.35. The solid, dashed, and dotted curves correspond to solvent-solute diameter ratios of 1:1, 1:10, and 1:100, respectively. Chapter 3. Results and Discussion 41 As expected, the functions are qualitatively different from those calculated at con-stant volume (cf. F ig 3.2), but the trend that 8g„u(r)/di becomes independent of solute diameter for large solutes is again observed. Chapter 3. Results and Discussion 42 3.3 Thermodynamics The excess energy, entropy, and Helmholtz free energy of solvation per particle were calculated for each system with increasing solute, diameter. Figs. 3.15-3.17 depict these thermodynamic, quantities as a function of the square of the solute diameter. The choice of the function to be plotted is again dictated by the aim to separate the terms which are cubic in solute diameter from those which are quadratic. In other words, it is assumed that the thermodynamic functions are of the form, < Aex > + Asurd'2 < Uex N, > = uuol<P + Usur d2 < Sex > = 8vot(P + Ssurd2 . (3.65) The slope and intercept in Figs. 3.15-3.17 give the volume and surface terms explicitly. Physically, the terms scaling as the volume and surface of the solute are associated with finite size effects and surface excess quantities respectively. Mathematical ly it is clear that these terms should separate for the internal energy. In the. previous section it was shown that at short range 8 g „,,{!') can be separated into two terms; one which scales as the volume and the the other as the surface area of the macroparticle. For the systems studied here, 8gul/(r) completely determines the value of the internal energy, as the first term in the energy expression E q . 2.18, is zero for hard sphere solutes. Also, since the major contribution of the potential is at short range, the short range structure of 8gl/u(r) determines the form of the internal energy. If the solute did not interact with the solvent through a hard potential, </[^ ](r) would contribute to the energy of solvation according to Eq . 2.18. As shown earlier, <7^(r) becomes independent of the diameter of sufficiently large solutes. Chapter 3. Results and Discussion Figure 3.15: Excess Energy of Solvation as a Function of Solute Diameter. Figure 3.16: Excess Free Energy of Solvation as a Function of Solute Diameter. Chapter 3. Results and Discussion 45 Figure 3.17: Excess Entropy of Solvation as a Function of Solute Diameter. Chapter 3. Results and Discussion 46 The only effect of the solute is to shift the function along the r axis, so that it becomes non-zero for 7" greater than solute-solvent contact. Since g\°l(r) is multiplied by r 2 and integrated over ?•, the shift introduces a third order dependence on the position of contact, known generally as the excluded volume effect. A surface dependence would arise from the term in g\°\{r), since the number of solute-solvent interactions is proportional to the surface area of a macroparticle. Tables 3.2 and 3.3 contain the values of the two contributions to the internal energy of the L-.J solution obtained via two routes. The internal energy of a hard sphere solution is identically zero. The numbers in Table 3.2 come from numerically integrating Eq . 2.18 for a number of solutes and obtaining the slope and intercept by extrapolating the linear region of plots such as F ig . 3.15 for each system studied. In Table 3.3, the volume and surface terms are calculated by integrating the equations, Uvot = J47rr2dr{Ul/l/{r)6gl(r)} Usur = J4irr2dr{uu,(r)tg2(r)} . (3.66) The only true difference between the two methods of calculation is that when using Eqs. 3.66, graphical analysis is not needed. The discrepancy in the numbers obtained for the contributions to the excess internal energy by the two methods is a result of numerical inaccuracies. If the graph of < Uex > / A r / t is extrapolated from sufficiently high values of d, the slope and intercept approach the values obtained by integrating Eqs. 3.66. The volume term of the internal energy is negative for all systems studied, indicating that when configuration space is reduced, particles rearrange in a way that lowers energy. The magnitude of this component increases from a density of 0.50 to 0.60, but decreases substantially from a density of 0.60 to 0.80. A t high density, the repulsive part of the potential is dominant, and the fact that the volume term becomes less negative can be attributed largely to the increase in effective density when the solute is added. Chapter 3. Results and Discussion 47 Table 3.2: Excess Energy of Solvation from extrapolation of < Ucx > /N^d2 vs. d plots. System />* T* p < Uex > /JVM P < Uex > /N„ Volume Term Surface Term H S / H S 0.60 0.0000 0.0000 H S / H S 0.80 0.0000 0.0000 H S / L J 0.50 1.35 -0.5446 -1.3334 H S / L J 0.60 1.35 -0.7201 -0.5979 H S / L J 0.80 1.35 -0.3563 2.2374 H S / L J 0.80 1.10 -0.7236 2.2700 Table 3.3: Excess Energy of Solvation from integration of 6 g\(r)and 8g2(v). System P* T* p < Uex > /Nfl p<Uex> Volume Term Surface Term H S / H S 0.60 0.0000 0.0000 H S / H S 0.80 0.0000 0.0000 H S / L J 0.50 1.35 -0.5438 -1.3998 H S / L J 0.60 1.35 -0.7189 -0.6938 H S / L J 0.80 1.35 -0.3554 2.1816 H S / L J 0.80 1.10 -0.7215 2.0572 Chapter S. Results and Discussion 48 This conclusion is supported by the values obtained for the internal energy at constant total density, given in Tables 3.4-3.5. The magnitude of the volume term does not decrease as much at p = 0.80 in the latter case, because the effective density is not increased upon introduction of the solute. The slight decrease in the magnitude of the volume component of the internal energy of high density liquids calculated at constant total density implies that at high density the rearrangement of solvent particles is less effective in lowering the energy. The surface excess energy increases with solvent density, becoming positive at p = 0.80. The same trend is observed in the numbers calculated at constant total density. The surface excess energy at high density is greater in the former case because the actual solvent density is higher. When the bulk density is low, the increase of the number of particles near the surface results in a negative contribution to the internal energy. For a high density l iquid, repulsive forces become dominant as the number of solvent particles near the surface is increased, resulting in a positive contribution to the excess internal energy of solvation The Helmholtz free energy is also calculated using two different equations. Both methods involve obtaining the free energy for a number of solute particles with increasing diameters, and extrapolating the linear region of plots such as F ig . 3.16 to find the slope and intercept. The values in Table 3.6 were calculated by integrating E q . 2.26, while those in Table 3.7 come from using Eq . 2.61. The latter has the advantage that it is not restricted to hard solutes. Both the volume and the surface contributions to the Helmholtz free energy are found to be positive. For a given density, the contribution of the volume term is greater in a hard sphere than in an L - J solvent. The component which scales as the surface area is smaller for a hard sphere fluid at low density and comparable for the two solvents at high density. Chapter 3. Results and Discussion 49 Table 3.4: Excess Energy of Solvation at constant total density from extrapolation of < Uex > /N^d2 vs. d plots. System P* T* p < Utx > IN» p < Uex > /N^ Volume Term Surface Term H S / L J 0.50 1.35 -0.5437 -1.4158 H S / L J 0.60 1.35 -0.7261 -0.6005 H S / L J 0.80 1.35 -0.6292 0.7689 Table 3.5: Excess Energy of Solvation at constant total density, from integration of 6gl(r)and 8g2(v). System P* T* p < Uex > 1 p < Uex > /N„ Volume Term Surface Term H S / L J 0.50 1.35 -0.5442 -1.3626 H S / L J 0.60 1.35 -0.7191 -0.6546 H S / L J 0.80 1.35 -0.6284 0.7104 Chapter 3. Results and Discussion 50 Table .3.6: Excess Free Energy of Solvation, from Eq . 2.26. System />* T* /j < Aex > /.\„ P < Aex > /7VM Volume Term Surface Term H S / H S 0.60 1.5783 2.2080 H S / H S 0.80 4.0954 5.3757 H S / L J 0.50 1.35 0.2253 1.1695 H S / L J 0.60 1.35 0.4845 2.8352 H S / L J 0.80 1.35 2.06927 7.2559 H S / L J 0.80 1.10 1.89931 4.2404 Table • 3.7: Excess Fr ee Energy of Solvation, from Eq . 2.61. System P* T* < Aex > (3 < Atx > /N„ Volume Term Surface Term H S / H S 0.60 1.6483 2.9238 H S / H S 0.80 4.2935 5.2179 H S / L J 0.50 1.35 0.2342 1.4730 H S / L J 0.60 1.35 0.5071 2.6999 H S / L J 0.80 1.35 2.1667 6.4010 H S / L J 0.80 1.10 1.9824 6.6915 Chapter 3. Results and Discussion 51 The values obtained for the components of the Helmholtz free energy by the two methods are considerably different. The two routes to the free energy are not thermody-namically different, and within the H N C approximation are mathematically equivalent as well. We believe that the source of the discrepancy in the numbers obtained by the two methods is numerical. The difference in the Helmholtz free energy for each solute between the two techniques is approximately 4-5 %. The discrepancy is magnified when the numbers are used in graphical analysis to determine the slope and intercept. In Table 3.8, calculated values for the excess entropy or solvation are given. These are obtained by using the thermodynamic relation, given in E q . 2.27, and subsequently finding the slope and intercept of graphs such as F ig . 3.17. Both contributions to the entropy are found to be negative, and becoming more negative with increasing solvent density. The fact that the volume term is negative indicates that when configuration space is reduced with or without a corresponding increase in bulk density, the entropy of the system decreases. The negative contribution of the surface excess entropy is interpreted as the increased spatial ordering of the liquid near the surface. Chapter 3. Results and Discussion 52 Table 3.8: Excess Entropy of Solvation. System ^* T* fi < .S'"' > /kBN„ Volume Term 0 < Se* > /kBN„ Surface Term H S / H S 0.60 -1.5783 -2.2080 H S / H S 0.80 -4.0954 -5.3757 H S / L J 0.50 1.35 -0.7699 -2.5029 H S / L J 0.60 1.35 -1.2046 -2.8912 H S / L J 0.80 1.35 -2.4256 -5.0185 H S / L J 0.80 1.10 -2.70438 -4.5741 Chapter 4 Conclusions The Ornstein-Zernike (O-Z) equation for a two component mixture consisting of parti-cles wi th spherically symmetric potentials was solved with the hypernetted-chain ( H N C ) approximation to yield solute-solvent and solvent-solvent pair correlation functions at infinite dilution. The O-Z equation for the derivative of the solvent-solvent radial distri-bution function êgul/(r), is obtained by expressing the pair correlation functions as power series in solvent density. The function 6 </„„(»•) contains information about the rearrange-ment of solvent particles upon introduction of a solute particle. Solutions to both sets of coupled equations were obtained via an iterative procedure. Results were presented for solutes with diameters ranging in size from 1 to 150 times that of the solvent. The mixtures studied consist of hard solutes of varying sizes in hard sphere and Lennard-Jones (L-J) solvents, at constant volume and a number of temperatures and solvent densities. A set of results was also presented for mixtures at constant total density. Each set of conditions defines a state well away from phase transitions. It is well known [9] that for large solutes the shape of the solute-solvent radial distri-bution function (r), becomes independent of the macroparticle diameter. Therefore, the large solute l imit of a mixture at infinite dilution can be used to represent a bulk fluid in contact with a flat surface. In this thesis, the dependence of Sg„„(r) on solute diameter was investigated. The function Sgul/(r) was shown to be a sum of two dominant terms; one cubic 53 Chapter 4. Conclusions 54 and one quadratic in solute diameter. The two terms were separated, both graphically and algebraically. The cubic term was associated with the change in solvent-solvent correlations due to the exclusion from configuration space of the volume of the solute particle. It was found to be a function which has a short-range oscillatory structure and which decays to zero only at long-range. The term which scales as the surface area is interpreted as the function describing solvent rearrangement upon insertion of a flat surface into the pure solvent. The latter function also exhibits short-range oscillations, and contains a negative contribution to the total change in solvent correlations. The radial distribution functions and their derivatives were used to evaluate excess thermodynamic quantities of solvation. The excess internal energy, Helmholtz free energy, and entropy of solvation were shown to lie separable into cubic and quadratic terms in solute diameter. The excess internal energy of solvation at infinite dilution can be obtained from in-tegrals over the solvent-solute radial distribution functions and the derivatives of the solvent-solvent radial distribution functions. In the present calculations, only hard so-lutes were considered. Hence, the excess internal energy of solvation is determined entirely by the change in solvent structure due to the presence of the solute, as the contribution of solvent-solute correlations is zero. At low densities, both the volume and the surface terms of the internal energy were found to be negative. At higher density, the volume term becomes more positive. This effect is less pronounced when the total density is kept constant, and was therefore attributed in part to the increase in effective density upon insertion of a macroparticle into a l iquid at constant volume. As the density is increased, the surface excess energy also becomes more positive, and is greater than zero at high density. This trend was interpreted as a result of increased solvent density near the surface. The excess Helmholtz free energy of solvation was calculated via two routes. One Chapter 4. Conclusions 55 technique is exact for an exact distribution function, but is restricted to hard solutes. The second method, which relies on the H N C approximation, is applicable to general solutes. The excess free energy was found to be positive for all solutions considered, indicating that inserting a macrosphere into a hard sphere or a L - J fluid increases the chemical potential of the l iquid. From < Uex > /7VM and < Aex > jN^, the excess entropy of solvation was found by evaluating the thermodynamic relation connecting internal energy, Helmholtz free energy, and entropy. Both the surface and volume terms of the calculated excess entropies are negative, indicating that both finite size and surface el feet s involve increased spatial ordering of the l iquid. The present approach can be extended to charged systems and liquids with orienta-tionally dependent potentials, such as water. In particular, contribution of orientational ordering to the entropy of liquids or liquid crystals near a flat surface can be studied. In its current form, the procedure can be used to obtain information about liquids interact-ing with macroparticles and surfaces through hard and soft potentials. Bibliography [1] J .L .Lebowi tz , E .M.Waisman , The Liquid State of Matter: Fluids Simple and Com-plex E . W . M o n t r o l l and J .L.Lebowitz , éd., North Holland Publishing, (1982) [2] J -P. Hansen, I. R. McDona ld , Theory of Simple Liquids, 2nd. ed., Academic Press (1986) [3] D . Henderson, M . Plischke, J . Phys. Chem. 92, 7177 (1988) [4] R. Zwanzig, J . Chem. Phys. 22, 1420 (1954) [5] L . E . Reichl , A Modern Course in Statistical Physics, University of Texas Press,(1980) [6] T . L . H i l l , Statistical Mechanics, New York: Dover Publications Inc., (1987) [7] D. A . McQuarr ie , Statistical Mechanics, Harper and Row, (1976) [8] J . W . Perram, L . R. Whi te , Faraday Discussions o the Chemical Society, 59, 29 (1975) [9] D . Henderson, J . A . Barker, F . F . Abraham, M o l . Phys. 31, 1291 (1976) [10] L . S. Smith , L . L . Lee, J . Chem. Phys. 71, 4085,(1979) [11] G . M . Torrie, P. G . Kusal ik , G . N . Patey, J . Chem. Phys. 88, 7826, (1988) [12] D . R. Berard, G . N . Patey, J . Chem. Phys. 95, 5281 (1991) [13] R. Kjellander, S. Sarman, Molec. Phys. 70, 215 (1990) [14] D . R. Berard, G . N . Patey, (to be published) [15] F . Garisto, P. G . Kusal ik , G . N . Patey, J . Chem. Phys. 79, 6294 (1983) [16] S. Singer, D. Chandler, M o l . Phys. 55, 621 [17] L . Frisch, J . L Lebowitz, H . Reins. J . Chem. Phys. 31, 369 (1959) [18] H - A . Y u , M . Karplus , J . Chem. Phys. 89, 2366, (1988) [19] F .Lado, M o l . Phys. 47,283 (1982) 56 Appendix Helmholtz Free Energy for a Mixture In this section, the Helmholtz free energy for a two component mixture at finite con-centration is derived. The details closely follow the analysis of fluids with non-spherical potentials by Lado [19], and the work of Singer ct.al. [16] with the extended R I S M integral equations. To evaluate the Helmholtz free energy arising from the interaction of particles, the coupling parameter integration method is again followed. This time, however, the param-eter, £, turns on and off all interparticle interactions. The total potential of the system is defined as a function of £ by LIN{ri,riV; 0 = £ "iA'u) • (A.l) i<j The subscript is included in the pair potential to emphasize that in a mixture the type of interaction depends upon the nature of particles /' and j . To continue, the Helmholtz free energy as defined in the canonical ensemble is evalu-ated. The configuration integral now has the form, ZN(Ç) = Jvexp[-!UrH(ru...,rN;Z))drN . (A.2) As Eqs. 2.21 and 2.22 indicate, the derivative of the Helmholtz free energy with respect to the coupling parameter has to be found. W i t h the above configuration integral, this becomes, dp < A{Ç) > 1 dZN(Q --zW)~lîT- (A-3) 57 Evaluating E q . A.3 explici t ly yields, which is the canonical average of the potential energy. The above equation, therefore, can be written in terms of the radial distribuiou function as, The notation of E q . A.5 is specific to an isotropic, homogeneous, two component mixture. The integrand can now be manipulated to a form which is an exact differential with respect to £. To achieve this, it is convenient to introduce two new functions, 59 Upon rearrangement of E q . A.10 the integrand of Eq . A.8 can be expressed as a sum of an exact differential, and another term which needs to be manipulated further. The expression in its useful form is df d / 1 7 2 . dc The total contribution of the second term to the Helmholtz free energy per particle of a mixture is, 2^Ex.-.Yi \\h%Vàr. (A.12) This sum can be expressed as an exact differential of £ if O-Z equation for a mixture in fe-space is written in a matr ix form, ( A . l l ) h = c + h • p • c , (A.13) where (A.14) and c = (A.15) (A.16) Pu 0 0 P» M u l t i p l y i n g both sides of Eq . A.13 by the matrix p, and adding the identity matrix gives, after some rearrangement, the relation, [I + ph]-' = [ I - p c ] . (A.17) we now have the integrand of A.8 as an exact deferential with respect to the coupling parameter. Therefore, the integration with respect to £ is t r iv ia l . The final form of the Helmholtz free energy of a two component mixture is
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Statistical mechanics of solvation of macroparticles
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Statistical mechanics of solvation of macroparticles Vasarhelyi , Reka Z. 1992
pdf
Page Metadata
Item Metadata
Title | Statistical mechanics of solvation of macroparticles |
Creator |
Vasarhelyi , Reka Z. |
Date Issued | 1992 |
Description | Integral equation theories are employed to study the solvation of large particles. The solutions studied consist of spherically symmetric solvent particles which interact through hard sphere and Lennard-Jones potentials, and hard sphere solutes. In particular, the Ornstein-Zernike (O-Z) equation is solved with the hypernetted-chain (HNC) closure, to obtain the pair correlation functions of mixtures at infinite dilution. The pair correlation functions in the O-Z equation and the H N C closure are expressed as power series in solute density to yield a pair of coupled equations which determine the derivatives with respect to solute density of the solvent-solvent pair correlation functions. The latter describe the perturbation of the solvent upon addition of a single solute particle. The derivatives are analysed to yield components that scale as the volume and surface area of the macroparticle, and which are then identified as changes in solvent structure due to the presence of a finite size particle and a flat surface respectively. From the pair correlation functions and their derivatives the excess internal energy, Helmholtz free energy, and entropy of solvation are calculated. The excess quantities are also separated into contributions from finite size and surface effects. Both components of the excess internal energy are negative at low densities, and become less negative for high density liquids. The magnitude and sign of the two contributions to the energy depend on physical conditions such as temperature and pressure. The excess entropy of solvation is negative for all systems studied, indicating that introduction of a macroparticle or a flat surface increases the spatial ordering of a bulk liquid. |
Extent | 1458389 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-12-18 |
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.0061739 |
URI | http://hdl.handle.net/2429/3131 |
Degree |
Master of Science - MSc |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1992-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1992_fall_vasarhelyi_reka_z.pdf [ 1.39MB ]
- Metadata
- JSON: 831-1.0061739.json
- JSON-LD: 831-1.0061739-ld.json
- RDF/XML (Pretty): 831-1.0061739-rdf.xml
- RDF/JSON: 831-1.0061739-rdf.json
- Turtle: 831-1.0061739-turtle.txt
- N-Triples: 831-1.0061739-rdf-ntriples.txt
- Original Record: 831-1.0061739-source.json
- Full Text
- 831-1.0061739-fulltext.txt
- Citation
- 831-1.0061739.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-0061739/manifest