J. Phys. Chem. 1988, 92, 7007-7012
7007
Numerical Analysis of Dlffusion in the Presence of Nonlinear Reverslble Trapping W. D. Murphy, Science Center, Rockwell International, Thousand Oaks, California 91 360
H. M. Rabeony, and H. Reiss* Department of Chemistry and Biochemistry, University of California, Los Angeles, Los Angeles, California 90024 (Received: July 18, 1988)
This paper provides accurate numerical solutions of “stiff“, nonlinear, diffusional boundary value problems involving diffusant which can be bound in a thermodynamically reversible manner to “traps” within the diffusion field. The physical system to which this model corresponds could, for example, be a film of conducting polymer within which the diffusant is an electroactive dopant. This diffusant is also distributed in a thermodynamically reversible manner between the external phase and the polymer at the film surface. Thus full thermodynamic reversibility is possible for the equilibrium absorption (isotherm) of the dopant. Nevertheless, the numerical solutions indicate, as was suggested and demonstrated approximately in an earlier paper, that great disparities can arise between the times required for diffusion-in and diffusion-out. Furthermore, sharp steps and points of inflection occur in the concentration profile of diffusant so that these artifacts (actually observed in conducting polymers) cannot be taken as evidence for dopant absorption which is not thermodynamically reversible.
Doping of Electroactive Polymers A recent paper’ addressed the problem of doping conductive polymers in a thermodynamically reversible manner. In particular, the isotherm for the absorption of iodine vapor by polythiophene, at 50 OC, when measured during desorption showed no difference from the isotherm measured during absorption. This absence of hysteresis, coupled with other observations,2 indicated the achievement of thermodynamic equilibrium with respect to the distribution of iodine between the vapor and the polymer. However, the times required for desorption were dramatically longer than those for absorption. In many cases where observation is p ~ s s i b l ethe , ~ concentration profiles of electroactive dopants, diffusing into conductive polymers, exhibit a sharp, slow moving step or a pronounced point of inflection. These various phenomena have sometimes been taken as an indication that thermodynamically reversible equilibrium could not be achieved. In ref 1, in addition to the confirmation of reversibility implicit in the observed absence of hysteresis, some theoretical analysis was performed to demonstrate that, even with reversibility, steps, points of inflection, and disparities between absorption and desorption times could occur. This theory was approximate and based on the existence of reversible “traps”. For example, a possible but by no means demonstrated model for the trapping of iodine in polythiophene could be 312 + 2e-.PT 213-.PT (1) where I2 refers to free iodine dissolved in polythiophene, e-.PT to an electron in the valence band of polythiophene (PT), and IFPT to a triodide ion bound (in some way) to the polythiophene (PT) chain. Iodine in this complex could be immobilized (trapped) so that it could not diffuse while the free Iz could diffuse. The equilibrium constant relationship corresponding to eq 1 may be expressed as [12] [e-.PT] [I3-.PT]
= K
where the square brackets denote concentration and K is the equilibrium constant. The smaller the value of K, the smaller will be the concentration of the I2 diffusant, and the slower the overall diffusion process. (1) Kim, Dai-uk; Reiss, H.; Rabeony, H. M. J . Phys. Chem. 1988, 92, 2673. (2) Reiss, H.;Kim, Dai-Uk. J . Phys. Chem. 1986, 90, 1973. (3) Handbook of Conducting Polymers; Skotheim, T. A., Ed.; Dekker: New York, 1986; pp 63, 64.
0022-365418812092-7007$01 S O / O
In ref 1, the theoretical analysis was not performed specifically on the system iodine-polythiophene (and therefore not exactly with eq 1 and 2) but rather on a more simple model system, capable of demonstrating the various effects thus far mentioned. In this system eq 1 is replaced by A+T*AT
(3)
where A is the free diffusant, T the trap, and A T the diffusant-trap complex. Equation 2 is replaced by ii(N - 8 ) -a (4) B
--
(This is eq 2 of ref 1.) In eq 4, ii, 8, and N are the concentrations of A, AT, and total traps, respectively, and N - 8 is therefore the concentration of T. a is the equilibrium constant. Eventually ii, 0, and a were normalized by division by N so that eq 4 was converted into u(1 - v ) --a V
(eq 8 of ref 1) where u, v, and a are the normalized quantities. When N greatly exceeds the total concentration of diffusant, v must be neglible in comparison to unity, and eq 5 becomes u =
ffu
(6)
This linear relation between u and v, Le., between free and trapped diffusant, leads essentially to a “chromatographic” mode of migration in which the effective mobility or diffusion coefficient is proportional to u / ( u + v ) and cannot account for steps, disparities in times for reversal, etc. The occurrence of these phenomena requires the full nonlinear relationship, eq 5 . The analysis in ref 1 was analytical, but only approximate. For example, it made use of the physically obvious fact that in the presence of very strong trapping (small a),diffusant would have to “fill” all local traps before it could pass through a region containing unoccupied traps. In absorptive diffusion this gave rise to steps, except in the limit that eq 6 held. Again, appealing to the physically obvious, in the case of desorptive diffusion, diffusant would have to traverse a zone, adjacent to the surface of a polymer film, low in diffusant content so that trapping would dominate and the effective diffusion coefficient would be very low. Thus the “skin”of the polymer would constitute a slowly permeable region and diffusion-out would be accordingly slow (much slower than diffusion-in). Here the authors of ref 1 also envisioned a step and based their approximate analysis on its presumed existence. However, as we shall see below, such steps do not occur during desorptive diffusion. Indeed, more careful physical con0 1988 American Chemical Society
7008 The Journal of Physical Chemistry, Vol. 92, No. 24, 1988 sideration also indicates that they should not occur. We discuss this and other issues below, where we provide accurate numerical solutions to these difficult (and “stifP4) boundary value problems. The Boundary Value Problems
Assuming that Fick‘s law applies to the diffusion of the free species, the relevant partial differential equation governing diffusion is (7)
where D is the constant diffusion coefficient, x and t represent position and time, respectively, while
c=u+v
(8)
denotes the rota1 local concentration of diffusant. We are of course considering one-dimensional, plane parallel diffusion; i.e., the x direction is normal to the surface of the film into which the dopant is diffusing. In ref 1, eq 5 and 8 were used to reexpress eq 7 in terms of u only. For our purpose it proves convenient to express the equation in terms of c only, again with the help of eq 5 and 8. It is also convenient to introduce the dimensionless variables, and T , in place of x and t , i.e.
5 =x/L T
(9)
= Dt/L2
(10)
where L is the thickness of the film. Then eq 7, for c ( ~ , T )becomes , ac _ AT
For the case of absorptive diffusion (diffusion-in from the film surface at = 0) the initial condition is c ( 5 m = 0,
5 >0
(12)
while the boundary conditions are C ( 0 , T ) = Co
(
=0
For desorptive diffusion (diffusion out of the film) the initial condition is C(5,O)
= co, 5
>0
(15)
and the boundary conditions are C(0,T)
=0
(16)
The boundary value problem, for diffusion-in, thus consists of eq 11, 12, 13, and 14 while, for diffusion-out, it consists of eq 11, 15, 16, and 17. We have solved both of these problems, for various sets of the parameters a,co, by a numerical method described in the Appendix. These solutions are discussed in the next section. Solution and Discussion
Before discussing particular solutions, it is appropriate to estimate what we might expect to see on the basis of physical considerations. It is useful to itemize these expectations as follows. A . Diffusion-Zn. (1) When co > 1, Le., when there can be more diffusant locally than is required to fill all the local traps, there will always be a considerable excess of free diffusant. If a is very (4) Gear, C . W . Numerical Initial Value Problems in Ordinary Differential Equutions; Prentice-Hall: Englewood Cliffs, NJ, 1971.
Murphy et al. small, traps will always be filled locally until only a very small number of unoccupied traps remain, and no diffusion past such a position will be possible until this degree of filling has occurred. Thus the concentration profile is expected to advance as a sharp step, movement of the step taking place when the traps behind it have been filled. (2) When co < 1, there can never be an excess of local diffusant over traps, even when a is very small. Thus the amount of free diffusant allowed by eq 5 will be very small, although larger a t smaller values of 5, where the concentration of total diffusant is larger. Thus diffusion will occur more readily at smaller values of 5 than at larger values. The result will be an accumulation (a backing-up) of diffusant at smaller t so that the concentration profile will exhibit a point of inflection. (3) In general diffusion will be slower when a is smaller, since there will be less free diffusant. B. Diffusion-Out. (1) Because of eq 16 the region of the film adjacent to its surface will always have a total concentration (free plus trapped) of diffusant, low in comparison to the concentration of traps. As a result, and especially with small values of a , most of the diffusant in the region will be trapped and immobilized. Thus the region will behave as a slowly permeable membrane in comparison to regions at larger values of t. To a first approximation, at larger values of t, all gradients will have time to level out, so that the only remaining gradient will be a steep one in the slowly permeable zone near the film’s surface. The diffusant then “leaks out” in such a manner that the concentration profile is always horizontal behind the slowly permeable zone while its level is uniformly lowered with time in the manner of reservoir behind a dam whose floodgates are slightly opened. (2) Again, in general, diffusion will be slower when a is smaller, since there will be less free diffusant. (3) In general, the times for diffusion-out will exceed the times for diffusion-in. The reason for this is the following. Consider first the case where co > 1.O. In this case, and for diffusion-in, there is a “perpetual” supply of diffusant at t = 0 at a level exceeding the number of traps. Thus, even when the traps are filled, there will be free diffusant available to diffuse “over” them. Thus diffusion-in can be relatively fast. On the other hand, even though co initially exceeds the concentration of traps throughout the film, for the case of diffusion-out, since there is only a finite supply of diffusant its concentration must eventually fall below that of the number of traps. Then, a large fraction of diffusant will be bound and immobilized so that the times for diffusion-out will be relatively long. (4) As c, falls below unity the times for absorption and desorption tend to approach one another, since there is little free diffusant in either case. However, the perpetual supply at 5: = 0, in the case of absorption, tends to keep the time for absorption shorter than that for desorption. We now turn to the various numerical results which provide examples of the above items, A (1, 2, 3) and B (1, 2, 3, 4). Figures 1, 2, and 3 illustrate concentration profiles, c versus 6, during diffusion-in (absorption) at reduced times T = 0.1, 1, 10, and 100, respectively, for cases in which the respective values of co are 1.2, 1.0, and 0.8. In each figure we show profiles (dot-dashed line) for two values of a,namely a = lo4 and lo-’. Curves (dot-dashed line) for v, the concentration of trapped, immobile diffusant are also shown. In Figure 1, for co = 1.2 > 1.O and a = lo4, the sharp step is clearly evident, in accordance with the discussion of item A.l, above. For the larger value of a,namely a = lo-’, a large enough fraction of diffusant remains untrapped so that a sharp step is not observed. The vertical distances between corresponding dot-dashed and full curves measure the concentration u of free diffusant. Figure 3 illustrates the situation when co = 0.8 < 1.0 and confirms the discussion in item A.2. Again values of a = and lo-’ are involved. In accordance with item A.2 the step has vanished in favor of a point of inflection. Again, the dot-dashed line shows the concentration profile for the trapped diffusant. However, for the a = lo4 case, that profile is almost indistinguishable from the full curve depicting the profile for total diffusant. Figure 2 represents an intermediate case, co = 1.0. In all three figures, item A.3 is confirmed, Le., the profiles for the smaller
The Journal of Physical Chemistry, Vol. 92, No. 24, 1988 7009
Numerical Analysis of Diffusion
P
7 --.-
O
0. 4 I
0.2
-
0.0 0.0
,
I
0. 2
0. 1
0. 3
"\ 0. 4
0. 5
E Figure 1. Profiles of concentration c(5,7) with a = lo4, a = lO-l, and co = 1.2 at some intermediate time 7 during absorption. The dot-dashed
line is the concentration 457) of trapped diffusant.
-.- ,
0.0
0. 1
0. 2
0.4
0.3
0.5
E Figure 4. Profiles of concentration 457) with a = lo4 and co = 1.2 at
some intermediate time 7 during desorption. The dot-dashed line is the concentration u(5,7) of trapped diffusant.
.-.
1 \
0.2 -
0.0
!
4 0.0 0. 1 0. 2 0. 3 0. 4 0. 5 0.0
0. 0
€ Figure 2. Profiles of concentration c ( [ , T ) with a = lo4, a = and co = 1.0 at some intermediate time 7 during absorption. The dot-dashed line is the concentration u ( [ , 7 ) of trapped diffusant.
0. 1
1
1
0. 2
0.3
I
0.4
0.5
E Figure 5. Profiles of concentration c([,7) with a = lo4 and co = 1.0 at
some intermediate time 7 during desorption. The dotdashed line is the concentration u ( [ , 7 ) of trapped diffusant.
a-10-4,r100
0.21
0. 0
\ 0. 1
'-\ 0.2
0. 3
-
0. 4
0. 5
E Figure 3. Profiles of concentration c ( f , 7 ) with a = lo4, a = lO-I, and co = 0.8 at some intermediate time 7 during absorption. The dot-dashed line is the concentration ~ ( 5 , s of ) trapped diffusant.
(a = lo-") always lag behind those for the larger a (a = lo-'), Figures 4, 5 , 6, 7, and 8 illustrate concentration profiles for diffusion-out (desorption). Again, dot-dashed line shows the concentration profiles for trapped diffusant. In Figure 4, co = 1.2 and a = lo4. Profiles for the reduced times T = 1, 10, and 100 are shown and confirm, in accordance with item B.1, that at larger values of F, the gradient levels out and that the level tends to be lowered uniformly with the passage of time. Also, in accordance with B.l, there are no points of inflection. This is in contradistinction to the presumption in ref 1. Figure 5 provides a similar series of profiles for the case co = 1.0 and a = lo4, at time 7 = 0.01, 1, and 100. (Y
0.04 0. 0
I
I
0. 1
0.2
0.3
0. 4
0.5
E Figure 6. Profiles of concentration c ( f , 7 ) with a = 10-4 and co = 0.8 at some intermediate time 7 during desorption. The dot-dashed line is the
concentration u ( [ , 7 ) of trapped diffusant. Figure 6 again provides a similar series of profiles for the case co = 0.8 and a = lo4, this time at T = 0.1 and 10. Figures 7 and 8 show profiles for a = 10-I and co = 1.0 and
0.8, respectively, at T = 0.01 and 0.1. All profiles for diffusion-out have essentially the same shape, again in accordance with item B.1. These five figures confirm item B.2. For example, in Figure 4 where (Y = 10-4 we have to deal with values of T as large as 100 where as in Figure 7 and 8 where a = IO-', the largest values of T are only of the order of 0.1, Figures 9, 10, and 11 address the question of the disparities in the times for absorption and desorption. In these figures Q(T)
7010 The Journal of Physical Chemistry, Vol. 92, No. 24, 1988
Murphy et al.
-0.
n
-1.
E?
m 0
-1.
-2.
y
0.0 0.0
0.2
-2.5
t
I
1
0.1
0. 3
0. 4
0. 5
E Figure 7. Profiles of concentration c(5,7) with a = lo-' and eo = 1.0 at
some intermediate time 7 during desorption. The dot-dashed line is the concentration u(&7) of trapped diffusant.
I
I
I
-5
I
I
I
6
-1
-3
I
I
5
3
1
log (r)
Figure 10. Amount Q of diffusant absorbed and desorbed versus 7 with co = 1.0; the upper curve for each pair, a = lo-' and a = lo4, corre-
sponds to absorption and the lower to desorption.
1
0.8
0.6
--
'i:
0.4
v
0
0.2
-2.5
0.0
0.2
0.1
0.4
0.3
0.5
E
:
-5
I
I
-3
I
I
-1
I
I
1
I
I
I
3
log (r)
Figure 8. Profiles of concentration c(5,7) with a = lo-' and c, = 0.8 at
some intermediate time T during desorption. The dot-dashed line is the concentration u(5,7) of trapped diffusant.
Figure 11. Amount Q of diffusant absorbed and desorbed versus T with co = 0.8; the upper curve for each pair, a = lo-' and a = lo4, corre-
sponds to absorption and the lower to desorption. TABLE I: Reduced Time T Required To Achieve C ( ] / * , T ) = 0 . 9 for ~ ~ Different Values of co and a during Absorption Ly ~n = 0.8 cn = 1.0 Cn = 1.2 lo-' 1 1 1 10-2 10 10 1 10-3 100 10 1 10-4 1000 100 10
n
8 m 0
3
-2.5
1
I
-5
-3
I
I
-1
I
1
1
,
I
3
I
5
log (r)
Figure 9. Amount Q of diffusant absorbed and desorbed versus 7 with co = 1.2;the upper curve for each pair, a = lo-' and a = lo4, corre-
sponds to absorption and the lower to desorption. is either the amount of dopant absorbed in time 7 (if absorption is under examination) or the amount of dopant desorbed in time 7 (when desorption is under consideration). The curves plot log Q versus log 7 . In each figure, two pairs of curves are shown, Figures 9, 10, one pair for a = lo-' and the other for a = and 11 differ in having co = 1.2, 1.0, and 0.8, respectively. In each pair of curves, the upper member corresponds to absorption and the lower to desorption, thus confirming item B.3 which states that times for absorption should be smaller than those for desorption. Furthermore, as expected, the disparity increases as a decreases, but decreases (the members of a pair become closer) as co is decreased. This confirms item B.4.
TABLE I 1 Reduced Time T Required To Achieve C ( ' / ~ , T ) = 0 . 1 ~ ~ for Different Values of co and a during Desorption Ly Cn = 0.8 Cn = 1.0 Cn = 1.2 10 10 10-1 10 10-2 100 100 100 lo-' 1000 1000 100 10000 10000 1 Od 10000
Tables I and I1 provide sharper views of the trends in the times for absorption and desorption. Since the processes are asymptotic, the time for absorption is defined as that for which c at 5 = 1 / 2 attains the value o.9c0, while for desorption the criterion is c at = 1 / 2 attains the value 0 . 1 ~ ~Table . I refers to absorption and Table I1 to desorption. The numbers in these tables are the times according to the above criteria. It will be noticed that, for a = IO4 and co = 1.2, T for absorption (Table I) is only 10 as compared to 10000 for desorption (Table 11). On the other hand, if co is reduced to 0.8 while a remains at lo4 the time for absorption is increased to lo00 while that for desorption remains at 1Oo00. Thus the ratio of times is 1000 for co = 1.2, but only 10 when co = 0.8. This tendency for the two times to approach one another, as co is decreased, is in accordance with item B.4, and the greater times required for desorption confirm item B.3.
The Journal of Physical Chemistry, Vol. 92, No. 24, 1988 7011
Numerical Analysis of Diffusion
Appendix
TABLE III: Diffusion-In Cases“ ff
1 (-1)
co
T
NQ
AT
0.8
l (-14) 1 l (-14) 1 1 (-14) 1 1 (-14) 10 l (-14) 10 1 (-14) 1 1 (-14) 100 l (-14) 10 l (-14) 1 l (-14) 1000 l (-14) 100 1 (-14) 10
1 3 1 4 1 4 1 4 1 2 1 4 1 4 1 2 1 4 1 1 1 1 1 1
1.00 (-14) 6.03 (-2) 1.00 (-14) 7.64 (-2) 1.00 (-14) 8.30 (-2) 1.00 (-14) 5.60 (-1) 1.00 (-14) 4.78 (-1) 1.00 (-14) 4.89 (-2) 1 (-14) 9.77 l (-14) 4.65 l (-14) 3.69 (-2) l (-14) 654 l (-14) 212 l (-14) 10
1 1.2 1 (-2)
0.8 1 1.2
1 (-3)
0.8 1.0 1.2
1 (-4)
0.8 1.0 1.2
NSTEP NFE NJE 10 146 10 158 10 184 10 180 10 247 10 317 10 190 10 390 10 1187 10 195 10 733 10 3313
23 279 23 309 23 379 23 371 23 636 23 866 23 430 23 1144 23 3597 23 462 23 2443 23 8899
4 31 4 33 4 35 4 36 4 48 4 59 4 40 4 81 4 290 4 40 4 191 4 919
The numerical procedure for integrating eq 11 is similar to that discussed in Reiss and Murphy? For completeness we give some of the details and indicate where care is necessary. Basically, eq 11 is discretized by the method of lines; that is, only the spatial derivatives are differenced. For example, the interval [0, 1/21 is divided into a sequence of NPTS points such that 0 = 51 < tz < ... < [NPm = 1/2
- ti.Using centered differences, we
and hi is defined as hi write
Ci = C(ti,T) For the boundary condition, eq 14 or eq 17, one-sided differences are employed, i.e.
‘The numbers in parentheses represent powers of 10, Le., 1 (-3) =
(CNPTS - CNPTS-I)
1 x 10-3.
TABLE I V Diffusion-Out Cases’ ff
Cn
T
NO
AT
1 (-1)
0.8
1 (-14) 10 1 (-14) 10 1 (-14) 10 1 (-14) 1 (-10) 100 1 (-14) 100 1 (-14) 100 1 (-14) 1000 1 (-14) 1000 1 (-14) 1000 1 (-14) 10000 l (-14) 10000 l (-14) 10000
1 4 1 3 1 4 1 2 3 1 3 1 2 1 3 1 3 1 2 1 3 1 3 1 3
1 (-14) 0.787 1 (-14) 0.846 l (-14) 0.738 1 (-14) 1 (-10) 8.87 l (-14) 8.60 1 (-14) 12.8 l (-14) 79.5 l (-14) 76.9 l (-14) 133 l (-14) 778 l (-14) 854 l (-14) 923
1.o 1.2 1 (-2)
0.8 1
1 (-2)
1.2
1 (-3)
0.8 1.o 1.2
1 (-4)
0.8 1.o 1.2
NSTEP NFE NJE 10 166 10 169 10 169 10 22 172 10 177 10 183 10 176 10 180 10 193 10 179 10 178 10 200
23 297 23
308 23 301 23 47 314 23 322 23 330 23 321 23 328 23 375 23 327 23 328 23 403
4 35 4 38 4 35 4 8 39 4 40 4 41 4 40 4 41 4 44 4 41 4 42 4 50
“The numbers in parentheses represent powers of 10, Le., 1 (-3) = 1 x 10-3.
Conclusion The model to which the above boundary value problems correspond is a thermodynamically reversible one. The traps establish instantaneous local equilibrium with diffusant, and the fKed values of c(0,T) correspond to reversible solubility coefficients. Yet it is clear that (1) longer times for desorption than for absorption and (2) sharp steps and inflection points for concentration profiles occur. These features cannot therefore be used as evidence for irreversibility. The accurate numerical analysis of the present paper confirms the conclusions arrived at by means of both experiment and a more approximate theory in ref 1.
Acknowledgment. This work was supported in part by the Air Force Office of Scientific Research under Contract F49620-86C0060.
The boundary condition at t; = 0 gives C1 = 0 for eq 16 and C1 = c,, for eq 13. These are known values; however, to simplify our software, the dummy equation dCl/dt = 0 is employed. Substituting eq A1 and A2 into eq 11 yields the following nonlinear stiff4 system of ordinary differential equations: dCl/dr = 0 dCi/dr =fi(Ci-l,CirCi+l,T)
i = 2, 3,
..., NPTS - 1
(A3)
with initial conditions given by
Ci(0) = c(ti,O) (A41 The value of c([,O) in eq A4 is taken from eq 12 or 15 depending upon the type of diffusion being investigated. Unfortunately, standard Adams-Moulton-Bashforth predictor-corrector methods or RungeKutta methods cannot be used to solve eq A3 and A4 because the stiffness of this system will require time steps which are intolerably too small. Instead, a stiff integrator by Hindmarsh: which contains backward differentiation operators of order 1-5, is employed. In addition, this algorithm automatically changes the step size and the order by monitoring the error growth and comparing the differences between predictor and corrector results with the truncation errors for various orders. This integrator uses the Jacobian of the right-hand side of (A3) to form a pseudo-Newton’s method to converge the corrector equation. It makes use of the fact that this Jacobian is tridiagonal. More details of this algorithm can be found in ref 6 and 7. The main advantage of this approach over a more classical method such as Crank-Nicolson is that the step size and order are automatically adjusted to guarantee that a relative error tolerance is automatically satisfied. In the cases run for this paper the relative error tolerance was set to lo4. Since large gradients can occur near the left boundary, this automatic monitoring can be very important. Studies were made on the spatial gridding to determine if special care is needed near the left boundary. Our standard case consists (5)Reiss, H.;Murphy, W. D.J . Phys. Chem. 1985,89,2596. (6)Hindmarsh, A.C.Lawrence Livermore Report UCID-30059,Rev. 1, March, 1975. (7)Sincovce, R. F.;Madsen, N. K. ACM Trans. Math Software 1975, 1 ( 3 ) , 232.
J. Phys. Chem. 1988, 92, 7012-7015
7012
of 51 uniform points on the interval [0, 1/21, Le., h, = 0.01. A second grid consists of hl = h2 = ... = hlo = 0.001 and hll = h12 = ... hS9= 0.01. Both grids produced equally accurate results, and the nonuniform grid was found to be unnecessary in order to guarantee accurate results. The uniform grid has the additional advantage of being second order in the spatial direction. In the time direction, this algorithm is optimally between first and fifth order. Time integration results were printed at T = 10" for n = -14, ) -13, ...,-1, 0, 1, ...,4 or until the conditions on c( I / ~ , T(previously discussed) were satisfied. In Tables 111 and IV we show some of this data for the diffusion-in and -out cases, respectively. We use the notation 1 (-14) = lO-I4. Here N Q equals the current order of the time-step integrator, AT is the current step size, NSTEP equals the number of time steps, N F E is the number of functional evaluations (right-hand side of eq A3), and N J E is the number of Jacobian evaluations. Note that the Jacobian is not
reevaluated unless the previous Jacobian failed to iterate the corrector to convergence in the Newton's algorithm. It is obvious from Table IV that our integrator can easily handle all diffusion-out cases even though the maximum time for diffusion (not CPU time) is larger than for the corresponding diffusion-in cases. Notice that the most difficult problem to solve numerically is the diffusion-in case with parameters a = lo4 and co = 1.2. The difficulty arises from the fact that large gradients are present in the solution. Nevertheless, our integrator automatically adjusts the step size and order to track these large changes. This example clearly illustrates why a robust algorithm is required to solve this rather simple looking parabolic partial differential equation. Also, automatic adjustments in AT are essential for efficient integration because in some cases the maximum value of T was 10 000. All cases were run on the VAX 780 in double precision to minimize the effect of round-off errors when integrating for long periods of time.
Intramolecular Ll+ and Na+ Bonds wlth Fluctuating Charge in Dkarboxylates and In a Dialcoholate: A Fourier Transform Infrared Spectroscopic Study Bogumil Brzezinski, Department of Chemistry, A . Mickiewicz University, 60- 780 Poznari, Poland
Georg Zundel,* and Rainer Kramer Institute of Physical Chemistry, University of Munich, 0-8000 Munich, FRG (Received: January 25, 1988)
Li+ and Na+ complexes of bis(tetrabuty1ammonium) salts of tetramethyl-o-phthalic acid (1) and ethylmaleic acid (2) and of the Li+ complex of the bis(tetrabuty1ammonium) salt of o-bis(hydroxymethy1)tetramethylbenzene (3) were studied in solutions by FTIR spectroscopy. Li+ complexes with intramolecular Li+ bonds are formed. Continua in the far infrared demonstrate that with the complexes of 1 and 2 the Li+ ions fluctuate in -O-Li+-O- bonds. In these bonds a broad flat single-minimum potential well is present. In the complex of compound 3 the Li+ ions fluctuate probably in -OLi+-.O-O-.Li+O- bonds with a double-minimum potential. These bonds interact strongly with their environments, caused by the so-called Li+ polarizabilityof these bonds. In the case of the Na+ complexes a very broad absorption is observed in the region 150-25 cm-I. In -O-Na+-O- bonds the amplitude of the Na+ ion fluctuation is much smaller than that of Li+ ions in Li+ bonds due to the larger mass of Na+. Thus, the Na+ polarizability of the bonds is much smaller than the Li+ polarizability. The same is true with regard to the interaction of these bonds with their environments.
Introduction Caused by the motion of protons in hydrogen bonds with or with m ~ l t i m i n i m apotentials, ~.~ polarizado~ble-minimuml-~ bilities occur. These proton polarizabilities of hydrogen bonds are 2 orders of magnitude larger than the usual polarizabilities, caused by a distortion of electron systems. In the IR spectraG8 such hydrogen bonds cause continua. In the case of longer hydrogen bonds with double minimum such continua extend from (1) Weidemann, E. G.; Zundel, G. 2.Naturforsch., A : Phys., Phys. Chem., Kosmophys. 1970, 25, 627. ( 2 ) Janoschek, R.;Weidemann, E. G., Pfeiffer, H.; Zundel, G. J. Am. Chem. Soc. 1972, 94, 2387. ( 3 ) Eckert, M.; Zundel, G. J . Phys. Chem. 1987, 91, 5170. (4) Fritsch, J.; Zundel, G.; Hayd, A.; Maurer, M. Chem. Phys. Lett. 1984, 107, 65. ( 5 ) Eckert, M.; Zundel, G., J . Mol. Struct. (THEOCHEM), in press. (6) Zundel, G. In The Hydrogen Bonds-Recent Lkwlopments in Theory and Experiments; Schuster, P., Zundel, G., Sandorfy, C., Eds.; North Holland Amsterdam, 1976; Vol. 11, Chapter 15, p 681. (7) Zundel, G.; Fritsch, J. In The Chemical Physics of Solvation; Dogonadze, R.R.Kllmafi, E., Kornyshev, A. A,, Ulstrup, J., Eds.; Elsevier: Amsterdam, 1986; Voi. 11, pp 21ff. (8) Zundel, G. In Biomembranes, Protons and Water: Structure and Translocation ( A Volume of Methods in Enzymology): Packer, L., Ed.; Academic: New York, 1986; Vol. 127, Part 0, p 439.
0022-3654/88/2092-7012$01.50/0
3000 cm-' toward smaller wavenumbers, whereas short hydrogen bonds with a broad flat potential well in which the proton fluctuates cause continua in the region 1700-600 ~ m - ' . ~ ,Such ~ continua are found with charged B+H-B * B-H+B bonds as well as with noncharged AH-aB == A--.H+B bonds. If Li+ ions fluctuate in a Li+ bond, large so-called Li+ polarIn the case of NLi+-.N izabilities of these bonds are * N-.Li+N bonds, for instance, these polarizabilities are indicated by continua extending from 450 cm-' toward smaller wavenumbers. Until now only Li+ bonds with noncharged acceptors were studied in which the ion fluctuates. In the following, bonds with charged acceptors are studied in which the cation fluctuates, whereby Li+ and Na+ bonds are compared.
Results and Discussion Li+ Bonds. If salts are dissolved in organic solvents, the so-called ion motion band is usually observed in the region 500-100 (9) Brzezinski, B.; Zundel, G. J . Mol. Strucr. 1981, 72, 9. (10) Brzezinski, B.; Zundel, G. J . Chem. Phys. 1984, 81, 1600. (11) Bnezinski, B.; Zundel, G. J . Chem. Soc., Faraday Trans. 1 1985,8J, 2380.
0 1988 American Chemical Society
