Molecular dynamics simulations of highly concentrated salt solutions: Structural and transport effects in polymer electrolytes

SOLID STATE

Solid State Ionics 53-56 (1992) 1011-1026 North-Holland

IOHICS Molecular dynamics simulations of highly concentrated salt solutions: Structural and transport effects in polymer electrolytes M. Forsyth ~, V.A. Payne, M.A. Ratner, S.W. de Leeuw 2 and D.F. Shriver Department of Chemistry and Materials Research Center, Northwestern University, Evanston, IL 60208, LISA

Structural, thermodynamic and transport properties have been calculated in concentrated non-aqueous Nal solutions using moleculardynamics simulations. Althoughthe solvent has been represented by a simplistic Stockmayerfluid (spherical panicles with point dipoles), the general trends observed are still a useful indication of the behavior of real non-aqueous electrolyte systems. Results indicate that in low dielectric media, significant ion pairing and clusteringoccurs. Contact ion pairs become more prominent at higher temperatures, independent of the dielectric strength of the solvent. Thermodynamic analysis shows that this temperature behavior is predominantly entropicallydriven. Calculation of ionic diffusivitiesand conductivities in the NaI/ether system confirms the clustered nature of the salt, with the conductivities significantlylower than those predicted from the NernstEinstein relation. In systems where the solvent-ion interactions increase relative to ion-ion interactions (lower charge or higher solvent dipole moment), less clustering is observed and the transport properties indicate independent motion of the ions, with higher calculated conductivities. The solvent in this system is the most mobile species, in comparison with the polymer electrolytes where the solvent is practically immobile.

1. Introduction Solvent-free solid polymer electrolytes have been pursued as potential substitutes for fluid electrolytes in a variety of electrochemical applications. Since the original observation that polyethylene oxide ( P E O ) / salt complexes have substantial ionic conductivity, significant advances have been made experimentally in the preparation and characterization of new materials. Advances in understanding the mechanism of conduction and ionic structure within these polymer electrolyte systems have been less rapid. From several spectroscopic and conductivity studies, it is generally accepted that significant ion pairing occurs in the polymer electrolytes. For example, McCallum et al. [ 1 ] and others [2 ] have measured conductivity as a function of concentration in low molecular weight poly (propylene oxide) (PPO) and poly(ethylene o x i d e ) ( P E O ) from very dilute sys-

tems to the more concentrated solutions more frequently used in polymer electrolyte applications. A broad conductivity m i n i m u m is observed between 0.01 and 0.i mol kg-1 salt concentration and then a m a x i m u m at much higher concentrations above 1 mol kg- J. The m i n i m u m was suggested to be a result of the formation of pairs at moderate concentrations that do not contribute to the conductivity, and then to the formation of larger, charged aggregates at higher concentrations that may transport charge. With increasing salt concentration, the overall viscosity of the solution increases and the mobility of charge carriers decreases. Another explanation proposed for the occurrence of the conductivity minim u m , and one that is often seen in explanations of this m i n i m u m in other electrolytes [ 2 - 4 ] , is that at higher ionic concentration, ( > 0 . 1 M) redissociation of ion pairs and aggregates occurs, resulting in free ions.

Current address: P.O. Box 50, Ascot Vale, Victoria 3032, Australia. :'- Laboratoryfor Physical Chemistry,Nieuwe Achtergracht 127, 1018 WS Amsterdam, The Netherlands.

Two extreme conduction mechanisms have been postulated in the literature. Cheradame and others [5 ] have suggested that it is p r e d o m i n a n t l y the free ions that carry the current and that at higher temperatures, more free ions will result. MacCallum et

0167-2738/92/$ 05.00 © 1992 ElsevierScience Publishers B.V. All rights reserved.

1012

M. Forsyth et al. / Molecular dynamics simulations

al. [ 1 ] suggest that triple ions and larger aggregates carry the current since no free ions will exist in such concentrated solutions. Recently, a paper by Boden et al. [6] discussed the use of pulsed field gradient (PFG) NMR techniques to measure the diffusion coefficients of both cation and anion in PEO/salt complexes at elevated temperatures (well above the melting point of PEO). From these measurements, the expected conductivity that would result if the system consisted of uncorrelated motions of free ions can be calculated from the Nernst-Einstein equation. These authors found that at very high salt concenlrations (2 to 4 tool kg -~) and lower temperatures, the Nernst-Einstein fit to experimental conductivity measurements was excellent, however the fit was not so good at lower concentrations ( ~ 1 mol kg- ~~ 2 mol kg- 1). It was postulated by the authors that it is indeed the single-ion motion which is responsible for ionic conduction but that the equilibrium picture may consist of larger clusters and may be approximated by the structure of a molten salt. Their suggestion is that the single-ion charge carriers result from a cooperative configurational rearrangement of clusters of ions. There are also conflicting interpretations of the structural changes in these polymer solutions as a function of temperature, Cheradame et al. [ 5 ] have suggested that increased conductivity results at higher temperatures due to the increased activated concentration of free carriers ( N ~ c e x p ( - E B / k T ) ). As the temperature increases the binding energy between ions can be overcome and more free ions result. Vibrational spectroscopy has been used extensively to observe ion pairing in both liquid and polymer electrolytes [ 7-15 ]. Kakihana et al. [ 7 ] and Teeters and Frech [15] have observed by vibrational spectroscopy, that i n c r e a s e d ion association occurs as the temperature increases, contradicting Cheradame's suggestions. The explanation given is that the decrease in entropy which occurs when the salt is dissolved in the polymer (resulting from "virtual" crosslinks that restrict the segmental and translational modes of the polymer) is recovered at higher temperatures by the ions associating less with the solvent and more with each other. Eventually this can lead to a "'salting out" or precipitation of the salt at very high temperatures [3]. Nitzan and Ratner [ 16] have presented a simple configuration discus-

sion that rationalizes this picture in terms of solvent configurational degrees of freedom. Despite these extensive experimental investigations, the theoretical understanding of ionic conductivity, the precise nature of the ionic species and the types of solvent-ion interactions that control the structure and conduction in these systems, is still incomplete. In this work, we present the results of molecular dynamics simulations (MD) of concentrated ionic solutions. The solvent in these systems is modelled by Stockmayer particles (i.e. spherical Lennard-Jones particles with a point dipole). Although this is not a realistic model for polymer electrolytes, we are primarily interested in the behavior of the ion in a dipolar fluid with dielectric properties close to that of the host polymers used in polymer electrolytes. Work is also underway to model the solvent by more realistic ether molecules and oligoethers. Molecular dynamics computer simulations allow control of physical variables such as temperature, ion concentration and the solvent dielectric constant. In this work a systematic investigation of the dependence of structural, thermodynamic and transport properties on these variables is discussed and a microscopic picture of concentrated ionic solutions is developed. Radial pair distribution functions and visual snapshots of the particles in the MD simulation box are presented in the following sections and allow the determination of the structure in the concentrated ionic solutions as a function of the important experimental parameters. The free energy of an ion pair in these solutions as a function of interionic separation (also called the potential of mean force) [17,18] is calculated, From the temperature dependence of this free energy, the entropic and energetic contributions to the formation of contact and solvent separated ion pairs is discussed. Finally, the diffusion coefficients of the ions and solvent molecules are determined using both mean-squared displacement and velocity autocorrelation function calculations. The conductivity is obtained from the Nernst-Einstein equation using these computed diffusion coefficients and compared to the conductivity calculated from the current autocorrelation function. Insights into the transport mechanism can be derived from such comparisons.

M. Forsyth et al. / Molecular dynamics simulations

2. Computational details M D simulations were carried out on a system consisting o f N Z ions and N D solvent molecules at fixed t e m p e r a t u r e and density. The physical p a r a m e t e r s presented in table 1 a p p r o x i m a t e a N a I / e t h e r system. The potential energy o f the system, taking into account Lennard-Jones, i o n - i o n , i o n - d i p o l e and dip o l e - d i p o l e interaction, is given by eq. ( 1 ): V ( R j .... , R , , p l , -..,/thD, Zl .... , ZNZ) = l (ND+NZ) (ND+NZ)

Z

Z

t=l

Jv~i

OcJ ( IR,s [ )

1 ND ND

"~- 2 i~=l j~i ODD(Rij,~'gi,tgJ) NZ ND

+ E Z 0,D(R0, Z,,~,) i=lj=l

1

NZ NZ

+~ i=~.I j~i~ ¢)~ou,(IR,~l,Zz,Z~).

(1)

The Ewald sum was used to evaluate the coulombic interactions using " t i n f o i l " conditions [19,20]. Reduced units were throughout, thus all lengths are measured

long range boundary employed in units o f

1013

the Lennard-Jones size p a r a m e t e r for the solvent ac and all energies are calculated in units o f the LJ well depth o f the solvent, tc. Reduced t e m p e r a t u r e is defined as T * = k B T / e c and the unit o f time is given by r = ( m c a ~ / e c ) '/2 .

Reduced density p * = 1 was used throughout. This is equivalent to a real density o f 2.7 g / c m 3. Simulations were run with 256 panicles, however, the number dependence o f the results was tested by carrying out simulations with 500 and 1098 particles. The resuits a p p e a r e d to be i n d e p e n d e n t of box size. Table 2 summarizes the computer experiments that will be discussed. The dielectric constant, t, was calculated in our systems using the K i r k w o o d formula [ 19,20] with tinfoil b o u n d a r y conditions, e = 1 + 3yg where y = 4 n p o l F / 9 k T , / t = dipole m o m e n t , Po = solvent density and g = ( M 2 ) / N D I t 2. ( M 2) is calculated throughout the simulation using M=ENO/u,. Note that this form can be used to estimate the dielectric constant contribution from the solvent even for a conductive system, where the overall ~ (¢o= 0) is not measurable. The potential of mean force, which is defined as the free energy o f an ion pair as a function of interionic distance, was calculated from the

Table 1 Physical parameters used in simulations. Particles

a (A)

~ (K)

q (e)

p (D)

mass (a.m.u.)

Na I solvent

1.9 4.2 3.07

58.2 420 90

+½; + 1 - I; - 1 0

0 0 1.3; 2

23 126.9 44

I

0.05

Table 2 Systems studied. Exp.

No. of solvent molecules

No. of ion pairs

Conc. (M)

T*

P*

/~ (Debye)

Iq l (e)

time steps

A B C D E F a)

240 240 224 248 232 240

8 8 16 4 12 8

1.8 1.8 3.6 0.9 2.7 1.8

2, 2, 3, 4 4 2,

l 1 I 1 1 1

1.3 1.3 1.3 1.3 1.3 2

~ l 1 1 1 1

100000 100000 100000 100000 100000 20000

a) With higher dipole moment only static properties were obtained.

3, 3.5, 5 3, 4 4

3, 4

l 014

M. l"orsyth el al. / Molecular c!vnamtcs simulation,s

pair distribution function for 8 ion pairs and 240 solvent molecules.

(rdf's) for N a - N a , Na-I, I-1, N a - s o l v e n t and l-solvent interactions. The total number of particles at r', d e t e r m i n e d b y the r u n n i n g i n t e g r a t i o n , is also s h o w n

for selected temperatures. These arc defined as 3. Structural properties

N,,b(r') =4~pb i

Figs. 1-4 present the radial distribution functions

gab( r)r2dr"

(2 )

(I

5

10

S

8

4

(a)

4

(b)

8

----

"I'-2 3

- -

1"=3

. . . . .

I'=5

6

2

4

../ .~""

....

I=

6 N

5

2

. ..../

0

0

0 0

05

1

1.5

R

2

25

(~ ~,)

10

35

3

f

",

12

/'

- -

1"=3

-

1'=5

....

15

2

25

3

35

20

(d)

3

2

0 05

15

- - 1 ' = 2

0 0

1

16

•

4

05

20

(c)

/

6

0

R 1o /~l

,

/ Z

I *~2

I* 3

3 N

z

10

2

25

35

3

0

0 0

05

1

15

2

25

3

35

R (a. ~,l

R (o c A) 15

6

(e)

5

/

I0

4

3

.>7

5

¸,..=2 2

1

O

0 0

05

1

15

R

2

(o<. ~,)

25

3

35

Fig. 1. Radial distribution functions and running integration as a function of temperature (7"*=2 ( 180 K), T*=3 (270 K), T*=5 (450 K)) in experiment A: 8 ion pairs, 240 solvent molecules, :L= 1.3 Debye, p*= 1 and q= ~-e. Distances are scaled by the Lennard-Jones parameter of the solvent, ~ . = 3.07,4. The relatively larger noise level for the ion-ion curves is due to fewer particles, hence to poorer statistics.

M. Forsyth et al. / Molecular dynamics simulations

20

!

16

~

12

z

7

(a)

..'~J

6

24

5

20

4 N

&

z --

:'.

4

3 --Ti=2

"F*=3

8

i

0,5

1.5

1

R

2 (a c

2.5

3

0

0 05

0

20

16

1.5

2.5

35

3

20

16

4

=

12

3

tI . T*=2

,"

N i

--

--

8

T*=3

0

1

2

(d)

~--~

2

I *=2

- -

-

....

o.5

1.5

5

-

0

1

(C)

•

2

0

1

R (c, A)

l!

J

z~

4

T*=4

A)

,

3

1

35

I,'

~

2

__ Ti=3

1"=4

~ 0

6

5

....

0

.~f~"

~16 'T

2 . . . .

. (b)

3

T*=2 --

7

28

./

.-2 '.

1015

2

2.5

3

4

1

0

0

3,5

T*=3 1"*=4

'

~ 0

0.5

1

o 1.5

2

2.5

3

35

R (o A) 8

80 70 60 5O

5

& 4o

a

p. ½

30

N

3

20

--

--

--

-- I"*=3

10

1":2

2

i??i

R(~ A)

Fig. 2. Radial distribution functions and running integration as a function of temperature ( T* = 2 ( 180 K ), T* = 3 (270 K ), T* = 4 (360 K ) ) in experiment B; 8 ion pairs, 240 solvent molecules, /1= 1.3 Debye, p * = 1 and q = 1 e.

where Pb is the number density of species b, gab(r) is the rdf at distance r. These results are most clearly understood with the help of visual snapshots obtained from the MD run at T = 2 7 0 K for 8 ion pairs with ( q = ½,/~= 1.3 D), ( q = 1,/z= 1.3 D) and ( q = 1,

p = 2 D), (fig. 5). The effect of strong ion-ion interactions (fig. 5b) compared to stronger ion-solvent interactions (fig. 5a and 5c) are clearly shown here. In fig. 5a the ions are fully solvated despite the relatively weak dipole strength of the ether-like sol-

o

0 0

0.5

1

1.5

2

25

3

35

1016

M. Forsyth et al. /Molecular dynamics simulations

15

7

// ,/ / ' /

&-

10

/

z z

/

/

/

/

(a)

',

6

l',,'

5

--~.=3 -- - 1.=4

10

7

/ /

II

/

/

(b)

__::::::

/

s

4

3 5

5

,

2 1

0

0 0

0.5

1

15

R

2

25

3

0 ~)

05

1.5

(o A)

R

il

(c)

2.5

3.5

3

(o, ,&) (d)

15

20

1"*=3 //

'/"

2

2

~lII

20

l

I

& Z

0

35

--

--

I6

T*=4

3

8

12

2

N

,8

4

0

0 0

05

1

15

2

25

/•

3

3.5

1

t5

2

(o ~

25

3

35

.&)

(e)

iI

I

40

05

R

R IOc ,~)

50

0 0

It 3O i

T.=3 -- "1*=4

20

10 1

l

o

" o

o.s

~ ~

~ R

vent molecules tween the ions charge is le, as are very strong weaker dipoles

o

~

~.s

2

2.s

3

3.s

(o c ~x)

since the coulombic interaction beis also relatively weak. When the in N a + I - , the i o n - i o n interactions and are not fully o v e r c o m e by the of the solvent molecules. These ef-

Fig. 3. Radial distribution functions and running integration as a function of temperature ( T * = 3 (270 K), T*=4 (360 K)) in experiment C; 16 ion pairs, 224 solvent molecules,/~= 1.3 Debye, p * = l , q = l e.

fects are more quantitatively illustrated by the rdf's. In fig. 1 the N a - N a and I-I distribution functions do not show any significant structure but just level off to g(r)= 1. A peak is observed in the cation-anion distribution function at r--0.9 and another between

M. Forsyth et al. / Molecular dynamics simulations

1017

5

5

(a)

(b) 4

' --T*=3

T'=3 - -

T'=4

z

- -

T'=4

....

T" =5 3

3

-T*=5

;=

N

N

z 2

2

2

•

1

0

0

0

0.5

1.5

R (a

2

2.5

3.5

3

"

"

a

o

o.s

I

2

2.s

3

3.5

2O

5

20

(c)

10

1.s

R (o'¢ A)

A)

15

e.

0

.

(d)

T'=3 - -

16

4

12

3 •

T'=4

16

" T ' =5

' --T'=3

u~

- -

T'=4

....

T =5

N

12

N

u~ ,

z

8

8

4

0

0 0

1.5

05

2

2.5

3

0

35

0.5

1

1.5

2

2.5

3

3.5

7

6

- -

,F~ = 3

- -

T

. . . .

T " =S

(e)

=4

4

5 3

4

N

z 3

~"

•

2

1

R (oA)

Fig. 4. Radial distribution functions and running integration as a function of temperature ( T *= 2 ( 180 K), T* = 3 ( 270 K), T *= 4 (360 K)) in experiment F; 8 ion pairs, 240 solvent molecules, /1=2 Debye, p*= l and q= l e.

r = 1.8 and 1.9, w h i c h correspond to contact ion pairs and solvent separated ion pairs respectively. T h e r u n n i n g integration s h o w s that overall ion pairs are few in n u m b e r . The two peaks o b s e r v e d in both the c a t i o n - s o l v e n t and a n i o n - s o l v e n t rdf's correspond to

the first and s e c o n d s o l v a t i o n spheres. T h e effect o f temperature on ion pairing in this system is also interesting. It appears from fig. 1c, that ion pairing is reduced as temperature decreases. From pure energetics, one m a y expect that as the temper-

0

'

0

0.5

1

1.5

2

.

.

.

2.5

.

.

3

0

3.5

1018

M. k~)rsyth et al. / Molecular ctvnamics .simulatiml,S . ~

100 I 8M, ~ - Z D

80

36M,~I

3D

0 2 M , ~ - I JD

•

--~

] 8M, ~ . 1 3 D , p ° i 5

60

40

\

•

(a)

~ A

150

232 5

- -

___

315

397 5

480

Temperature (K) Fig. 6. Static dielectric constant as a function of temperature, salt concentration, and dipole moment.

{b)

Fig. 5. (a) Snapshot of ionic structure, 7"= 270 K, q= ~ e, ll= 1.3 D. The larger spheres represent the anions and the smaller spheres represents the cations. Two periodic MD boxes are shown here to clearly show clustering. The solvent has been left out for clarity: (b) snapshot of ionic structure: T=270 K, q= 1 e, ll= 1.3 D, 8 ion pairs; (c) snapshot of ionic structure; T=270 K, q= 1 e, I~= 2 D, 8 ion pairs.

ature increases, ion pairs w o u l d be less f a v o r e d as suggested by C h e r a d a r n e . O n e possible e x p l a n a t i o n for this decreasing ion pairing with d e c r e a s i n g t e m p e r a t u r e m a y be f o u n d in the b e h a v i o r o f the dielectric constant. Fig. 6 illustrates the t e m p e r a t u r e be-

h a v i o r of e as calculated u n d e r various c o n d i t i o n s and clearly shows an inverse t e m p e r a t u r e d e p e n d e n c e o f e. H e n c e one possible e x p l a n a t i o n tot the b e h a v i o r o b s e r v e d in fig. I may simply be that the dielectric screening is m o r e effective at lower temperatures. Fig. 2 shows results for larger q. This enhances i o n ion interactions twice as m u c h as ion-solvent. Fig. 2c shows that the c a t i o n - c a t i o n c o n t a c t pairs do not appear to decrease with d e c r e a s i n g t e m p e r a t u r e . T h e r e is actually ve~, little effect o f t e m p e r a t u r e in this case. T h e o t h e r r d f ' s in fig. 2 show that clustering is occurring in this system at all t e m p e r a t u r e s , with N a N a peaks, and l - I peaks present at the distances app r o p r i a t e for an fcc lattice. T h e r u n n i n g integration suggests that, on average, two to three anions and two solvent m o l e c u l e s are c o o r d i n a t e d to each cation. Increasing the c o n c e n t r a t i o n in this system with q = 1 and i t = 1.3 D, does not significantly change the structural b e h a v i o r . Significant clustering will still o c c u r as is seen f r o m the i o n - i o n r d f ' s in fig. 3 and the snapshot in fig. 7; h o w e v e r , the t e m p e r a t u r e dep e n d e n c e in this m o r e c o n c e n t r a t e d system is m o r e d r a m a t i c than that in fig. 2. O n c e again increased ion association is o b s e r v e d with increasing t e m p e r a t u r e , with the a n i o n c o o r d i n a t i o n about the s o d i u m increasing by one as the t e m p e r a t u r e is raised from 270 K to 360 K. Finally, in o r d e r to o b s e r v e the effect o f solvent strength on the ionic structure, the d i p o l e m o m e n t was increased to 2 D. T h i s c o r r e s p o n d s (fig. 4) to a

1019

M. I~brsvth et al. / Molecular dynamics simulations

Fig. 7. Snapshot of ionic structure; T=270 K, q= 1 e, It= 1.3 D, 16 ion pairs.

dielectric constant of 80 at T = 270 K. The effects of increasing the size of the dipole moment are quite dramatic. The ions are fully solvated at the lowest temperatures and no clustering is apparent. However, as the temperature is raised, the N a - I rdf shows that the contact ion pair peak becomes more prominent and the magnitude of the solvent separated peak at an ion separation of 1.8 ac is decreased. The integration shows that the actual proportion of contact pairs is still low, nevertheless it is interesting that again ion association increases with increasing temperature. The dielectric constant in this system is still high even at 360 K. Experimentally, increase of ion pairing with increasing temperature is not only observed in polymer electrolytes, as discussed by Kakihana et al. [7,8], but also in simple solvent systems. For example, Nal in 2-butanone reaches a maximum in solubility at - 30 °C and then decreases at higher temperatures [21 ]. Similar effects are observed in acetone and ethanol. It has been argued that in polymer electrolytes, the increased ion association at higher temperatures is an entropic affect due to the loss of low-frequency modes of the polymer when it is coordinated to the salt [ 18 ]. It is equally probable that entropic considerations may explain the increased association in simpler solvents and in our simple Stockmayer model. The ions can structure the solvent by causing the dipoles to orient themselves; while dipolar ion pairs also structure the local solvent, they do so less effectively due to the weaker (effectively dipole-dipole) pair/solvent interaction [22]. This then decreases the degrees of freedom of

the solvent molecules and so at higher temperatures, it is more favorable for the ions to associate and allow the dipolar solvent molecule to translate and rotate freely. An alternative explanation, based on the effect of temperature on the dielectric constant is also likely. In systems containing ions of lower charge or higher dipole strength solvent molecules, the solvent-ion interaction will be comparable to the i o n ion interactions. As the temperature is increased and the dielectric medium is weakened, ion association becomes more probable. In the next section we investigate the thermodynamics of solvation and determine the entropic and energetic contributions at different temperatures and dipolar strengths.

4. Thermodynamic

properties

The interionic potential of mean force ( p m f ) can be defined through the pair distribution function by gab(r)=exp

-- k T , / '

(3)

where W~b(r) is the p m f at interionic separation r. In our systems, the p m f is equivalent to the Helmholtz free energy, A ( r ) , of the system as a function of average separation of a given unlike ion pair, i.e. it corresponds to the work that would be required to separate a pair from a given separation, r, to infinite separation. From the temperature dependence of the p m f w e can therefore determine the entropic and energetic contributions. That is,

0A(r)] ~J~,v

= - S(r)

and A(r)=

U(r)-

TS(r) .

(4)

Fig. 8 is a comparison o f p m f ' s for 8 ion pairs with differing ion charges and dielectric constant. The temperature dependence in each case is also shown. In the case where dipole-ion interactions are strong relative to the interionic interactions (fig. 8a), two minima are clearly seen at the contact ion pair and solvent separated ion pair distances. In the case of a higher solvent dipole moment, fig. 8c, the first minimum does not appear energetically favorable. In-

1020

M. Forsyth el al. / Molecular dynamics simulation~

20

(a) 15

10 O

E

O.

0

i

~..,e -5

,,e

i

It/ -10

0.5

1

1.5

2

2.5

3

3.5

distance (~c) fi 7-

t

to

•

•

(b)

11':

- i

i;

lo

-15

t!

~}-0

- ~ - -

0.5

1

I

I

1.5

2

2.5

3

35

distance (~) 10

(c)

I

8 6 0

,I

'gr-*d

I--

,iX ¢",:N

4

2

I

T1450K

,~t •

i

0

T.27oKI

L--,-,.0,1 -- - -

¢,~

i

•

"

-2

, i

-4

I I II

lo

e

l

-,'." %/I "~_

! i B i i

--4

"6

0.5

I

1.5

2

2.5

3

3.5

distance ( ~ ) Fig. 8. Potential of mean force for 8 Na-I pairs as a function of temperature; (a) q =-~ e,/z = 1.3 D; (b) q= 1, # = 1.3 D; (c) q= 1, It = 2 D. Calculated from the pair distribution function, gNA-~,for systems with 8 ion pairs.

deed the ~econd m i n i m u m has a lower free energy in this case, when T = 2 7 0 K ( - 2 k c a l / m o l at r = 1.9 c& c o m p a r e d to + 1 k c a l / m o l for the contact ion pair distance). As the temperature increases, the contact m i n i m u m and the solvent separated m i n i m u m occur at shorter interionic separations and the contact pair becomes energetically more favorable, with the p m f d r o p p i n g to - 5 . 5 k c a l / m o l at T = 4 5 0 K. A third m i n i m u m is also observed in this system, probably indicating a second solvation sphere about each ion. The corresponding snapshot for this system, shown in fig. 5c, clearly indicated that free ions arc more a b u n d a n t than ion pairs. In the case of q = 0 . 5 e and weaker dielectric medium, the same temperature dependence is observed, with contact pairs becoming more favored at higher temperatures and the barrier to the formation of a solvent separated ion pair increasing. When the charge is increased to q = 1 e in the system, the contact pair becomes much more favorable and indeed at higher temperatures, instead of a m i n i m u m being observed at about 1.8 a~, (solvent separated pair) two m i n i m a corresponding to ionic positions in the fcc lattice appear (fig. 8b). Ions are obviously strongly attracted to one another in this case as was observed from the snapshots in fig. 5b. The entropic and energetic contributions to the p m f were d e t e r m i n e d from the temperature dependence as discussed above (eq. ( 4 ) ) and are illustrated for the case of q = 0 . 5 e, # = 1.3 D in fig. 9. It is evident that at all temperatures, the entropic contribution favors increased ion pairing. The entropy gain of the solvent when a contact pair is formed overrides the unfavorable energetic term. Indeed the b a r r i e r to the formation o f a contact ion pair appears to be largely energetic in origin. Similar behavior is observed in both other cases with the charge set to q = 1 e, however, the behavior displays slightly more noise. The above analysis has confirmed the ideas proposed by [7,8] that enhanced ion pairing at higher temperatures is a result of the entropic gain of the system. It is also in agreement with analytical results presented by Pettitt and Rossky [23] for the case of a single ion pair at infinite dilution in water. This behavior is thus not unique to polymer electrolytes. Even in a simple Stockmayer fluid, the dissolution of ions can result in ion-solvent coordination thereby removing translational and rotational degrees of

1021

M. Forsyth et al. / Molecular dynamics simulations

freedom and so decreasing the entropy of the solvent. In more complicated solvents there are some additional losses in entropy resulting from losses of librational and vibrational modes. At elevated temperatures, the drive to increase the solution entropy favors the disruption of ion-solvent interactions and results in an increasing number of contact ion pairs.

(a)

o

° o ° oo°°

E oo

5. Transport properties I

*5 0

~1

0.5

1

I

I

I

I

1.5

2

2.5

3

distance

3.5

(~c)

6

(b) 4 o 2

o o

Diffusion coefficients and conductivity were calculated from simulations that were generally run for 100000 time steps. In order to improve the statistics, sliding average techniques were used to calculate the velocity and current auto correlation functions for the ions. The correlation functions are defined by



%



0

E

and -2

(j(t)j(O) > -4



-6

where

-8 0

I

I

I

I

P

I

0.5

1

1.5

2

2.5

3

all ions 3.5

j=

Y~ q, virespectively.

The latter function includes ion-ion cross correlations. The self-diffusion coefficient can either be obtained by integration o f the velocity auto correlation function

20

o

(c)

15

10

'i

D= ~

5 o

E

(5)

i=l

distance (~c)

oo

(v(l)v(O)>dt,

(6)

0

0

or

from

the

mean

squared

displacement > where t = 0 is reset every 2000 time steps in this work. The self-diffusion coefficients of the cation, anion and solvent using both of these methods are given in table 3. Fig. 10 is an example of mean squared displacement data for q = 1 and q = 1 at T = 2 7 0 K. The frequency dependent conductivity can be obtained by Fourier transforming the current auto correlation function

(r2(t) --r2(O) >, using 2tD= ~(r2(t)

-5

-10

-15

o

I

I

I

P

I

I

o.s

~

~.s

a

a.s

a

a.s

distance (a) Fig. 9. Thermodynamic decomposition of the potential mean force for 8 NaI pairs with q=0-5 and,u= 1.3 D; (a) T=180 K, (b) T=270 K, (c) T=450 K. Calculated from the temperature dependence of the pmf (fig. 8a).

1022

M. Forsyth et al. / Molecular d y n a m w s simulations

Table 3 Calculated diffusion coefficients. Exp.

{from eq. (6)}X 10-9m2 s

Temp. (K)

[from ( r Z ( t ) > l x l 0 ~m-~s '

D+

D

D,

D+

D_

D•

180 240 315 450

1.27 2.77 3.68 6.37

1.07 2.25 2.14 3.32

1.76 3.14 3.8 6.09

1.02 2.5 6.05

0.85 1.7

1.65 3.0

3

6.03

B

180 270 360

0.22 0.5 0.81

0.31 0.62 0.9

1.81 3.42 -

0.23 0.49 0.85

0.24 0.53 0.79

1.82 3. I 0 4.56

C

270 360

0.43 0.44

0.61 0.52

4.13

0.42 0.35

0.42 0.34

2.39 3.93

D

360

1.05

1.35

4.82

1.16

1.17

4.89

E

360

0.68

0.85

4.32

0.66

0.69

4.14

A

2

• 1.6

•

^

12

•

©

;°

©

~o

©

~

04

x×X

x

X XX x

XXX

X 0 0

1

2

3

4

time(x2.354picoseconds)

*• •

1.5

exp(io)t)dt,

(7)

w i t h t h e dc c o n d u c t i v i t y g i v e n at ( o = 0 . Fig. 11 pres e n t s t h e v e l o c i t y a u t o c o r r e l a t i o n f u n c t i o n s for a s y s t e m o f 8 i o n pairs, 240 s o l v e n t m o l e c u l e s at T = 2 7 0 K, in a d d i t i o n to t h e c u r r e n t c o r r e l a t i o n f u n c t i o n , a n d fig. 12 gives t h e F o u r i e r t r a n s f o r m s o f these. F r o m fig. 12 t h e f r e q u e n c y d e p e n d e n t c o n d u c t i v i t y a n d d i f f u s i o n c o e f f i c i e n t s c a n be calculated. O n l y t h e dc c o m p o n e n t s o b t a i n e d f r o m t h e int e g r a t i o n o f t h e c o r r e l a t i o n f u n c t i o n s are p r e s e n t e d in t a b l e 4. In a d d i t i o n , t h e c o n d u c t i v i t y c a l c u l a t e d from the Nernsl-Einstein relation and the individual ion d i f f u s i v i t i e s is also g i v e n in t a b l e 4 t o g e t h e r w i t h t h e d e v i a t i o n s , A, f r o m w h i c h t h e d e g r e e o f ion cross-correlation can be determined [6,24].

×x xX

•

(j(t)j(O)) 0

©

• 0 0

0.8

1 3 V/% T

©

• ©

~((o) =

©

• • •

(a)

(~)

A

'1.

1 Ne 2

v

cr= - - -

(D+ +D

)(1-/1).

(8)

2 VkBT 0.5

0

I

2

3

4

5

time (x2.354 piceseconds)

Fig. 10. Comparison of mean squared displacements as a function of time at T=270 K, for system with ( a ) q = 1/2 e and

(b) q= 1 e, at an ion concentration of 1.8 M. (Experiments A and B ).

T h e f r e q u e n c y d e p e n d e n t F o u r i e r t r a n s f o r m c a n be r e l a t e d to t h e v i b r a t i o n a l m o d e s o f t h e s y s t e m w i t h t h e p e a k at ~ 2 4 0 c m - ~ in t h e F T o f t h e s o d i u m a u t o correlation function reflecting the vibration of this c a t i o n in its s o l v e n t or a n i o n " c a g e " . T h i s is in agreem e n t w i t h low f r e q u e n c y v i b r a t i o n s f o u n d in N a l solutions [25]. In c o n c e n t r a t e d s o l u t i o n s , s t r o n g i o n - i o n a n d i o n s o l v e n t c r o s s - c o r r e l a t i o n s will be p r e s e n t a n d t h i s will

1023

M. Forsyth et al. / Molecular dynamics simulations 1.5

T

velocity

I

I

eutocorreletlon

time

W

for

i

1.2

cetJonc

(;q)

velocity

eutocorreletlon

time

for

(b)

onions

1 I

0.8

A o

>-

0.5

A

0

~.

~v

0.6

"A

0.4

~ o.2

V -0.5 -0.2 i

I

I

I

0.1

0.2

0.3

04

i

i

-0.4 0.5

0.1

time (x2.354 picoseconds)

1.5

i

¢u~rlnl

i

eulocorrllltlOn

v

0.2

i

l

0.3

0.4

05

time (x2.354 pico~eonds)

12

i

(c)

function

,

1

~

veloclty

v i eutocorreletloft

time

r for

eolveflt

, moleculee

(d)

0.8 0

O.S

"%

06

0.4

o

g

~V 0 2 -0.S

-1

-02

0

0.1

02

0.3

04

0.5

I

I

I

i

o.1

0.2

03

0.4

0.5

time (x2.354 picoseconds)

time (x2.354 picoseconds)

Fig. 11. Velocity autocorrelation functions for the cation, anion and solvent molecules at T=270 K, q = 1 e and 8 ion pairs. The current correlation function is also shown.

25I

20

40

(a)

(b)

35

II 3O

"~

is

t 20

"~

15 10

5

0

100

200

300

400

co (cm-l)

500

600

700

800

0

too

200

3oo

4oo

see

I

1-

600

700

800

re (¢m-l)

Fig. 12. Fourier transform of the autocorrelation functions at T = 2 7 0 K as a function of frequency for 8 pairs with (a) q=0.5 e and (b) q = 1 e. In both cases the current response reflects that of the cation.

M. Forsyth et al, /Molecular dynamics simulations

1024

Table 4 Exp.

Temp. (K)

a (S/m)

A

180 270 315 450

4.64 6.96

B

180 270 360

0.72 1.81 1.64

C

270 360

4.45 2.35

12.5 7.7

0.64 0.7

D

360

1.50

6.5

0.77

E

360

2.26

I 1.3

0.8

8.88

a,.~]c (S/m) 5.22 7.81 10.1 5.24 7.59 9.15

J

0.11 0.11 0.12 0.86 0.76 0.82

the clusters moving rather than individual ions. To investigate the mechanism of ion transport, i.e. to see whether the observed conductivity is related primarily to the m o v e m e n t of the charged clusters or whether there are "free" ions moving from one cluster to another, the ion coordinates were m o n i t o r e d as a function of time and saved every ten time steps. The visualization of these coordinates as a function o f time indicated that no hopping from one cluster to another occurs. Instead the cluster is rearranging and ions move a r o u n d and within the cluster. There arc several other ways to m o n i t o r the m o t i o n of the ions with time. This was done in the system with 8 ion pairs and q = 1 for several thousand steps. Fig. 13 is an example of the c a t i o n - a n i o n separation with

5.5

affect the mobility of the ions and the overall solution viscosity. This will be reflected in decreased diffusion coefficients. In addition to that however, crosscorrelations occur which will affect the conductivity of the solutions. This will be reflected in the deviation p a r a m e t e r A presented in table 4. If no cross correlations existed, A would be zero. As expected, few correlations are observed in the system with q = ½ since ions are not so severely clustered and move independently o f one another. This is further illustrated in fig. 10, by c o m p a r i n g the mean squared displacements as a function o f time in the two systems with q = ½ and q = 1. In the case where q = ½, the ions move at different rates and thus the diffusion constant o f the cation is higher than the bulkier anion, whilst when q = 1, the ions move (on average) at the same rate. Boden et al. [6] have d e t e r m i n e d the deviation p a r a m e t e r in concentrated polymer salt systems as a function o f salt concentration and temperature and have found that A increases with increasing temperature and decreasing ion concentration. The magnitude o f 3 that they report suggests strong correlations, probably as a result of clustering. Similar trends can be observed in the data presented here. In other words, the deviation from N e r n s t Einstein behavior is greatest at lower concentrations and higher temperatures. Thus it appears that in our model N a I / e t h e r system, a large degree of clustering is present and the high diffusivities possibly relate to

(a)

.~o< 4s

~

\

/

/

2.5 0

2

4

6

g

10

12

time (x 2.1345 picoseconds) 12

q

r

- - ~ - - - - T

.

.

.

.

.

.

(b)

~°'<

""P"i"

8

4

0

2

4

6

8

10

12

time (x 2.1345 picesecends) Fig. 13. An example of thc time dependence of ion pair separation in a solution of 8 ion pairs, H= 1.3 D at T = 2 7 0 K: (a) q= I e, (b) q=0.5 e. With the smaller chargc thc ion separation increases more rapidly and fewer pairs remain "'in contact" compared with q = 1 e.

M. Forsyth et al. / Molecular dynamics simulations

time, of pairs that start as contact pairs. Most such pairs vibrate but stay within the contact pair separation, occasionally an ion will move out of the first coordination shell and continue to move away. In comparison, when the charge is decreased to q = 0 . 5 e, ion pairs are very few in number and they are also short lived. Fig. 13b shows four cation-anion pairs (with q = 0 . 5 e) and their separation with time at T = 2 7 0 K. Another way to monitor ion motion is to observe the mean squared displacement of individual ions as a function of time. In the systems where q = 1 e a n d / t = 1.3 D, most ions oscillated and moved gradually with time, however occasionally an ion " j u m p s " a considerable distance in only one time step. In the present simulations this would just correspond to hopping from one position on the cluster to another. In actual polymer electrolytes, ion diffusion over distances large compared to local displacements occurs, but the solvent motions are limited by chain entanglements. Solvent "segmental motions" facilitate ion transport, but the solvent particles themselves show much smaller diffusion coefficients [26]. This is a substantial difference from the behavior shown in fig. 10, where the solvent is the most mobile species. The important physical features concerning interionic correlations and their dependence on solvent and temperature obtained in this study should still be relevant to polymer electrolytes.

6. Conclusions The present simplifications o f the solvent system make straightforward comparisons with structure and transport in polymer electrolytes difficult. General trends observed, however, are very useful. It appears that in a low dielectric medium, significant ion pairing and clustering will occur since the solvent-ion interactions are not strong enough to overcome the strong interionic attractions. This is observed even in liquid solvent systems: if we compare the solubility o f Nal in ethers and ketones to that in alcohols [ 12 ], a substantial difference is seen. It should be noted that even in the formation o f polymer electrolytes, salts of low lattice energies are required [27 ]. The temperature behavior of ion pairing in these simulations is in agreement with that observed from

1025

spectroscopic experiments in polymer electrolytes [7,8,15] and can be attributed to the increased entropy o f the solvent upon increased ion association. This is independent of whether the solvent is a strong or weak dielectric. It is also interesting that, although significant clustering occurs, there is no real sign of ion triples in these systems. The diffusion coefficients and conductivities measured in these simulations are relatively high and correspond, as would be expected, to values in a fluid system rather than polymers. Due to extensive clustering in the system with q = 1,/z= 1.3 Debye, the observed conductivity in this case is probably a result of either the cluster mobility or the intra-cluster rearrangement - there is no evidence as yet for intercluster hopping. It is also possible that the solubility limit has been exceeded in this system, particularly at higher temperatures. In the systems where the ions are relatively free, the observed conductivities are probably a result of the transport of single ions. The lack of ion-ion correlations is supported in this latter system by the low value of the correlation factor, 3.

Acknowledgements This work was supported by the N S F / M R C through Northwestern Grant M R L D M R 8821571, and by the A R O DAAL-03-90-G-0044. M.F. and V.A.P. gratefully acknowledge a Fullbright Postdoctoral Fellowship and an NSF Predoctoral Fellowship, respectively. We are grateful to J.W. Perram and to A. Nitzan for incisive discussions; and to M.D. Todd for graphical advice.

References [ 1] J.R. MacCallum, A.S. Tomlin and C.A. Vincent, Eur. Polym. J. 22 (1986) 787. [2] M.C. Wintersgill, J.J. Fontanella, S.G. Greenbaum and K.J. Adami6, Brit. Polym. J. 20 (1988) 195. [ 3 ] C.W.Davies, Ion Association (Butterworths, London, 1962) chapt. 10. [4 ] J. Ludwig, K. Duppen and J.K. Kommandeur, J. Chem. Soc., Faraday Trans. I 80 (1984) 2943. [5] H. Cheradame, in: IUPAC Macromolecules, eds. H. Benoit and P: Rempp (Pergamon, New York, 1982) p. 351.

1026

M. Forsyth et al. / Molecular cJvnamics simulations

[6] N. Boden, S.A. Leng and I.M. Ward, Solid State lonics 45 (1991) 261. [7] M. Kakihana, S. Schantz and L,M. Torell, J. Chem. Phys. 92 (1990) 6271; 94 (1991) 6296. [8] S. Schantz, L.M. Torell and J.R. Stevens, J. Chem. Phys. 94 (1991) 6362. [9J D.E. Irish and M.H. Brooker, in: Advances in infrared and Raman Spectroscopy, eds. R.J.H, Clark and R.E. Hester, Vol. 2 (Heyden, London, 1976) p. 212. [ 10] B.L. Papke, M.A. Rather and D.F. Shriver, J. Electrochem. Soc. 129 (1982)p. 1434. [ 11 ] B.L. Papke, M.A. Ratner and D.F. Shriver, J. Phys. Chem. Solids 42 (1981) 493. [12] D.N. Bhattacharyya, C.L. Lee, J. Smid and M. Szwarc, J. Phys. Chem. 69 (1964) 608. [13 ] C. Carvajat, K.J. Toelle, J. Smid and M. Szwarc, J. Am. Chem. Soc. 87 (1965) 5548. [14] R. Woolridge, M. Fisher, G. Ritzhaupt and J.P. Devlin, J. Chem. Phys. 86 (1987) 4391. [15] D. Teeters and R. Frech, Solid State Ionics 18/19 (1986) 271. [ 16 ] M.A. Ratner and A. Nitzan, Faraday Discuss. Chem. Soc. 88 (1989) 19.

[17] W.L. Jorgensen, in: Advances in Chem. Phys., Vol. LXX (Wiley, New York, 1988 ) p. 469. [ 18] O.A. Karim and J.A. McCammon, J. Am. Chem. Soc. 108 (1986) 1762. [19] S.W. de Leeuw, J.W. Perram and E.R. Smith, Ann. Rev. Phys. Chem. 37 (1986) 245. [20] S.W. de Leeuw, B. Smit and C.P. Williams. J. Chem. Phys. 93 (1990) 2704. [21 ] Nonaqueous Electrolytes Handbook, Vol. II (Academic Press, New York, 1973) Chap. 1. [22] Y. Marcus, Ion Solvation (Wiley-lnterscience, New ~'ork, 1985). [23] B.M. Pettitt and P.J. Rossky, J. Chem. Phys. 84 (1986) 5836. [24] G. Cicotti, G. Jacucci and I.R. McDonald, Phys. Rev. A 13 (1976) 426. [25] B. Guillot, P. Marteau and J. Obriot, J. Chem. Phys. 93 (1990) 6148. [ 26 ] M, Minier, C. Berthier and W. Gorecki, J. Phys. ( Paris ) 45 (1984) 739. [27] D.L. Papke, M.A. Ratner and D.F. Shriver, J. Electrochem. Soc. 129 (1982) 1694.