Molecular dynamics simulations of restricted primitive model 1:1 electrolytes

Chemical Fhysics 69 (1982) 155-163 Korth-Holland PubIishinr Company

MOLECULAR DYNAk4PCSSEWJLAT1QNS OF RESTRICTED PRIM~VE

MQDEL

1: 1 ELECTROLYTES

D.M. HEYES

Received 19 November 1981; in tinai form 22 March 1982

A new mokcula~ dynamics method is described which enables the simultaneous simulation of hard sphere collisions in the presence of “soft” forces. Calculations show that the Debye-Hiickel theory, DHX, describes well the equation of state of clzrged hard spheres up to an equivalent electrolyte concentration of 1 hf. This results from a cancellation of errors, rather than the validity of the underlying assumptions. Better agreement for higher concentrations up to 2 M is obtained with previous Mor.teCarlo simulationsand higher order theories such as MSA and HNC. A two-. three- and four-particleresolution of the specificheat revealsthat the four-ptiicle tern is mainly responsible for the poor statistin of this quantity. Below3 M s&diffusion is well describedby Enskog’stheory. However,there is an approximately 15% detiatiop. from the NernstEinstein equation.

1. Introduction The computer simulation techniques of molecular dynamics, MD, and Monte Carlo, MC, have been used to investigate the thermodynamic and dynamical properties of ionic solutions [1,2] . There are technical problems associated with the MD calculations that do not arise with the simulation of the pure solute (i.e. a molten salt) [3]. As the reorientational motions of the water molecules are much faster than simple translation of tha species, tie solute and solvent evolve on essentially two different time scales. As a result, the satisfactory exploration of phase space is computationally expensive. Heinzinger and co-workers [!] have simulated aqueous alkali halides by MD using a five-interaction centre model for the water molecule. The combination of light and heavy mass centres necessitated a time step of duration 2 X IO-l6 s, which is roughly an order of magnitude smaller than iS normally used in MD [3] _ The model solutions were evoived for only several picoseconds and consequently a limited region of phase space was explored. As the atomic masses do not enter expiicitly in the procedure for MC equilibration better statistics can be achieved [2]. The limitations on MD simulation are being overcome, however, with the use of array processors 141, which can achieve an order of magnitude in-

0301-0104/82/0000-0000/!$02.75

0 1982 North-Holland

crease in simulation time and the development of simpler effective three interaction centre models for water [S] , in which a time step equivalent to I X IO-l5 s can be used. A less rigorous model for ionic solutions has been implemented which permits a more acceptable MD time step of 4 X IO-l5 s by eliminating the atomic components of each solvent molecule [S] _The ionic species interacted together with soft short range potentials and in addition each solute ion was subjected to a stochastic force in order to represent the random buffetings of the solvent molecules. Although it is doubtful whether it was valid to use brownian dynamics for solute and solvent particles having similar masses. In this work the process of simplifying the solvent LTcontinued further 2s the aim is to assess only those effects produced by the solute-solute Coulomb interactions on the thermodynamic and dynamical behaviour of the ions in solution. The solvent is represented by a cont*buum of uniform dielectric constant and zero viscosity. It is obvvous that this gross over simplification invalidates any attempt to calculate absolute thermodynamic and dynamical quantities of the reai system. However, as the solvent-solvent and ion-solvent interactions are aIso present at Siiite dilution then, to a first approximation, they should not contribute appre-

ciably to those changes on the solution’s thermodynamits and dynamics which are introduced only by the solute-solute Coulomb interactions. Debye and Hiickel proposed a model for dilute eIeCuol_yte solutions which accounts satisfactorily for the deviations from ideal behavioru for concentrations below ~1 molar, mainly in terms of the electrostatic solute-soInte interactions [7]. They represented the complex system of ions and polar solvent molecules by a distribution of hard spheres with point charges at their centres which move in a continuum with a uniform dielectric constant. In order to express analytically the statistical mechanics of this restricted primitive model, RPM, system certain mathematical approximations were made [i] _ The effectiveness of these assumptions for this system has been tested by Monte Carlo calculations [8,9] which numerically determined the equation of state of the particles. This is continued here with the complementary method of molecuIar dynamics, which also enables an assessment of the effect of the Coulomb forces on the dynamics of the RPM model in the dilute regime to be made. Here an infinitely hard repulsive interaction is favoured as it ciearly distinguishes between the success of the mathemztical approximations in the Debye-Hiickel theory (and its related succe~ors [9]) and the uncertainties of the ion-solvent interaction. Any softer repulsive interaction would be arbitrary and would obscure the effect of the electrostatic contributions. In addition, the length of the time step is determined primarily by the form of the repulsive part of the potential, as the degree of energy conservation is used as the criterion. Hard spheres do not contribute to the potential energy. Thus a much larger time step of 25 X IO-l4 s becomes acceptable ar,d enables a greater region of phase space to be covered. Any more realissic solvent which accounted for its granularity would limit the number of model ions in the MD cell and would require a parallel dynamical evolution. Also the time step would necessarily be =10-i’, s. Such a solvent representation would produce poorer statistics for the solute behaviour and would be a costly indulgence. This is particularly important for ion self-cliffusion coefficients which decrease in real solutions by only severai percent over the concentration ranges studied below [lo,1 I] _ In order to perform these calculations, two methodologically different forms of MD have been combined in the same MD program. This allows the perfectly elastic

collisions ofhard spheres to take place While also eva& uating the forces and velocity changes due to a smoothly vargig continuous interaction potential. The details of the model and its computational implementation are discussed in the next section.

2. The model A cubic MD cell of volume Vcontained N hard spheres of diameter o; N/7, with charge +q and N/2. with charge -q. The RPM potential energy between particles i andj separated by a distance r is, q&(r) = Z,ziEor-1 , r>l7; 0)

=m >

r
(2)

where the well depth E = q’/4moero

in which er is the dielectric constant of the medium relative to that of free space, co_ For positive and negative charges, then, ZX_= +I and -1, respectively. Ta!-5ng parameters of E,, = 8.8542 X IO-l2 3-l C2 m-l, cr = 78356, a temperature of 298.16 K, u = 4.25 _%[9] and mass, 1~ = 100 amu (corresponding to four water molecules in the first hydration shell of a light ion) then e/kB = 501.78 # and the reduced unit of time, a(~+)~/~ = 2.081 ps. The Calculations were carried out at the reduced temperature, T= 0.59 (*O.Q3)ejkg and at number densities, p, rang-

ing from 8.474 X lOA (or 0.00911 M) to 0.2774Na3/V {or 3 M). The Coulomb interactions were calculated using the Evjen method 181. In contrast to the simulation of moiten salts, lattice summation techniques are not essential here because of the low density and high reduced temperature of this system [9,12] _ The continuously smooth component of the potential was used to increment the positions of the particles at successive times separated by a constant AZ, using the Verlet algorithm [13],

tixe interval,

Sri(r) = Fri(r - Ar) + r;;(r) At”@,

(3)

r&f f At) = r&f) f 6rj(r),

(4)

where Sri(t) is the position increment that is used to move particIe i from its position at time I, r#), to that at time, r f Ai, ri(t i At). The force on i at time r is defined to be I;i(t) and mi is the mass of particle i. Note that the position increment takes the place of a velocity in the Verlet scheme. The

velocity of particle i, Vi(t), is only of secondary importance and is evaluated at At interval5 using the form&

edits = Sri(t - At) Af’

+ ~~‘i(t) Atlm,.

The above sequence of operations

(5)

is iohowed

pro-

vided that there are no hard sphere collisions in the tie window from t to t + AL This is determined [ 141, by solving the quadratic equation below for each particle pair i and j, ‘rv i- t’rii - Fr,$Atl’

= 02,

or

where rii = r#

- ‘i(f),

(8)

and Srij = “r&f) - &$f)_

(9)

If t’ > At for alI particle pairs then eqs. (3)-(5) can be applied to eIl particles as the next predicted collision occur5 after the forces are next going to be evaluated at time t + At. However, if t’ = a,A.t, where 0 < CLI< 1 for at least one (ii) combination then a cycle of smaller time steps is commenced until there are going to be no

more collisions before time f t At. Between times c and t + At the particles undergo elastic collisions without velocity modification by the continuous forces. Thus if i is a particle which does not undergo a collision between t and t + At, but there are rzc collisions involving other particles within this time interval, then the following sequence of intermediate updates on its position are made, w (I I)

+ [ 1 - (aI + 9 +, __.mQ] Sri(f)_

(12)

However, if particle i collides at the end of small time-

step i, 6

k=l

ol,At ) --r. r( t f ‘2

k=?

a,Ar)

+9

Gr&t), (13)

then, Grj(t) --f -&r&f) for the next change in position of i at time, t + “YZJiolkAf_ The value of this position

increment (equivalent to a time &) changes sign at each collision with another particle,j. Computationally this series of collisions is quickly generated as follows. At each time step tables of next collision partners for each i and the times to these collisions are created. As

each collisionj,

takes place then 9Af is subtracted

from the times to the next collision for each particle. Only for each colliding pair is it necessary on a collision to search over all other particles to fIId their new collision partners and times. Then a single trip through aNdimensional loop gives the next collision partners. Note that two passages through a (two particle) double loop arc required each time step; first, to calculate the position increments [from the Fi(t)] and then immediately afterwards to set up the table of next collision partners for each particle using the position increments derived from the Fi(i). This scheme is perhaps less rigorous than a similar one proposed recently in which the new Verlet trajectory is recalculated on impact as well [ 151. The accuracy of the above algorithm for moving the particles is determined by the Verlet component for treating the soft forces, as the impulsive dynamics are exact to within the accuracy of the machine_ The error of thi5 operation is proportional to Afi. The majority of particles in these dilute systems proceed for a time step without colliding with another particle. It therefore Seem5 an umreces5ary refinement to update the Verlet trajectories more frequently than at times separated by the constant time interval, At. It should be noted that the Verlet algorithm gives poor energy conservation for single collisions between isolated particles. The success of the Verlet scheme here depends on the long range of the Coulomb interaction. Thus true single collisions do not take place even at these extreme dilutions. This would not apply for a van der Waals systern at these densities, for exarnrple. The computations were continued for approxzately 20000,500O and 2000 time steps for N = 64,216 and 512, respectively. The time ste duration, Ar, ranged from 0.025 to 0.0045 o@r/e)I P2 for the two densities quoted above. The average number of hard sphere collisions per time step was 0.4,l. and 2. for N = 64,2 16 and 512, respectively.

3. Resdts and discussion III keeping with a comparable ,ciC work [S] the MB simulation results are compared with thepredictions of a form cf the Debye-HCckel theory known as the DHX model [S] _Higher order theories 1161 such as the MSA, HN, mode coupling j17] and EXP [I8] approaches are also examined. Although these higher order theories are frequently better representations of the FLPblsystem, they are not easy to evaluate. In contrast, the DHX approach is analytically solvable and contains some important concepts which make it an approximation worth discussing in detail. me DHX model is not entirely self-consistent because in order tc solve the Poisson-Boltunarm equation for the potential of mean force it hearises the oair distribution functions, which means making them proportional to the potentiaI of mean force. However, the derived thermodynamic properties use a pair distribution function which is proportional to the exponential of the potential of mean force as given by the B~Itzmann equation in its original unlinearised form. Although there are snore theoretically justiRable models, [16], it is essential to determine the minimum requirements in a theory which adequately reproduces the simulation data. The DHX interionic pair radial distribution function,$HX(r), between species i and j is obtained from theDHX potential of mean force, qHX(r), as follows, p(r)

= exp(OgHx/k,T)*

@iH”(r) = zizjeo exp[li(a - r)] /(I + Ka)r,

(14) (15)

where K is the DHX reciprocal length, K =

(4ii&n”BT)‘J2.

(16)

fIere. (KjP)

= 0.33,903 (&QL’*,

(171

in which c is the concentration in mol dm-3 or M. The pair distribution functions of the 0.1 M and 2 K MD states are compared in fig_ 1 with those predieted by the DHX model. Even at a concentration of 0.1 hi there is more ion ptig than is predicted using DI-IX theory. At concentrations of O.:, 0.5,Z.O and 3.0 hl the MD (DHX) unlike ion pair distribution functions at contact,&(u), are 3.3 (3.2), 2.6(2.3), 2.2(1.S) and 2.3jl.6). The corresponding l&e ion values, g,Jo),

a,(r)

3

2

Fig. 1. The unlike ion,g,$), and like ion, g&), pai radial disrriiution functions for (a) c = 0.1 M: 0 MD, W= 64). A MC c!=zoo),~--

DHX,-.-.-CMSA[19],--hhlSA[19].

and @) c = 2.0 M: DMD CV= 64), A MC f.N= ZOO),--- DHX, esp [ 18]_ Note that c = 10.81495 p_

are 0.3 (0.31),0.46 (O-43), 0.66 (0.57) and 0.85 (0.61). The linearisation of th.e pair distribution functions in the Poisson-Bohzmann equation results in an underestimation of the attraction between unlike species at short range (r 5 2~) and consequently low values for g”(u) in this region result. The agreement with previous Monte Carlo calculations [S] is excel!ent even at 2 M, at which concentration the DHX approximation is clearly inadequate. Higher order theories, such as the EXP approach [la], which incorporate the finite size of all model ions (not just the origin particle as in D-FIX)also reproduce the simulation pair distribution functions more faithfully. However, Ehe MSA pair distriiution functions [19] have inadequacies at low density 1161. The configurational energy ratio, U~NZCBT,where, N

pL&

2i=l

i

(18)

and N ui = c I$.., i+j V

(1%

D.&f.Heyes/Simu&tionsof primitiveelectroljtes

also reflects this tendency. Fig. 2 shows that the MD C$VkBT increase more rapidly with density than the DHX value,

even allowing for the srnallhr-dependence in this region. The configurational energy, U, depends on the integration of g&) -g&r) over r. There is a tendency for both of ‘he hlD unlike and like species pair distribution functions to be weighted more to small r tharrthose of the DDX model. Consequently there is a fortuitously good agreement between the U,&R,T for c < 1 M produced by the two methods even though their pair distribution functions are qualitatively dissimilar from very low concentrations. An experimental U/Nk,T is obtained by dividing the integral heat of dilution [20] by -O.36’i9 [21], in order to correct for the temperature dependence of er for water. The experimental iJ/NkBT deviates at much lower concentrations (c = 0.1 M) from DH theory. Thus the increasing disagreement between DID: theory and experiment for concentrations greater than =O.l hi is due mainly to the simple representation of the solvent rather than the mathematical approximations involved

159

in treating the Coulomb interactions. This gives rise to a poor effective solute-solute pair potential. Ihe MSA and HNC equations of state are solutions of rhe (Fourier transformed) Omstein-Zemicke equation when subjected to specific closure relationships. The MSA and HNC approaches predict essentially the same U/McBT as the computer simulation values, as fig. 2 shows. The MD hard sphere component of the “instantaneous” pressure, P, was obtained from averaging during the time step the nc collision contributions to the virial which are evaluated on contact of the spheres [13],

The results of eq. (21) agree within statistical error with those calculated from the pair distribution functions at r=u

[S],

PV-/h’k,T = 1 + U/3h’kBT + $r~po~mu

+g&r)] . (22)

The DHX osmotic coefficients, were obtained using eq. (22) and the relevant pHX,gU(u) and gs(u) which in “real” units at T = 0.5964 efkg are: flHX/nrk,T=

-($-j4)

0.6 i

q&I = expC%

(24)

g&0 = exp(-A),

(25)

where A = l-68251/(1 f 13984 c”~)_

Fig. 2. The reduced configurational energies,-U/ARgT, as a

function of concentration, c, for (a) hfD: c, Aand o appIy to and 512, respectively.The vertical ties denote standard error estimates.(b) MC: 8 and A apply to N = 64 and 200, respectively:(c) DHX: -(d)HNG 1171 -.-.(e) h&W.,,[17] - - (f) experimental data on NaCi in water

N = 64,216

at 25°C [20,21] ---.

(26)

Fig. 3 compares the hfD, MC, DHX, HNC and MSA osmotic pressure coefficientsPV/!bT [8]. Again MD and MC values agree well. Also the HNC approximation obtained from the virial is far superior to the DHX, and in this case iLfSA,approaches. As fig. 1 would indicate the MSA pair distribution functions at contact are in poor agreement with simulation values. Consequently the virial route to the osmotic pressure coefficients for the MSA is not successful.

160

DX. fleyes/SSmulntions of priktitv electrolytes

PV EL$

much longer than predicted from the ranges of interaction; Izrgely due to spurious effects resulting from the periodic boundary conditions which enforce structural correlation between particles separated by = the MD box side length. As a first step to correcting for this artefact of the computing model, the formalism for resolving such fluctuations into their two-, three- and four-particle components is given below. This is applied to the evaluation of the confgurational energy component of the specific heat per particle, C,, which is obtained from *he fluctuations, b, in U for this microcanonical ensemble [25] _A decomposition of j3 into self,OS and cross, PC, terms w2s made as follows, C” = 1.5Ap/‘(l-&/kB),

F&. 3. The reduced osmotic pressures,PV/M~T. The key is as for fii_ 2, except that the hard sphere component of the pressure was obtained from the virial on contact [I?] for o, A and c whereas * corresponds to the MD pressuresobtained usinggo andgz(o) ir. eq. (22) for the smallest N at each concentration. The experimental data are qprapriate to aqueous NaCI at 2S=C and are taken from ref. [22] _.4lso, the mode coupI@ approximation [ 171 gives - - -.

There is a greater tendency ?o ion-pair with increasing density. Consequently PV/NkBT manifests a minimum at c = 0.1 M. At low concentration the coulombic ener,T is the main interaction contribution whereas the hard sphere collisions dominate at higher density. This phenomenon has a parallel in the Joule-Thompson effect in dilute gases [IO] . These is increasing interest in studying by computer simulation the time fluctuations of first order thermodynamic quantities as a method of obtaining second order thermodynamic properties (81 and transport coefficients [S] . Similarly structure and stricture related ensemble averages provide insights into the characteristics of phase changes [ 141 and the interpretation of scattering experiments 1231. However, it is now widely accepted that these fluauations have a greater X-dependence than the fust order properties from which they are derived. Progress has been made at attempting to understand this at a microscopic level by consider’;ng the specific example of fluctuations in the interactioninduced polarisability anisotropy tensor [24] _These fluctuations can be written in terms of a sum ofcomponent two-, three- and four-particle density correltition function ensemble averages. It was inferred that the three- and four-particle contributors have a range

(27)

where A = 2j3kg T21V

(28)

2nd p=wv(UP,

(29)

P=PS+PC,

(30)

where (31) and f=f

(32)

A further resolution offiS into two, /3; and three, p3, particle components and PCinto three, 05 and four, pi, particle terms was achieved, ps=Pi +pg,

(33)

where 4 = ;(Q$, - co,i’>,

i#j,

p; = $(c@ij~~il>- <@$($.1?),

i#j,i#

P”=L$ f$, where

(34) I,j#

1, (35) (36)

Note,

P; =o;>

(38)

p~=a((~ij~kkl>-‘~~)ii>(41’), i+j,k+l,i#k, i# List

I,kfj,(39)

Each index implies a sum over ail particles excepr for the pairs indicated. The evaluation of these particle components does not add signi.tlcantly to the computational time as they all can be obtained from ~3~.fl$, and ~3which are readily evaluated along with the forces. The DHX expression for the specific heat,

cHX

= d [SEX

] /dT,

was evaluated numerically and its values are compared with the MD C, in table 1. Interestingly, the simple Debye-Hiickel specific heat, CfH, where, C,DH/fiB =

d[-&Ku/(1

f KG)]

/dT,

W)

which at T = 0.5964 E/kg gives,

CFH/kB = 0.54275(4;;~~)~~*/[1

f 1.29485(4iipo)1’2]2, (41)

has a densi dependence which agrees rluite well with that of C,DE& . For example, at concentrations of

0.009 11,O. 1,0.49.5, 1.0 and 3 M the CFH/kB (em/ ~3~) values are 0.043 (0.045), 0.088 (0.091), 0.104 (O.i03), 0.101 (0.099) and 0.086 (0.083). Perhaps the most si,@ificant conclusion to be drawn from table 1 is that both ps and PC result from the small differences between two large numbers, which thus gives rise to poor statistical accuracy for C,. The twcami four-particle terms are positive whereas the threeparticle terms are negative numbers. These combine to give a 0s which is positive and a PC which is negative. As N increases PC grows in magnitude most rapidly. Most of the uncertainty comes from the cross term which despite being smaller in magnitude than the self term consistently has several times the statistical error. As the three-particle term ccntributes to both self and cross components then the four-particle term is the main cause of poor statistics in this fluctuation property. The coulombic contribution to the specilic heat, a.1 k, is small when compared with the kinetic component which has a value of 1.5 kB _ The CzHX follows the MD and MC specific heats quite closely for c 2 0.5 hi. However, the simulation values continue to rise for hi&her concentrations, because the DJ!IX model insufficiently accounts for the extent of structuring in these media. Some dynamical aspects of the MD simulated systems are considered below. The self-diffusion coefficients, DMD (where

Table 1 energy component of the specific heat per particle, C,,ikg, obtained from MD, MC [SI and DH theory [71. A resolution of the components of the Coulomb energy fluctuations, A& where A = 2/(3k&N), is given for the hID results. The standard deviation of the hiD C&B b denoted by o. p = 0.0924646 c -ChiD chfC CDHX C 4; 4: ‘W (0) ‘W 4 ” -x-K--.(mol kB kB kB kg kEr kB kg X-B kB lo3 m3)

The ccnfig~mioml

0.00911

0.1 0.495 1.0 2.0 3.0 3.0 =) 3.0 b,

0.151

0.631 1.589 2.160 4.147 4.884 7.713 11.15

-0.101 -0.565 -1.491 -2.056 -4.008 -4.665 -7.512 -10.95

0.020

0.067 0.097 0.105 0.140 0.140 0.201 0.195

0.006 -0.012 -0.033 -0.028 -0.048 -0.046 -0.149 -0.132

0.026 0.055 0.064 0.076 0.092 0.094 0.052 0.063

0.040 0.087

it.103 0.124 0.152 0.155 0.083 0.101

(0.020) (0.008) (0.011) (0.005) (0.010) (0.073) (0.945) (0.065)

0.051 0.114 0.136 _ -

0.045 0.091 0.103 0.099 0.090 0.083 0.083 3.083

ON.

162

Heyes/Simdations of primitiue electrolytes

structure about the ions have a neghgiiie effect on iOn self diffusion. The equiv&nt conductances, Alvrr, (where q%-%&o = 1937 $2~~cd eq-t) were obtained by applying an electric field E 2 1 ELT-~@ (or ? MV cm-l) along the direction of one,x, of the cube sides [28]. ~ompiiiirium MD was usedto obtain AMD in favour of a route via a Kubo relation using the current -current correlation function [13] because of the former method’s superior statistical properties. Care must be exercised in analysing the results cf this method, however, The charges drifted with an average speed in this direction, U, of up to 1 m-112P1~2 (or 204 ms-l): The AhID were obtained from the following expression, C / 103

mol

me3

narD = @i/E.

The concentration denendance of the inversediffusion - . ~oefrlclentsLT1/01-1/2,-1/2 obtaimd bj- hiD direcrly (symbols as for fii. 2) and as predicted by Enskog theory of eq. (42),

(43)

Fig. 4.

* using the MD&,(o) mdg&) for the lowest N for eac‘n density-.Thefiledin symbols(otherwiseas for fg. 2) are

(kS~~~~lD)-~/q-Zo’r’/~~-~12.

During this process kinqtic energy was removed from each particle at a rate, Q, by scaling the velocities (in excess of P for the x component) in order to maintain a constant average temperature. This leads to an alternative expression for AhfD as follows, A”D

um-1%1/2

cm2 s-l) were ca!culate.Z frcm the slope of the mean square displacements with time 1131 and are shown in fig. 4. Das 1263 has considered self-diffusion in 2 dilute system of charged particles in a medium of uniform non-zero viscosity using the Smoluchowski diffusion equation. Here the viscosity of the medium is zero 2nd consequently it is more convenient to snalyse the h!D results in terms of a theory apprcpriate to dilute gases. At such low number densities the modified Enskcg theory of hard sphere dynamics [27] gives 2 diffusion coefficient, DE, where (DE)-’

= &2&l

= 57 X IOe5

= [4P/3(kBT)‘/7-]I;

k”(U)

+ g&)1

cJ2rP)

(42) which is in good agreement with DhlD. Consequently, the Coulomb forces have little direct ir . I.zence on the spheres’ dynamics. This suggests that sclote-solute interactions in dilute ionic solutions have an exceedingly sma!l effect on the concentration dependence of this tr2nsport coefficient. When the solvent is “removed” the model system behaves dynarnicaliy as a dilute hard sphere gas and self-diffusion c2n be predicted from the equilibrium structure. As the experimental diffusion coefticients nre also nearly concentration independent this shows that solute indr?ced changes in the solvent

_

(441

As hhS decreased with increasing E, a linear extrapolation to zero E was made. The conductivities decreased with field increase_ For example for fields of 0.16 2nd 0.79 eaS1qS1 at i M (IV= 64) AXrD= 1.7 and 1.4 q2.e-1/2m-1~2a,respectively. Similarly for E = 0.45 2nd 2.27 at 3 M (N = 216) then AMD = 0.52 and 0.45, respectively. There is some non-linear field dependence but such a fit was not warranted as it would not give conductances significantly diiferent from the lowest field values at each density. The range of field used was sufficiently large that it limited the mnge of values for any reasonable extrapolation to zero field. Both formulae for AhiD gave values which agree to within 1% for most states. According to the Nernst-Einstein relationship [ 11, D-’ = q’/kBTA.

(45)

Fig. 4 reveals that the measured AbID are about 15% below those predicted using DLm and eq. (45). The failure of the Nernst-Einstein relation is to be expected. This is because the process of electrical conduction opposes the equilibrium structure established by the Coulomb interactions to a greater extent than tb2t of self-diffusion [7] . This becomes noticeable at j-r&$ field strengths. For example, fields of 0,O.z and 0.9

D,*f.

&ye&&u&ions

e&qS1 at2M(N=64)gavegu(o)of2.2,2.0and 1.7, -UINkgT of 0.72,0.65 and 0.61, and PV/Nk,T of 134, 135 and 1.48. The hamiltonian used in these simulations does not explicitly include ion-solvent interactions. The charges move in a frictionless medium. Solvent effects have only an indirect infh~ence by modifying the ion-ion interactions. Thus solute-somte excess properties are all that ten be derived from these calculations. Hence, the derived transport properties have no parallel in a real solution; only their relative values are of significance. In this context it is particularly interesting that the self-diffusion coefficients and electrical conductance are closely related_ It is widely accepted that an electric field in the Xdirection creates an egg-shaped distortion of the ionic atmosphere in solutions [7] _This can be represented by an anisotropic pair distribution function, &), which on expansion in terms of the equilibrium distribution, g(r) [3_9], is given by [30] _ 0g(r)= g(r) + 6x dg(r)/ax f ... .

W

(47) where ex is a constant. Let the angular and radial integrations be performed on (x/&o). If xB is the x component of rji, the relative positions of spheres i and j, then the ensemble average, l~j(X~/r~>/~ I, about charge i equals ex/S. The strain, eX, measured during the simulations was statistically zero. Thus because of the low number density above (or equivalently absence of solvent particles) the system is unable to sustain any measureable anisotropic distortion which is characteristic of the perturbation (e.g., as is observed when p = 0.7 Na3/V [30] ).

Acknowledgement Professor Dr. L.V. Woodcock is thanked for helpful discussions. Part of this work was performed at the University of Amsterdam, which is thanked for a generous supply of computer time.

of primitiw? electrolytes

163

References 111G.I. S&Z, K. H&zinger and G. P&rlo%, Chem. Phys. Letters 78 (1981) 194. [2] M. Mezei and D.L. Beveridge,J. Chem. Phys. 74 (1991) 6902. [3] D.M. Heyes and J.H.R. Clarke, J.C.S. Famday Trans. II 75 (1979) I240. [4] D.C. Rapaport and H.A. Scheraga, Chem. Phys. Letters 78 (1981) 491. [S] H.J.C. Bereadsen, J.P.M. Postma, W.F. van Gunsteren and J. Hermans, in: Jerusalem symposia on quantum chemistry and biochemistry, ed. B. PuIIman, (Reidel, Dordrecht, 1981). [6] P. Turq, F. LanteIme and D. Levesqoe, Mol. Phys. 37 (1979) 223. [7] J.O’M. Bockris and A.K.N. Reddy, Modern electrochemistry, Vol. 1 (Plenum, New York, 1970). [8] D.N.Card and J.P. VaIIeau, J. Chem. Phys. 52 (1970) 6232. [9] J-P. VaUeau and L.K. Cohen, J. Chem. Phys. 72 (1980) 593.5. [lo] W.J. Moore, Physical chemistry (Lon,gman, London, 1962). [ 111 SK. Jalota and R. Paterson, J.C.S. Faraday Trans. I 69 (1973) 1510. [12] B. Larsen. Chem. Phys. Letters 27 (1974) 47. I131 J.P. Hansen and I.R. McDonsld, Theory of simple liquids (Academic Press, New York, 1976). I141. J.H.R. Clarke, J.F. hia&re and L.V. Woodcock, Faraday _ Discussions &em. So; 69 (1980) 273. [IS] RX Stmtt. S.L. Holmgren and D. ChandIei, hfoi. Phys. 42 (1981) 1233. [16] J.P. VaIleau, L.K. Cohen and D.N. Card, J. Chem. Phys. 72 (1980) 5942. [ 171 E. Waisman and J.L. Lebowitz, I. Chem. Phys. 56 (1972) 3093. [ 181H.C. Anderscn, D. Chandler and J.D. Weeks, J. Chem. Phys. 57 (1972) 2626. [ 191 L. Blum, Theor. Chem: Advan. Perspec. 5 (1980) 1. [20] V-B. Parker, Thermal properties of aqueous m-d-univalent electrolytes (N.B.S. Report, US, 1965). 1211 J.C. Rasaiah, J. Chem. Phys. 52 (1970) 704. [22] R.A. Robinson and R.H. Stokes,_Electrolyte solutions (Butterworths, London, lS59) p. 483. [23] J.H.R. Clarke and J. Bruin& Chem. Phys. Letters 80 (1981) 42. 1241 J_H.R. Clarke and L.V. Woodcock, Chem. Phys. Letters 78 (1981) 121. [25] J.L. Lebowitz, J.K. Percusand L. Verlet, Phys. Rev. 153 (1967) 250. 1261 AX. Das, Chem. Phys. Letters 54 (1981) 249. [27] J.F. Clarke and M. McChesney, The dynamics of real gases (Butterworths, London, 1964). [28] B.R. Sundheiro, Chem. Phys. Letters 60 (1979) 427. [29] J.E. Enderby and G.W. Neilson, Rep. Progr. Phys. 44 (1981) 593. [30] DM Heyes, JJ. Kim, C-l. hfontrose and T.A. Litovitz, J. Chem. Phys. 73 (1980) 3987.