3454
J . Phys. Chem. 1992, 96, 344-3457
rate of its dissociation. The result for Einstein’s model is similar.
XlO” 6.00
5 00
T--
Concluding Remarks
4
In the real solid, the phonon spectrum is neither Einstein nor Debye distributed and is rather complicated, especially for the surface phonons. An appropriate model suggested by Leibfried and BrenigI2 can be constructed as a mixture of these two distributions
4 00 A
N
. 5 3 00
E(w) = k e T [ 3 w 2 / 2 ~ D 3 H- (W~) ~+ N E / N ~ ( w-Ew ) ]
6 2 00
0.00
r
0
200
400
r
T
s
r
600
n
s
r
,
-
t
,
1
800
,
u
~
1
~
1
,
1
1000
wD ( c m - ’ >
Figure 3. Preexponential factor, A, versus the Debye cutoff frequency of the solid, wD.
different ratios ub/wD.As is visible from both the figure and eq 10, the large wDlimit of A is ATST,and in a small region of wD around about 200 cm-I, A decreases more than five times, i.e. the solids with lower Debye cutoff frequencies give assistance to the energy dissipation from the adsorbed molecule and retarded the
However, it is shown in the present paper that the contributions of Einstein’s and Debye’s distributions in the preexponential term are similar. For this reason we conclude that the influence of the real phonon spectrum is of the same order as this obtained here. In this work the observed delay of the rate constant of adsorbed molecule dissociation generated by energy dissipation from molecule into the solid is only in respect to the TST constant. For a real comparison of the reaction rates in the gas phase and in the adsorbed state, the friction in the former one should be taken into account. Registry No. NO, 10102-43-9; Ni, 7440-02-0.
(12) Leibfried, G.; Brenig, W. Z . Phys. 1953, 134, 151.
A Variational Solution to the Hypernetted Chain Equations Applied to the Electrical Double Layer Scott E. Feller and Donald A. McQuarrie* Department of Chemistry and the Institute of Theoretical Dynamics, University of California, Davis, California 9561 6 (Received: November 22, 1991)
A variational method is described for the solution of the HNC/MSA equations for electrolytes near a charged planar surface. The variational method involves a much simpler computational scheme than has been presented previously for this problem and avoids some convergence difficulties associated with iterative methods. The method is demonstrated with calculations in the restricted primitive model for electrolytesof several concentrations and ionic charges near a surface with a fixed charge density. We also consider a system containing a mixture of monovalent and divalent counterions. The results are compared with the modified Gouy-Chapman theory which is the theory, because of its relative simplicity,most widely used in applications. The variational method may be useful to workers in the many areas where the electrical double layer plays a role.
Introduction
The electrical double layer is of great importance in the study of many problems in biophysics, physical chemistry, and colloid physics. The standard theoretical treatment has been that of Gouy and Chapman, which entails the solution of the Poisson-Boltzmann equation. If the Poisson-Boltzmann equation is solved with the condition that the ions may not come within the Stem layer, which is located one-half of an ionic diameter from the surface, the result is the modified Gouy-Chapman (MGC) approach. Due to its neglect of the finite size of the ions in the bulk of the solution as well as ion-ion correlations, the MGC method fails in many cases. For example, it fails in systems of high surface charge density and high electrolyte concentration, or for solutions containing divalent ions. A more rigorous statistical mechanical approach is to apply the Omstein-Zemicke relation to this system, by taking the limit where one type of particle is a charged plane of infinite extent, and then to use an approximate closure relation between the resulting integral equations.’ Author to whom correspondence should be addressed at the Department of Chemistry, University of California, Davis, CA 95616-8618.
0022-3654/92/2096-3454%03.00/0
The hypernetted chain/mean spherical approximation (HNC/MSA) has previously been applied to the restricted primitive model of a planar double layer with great su~cess.2~The calculated ion distributions show excellent agreement with those obtained from Monte Carlo computer simulations, which can be considered exact results for the model.49s The HNC/MSA approximation is a hybrid theory which utilizes the H N C closure for the wall-ion direct correlation function, and the bulk electrolyte MSA for the ion-ion direct correlation function. The MSA is chosen for the ion-ion interactions not only because it has an analytic solution which helps to simplify the problem, but also because it yields excellent results. A more detailed description of the HNC/MSA as applied to this type of system may be found in a review article by Carnie and Torrie.’ Calculations involving the HNC/MSA require the solution of coupled nonlinear integral ( 1 ) Carnie, S. L.; Torrie, G. M.Ado. Chrm. Phys. 1984, 56, 141-253. (2) Lozada-Cassou, M.; Saavedra-Barrera, R.;Henderson, D. J . Chcm. Phys. 1982, 77(10), 5150. (3) Lozada-Cassou, M.;Henderson, D. J . Phys. Chem. 1983, 87, 2821. (4) Torrie, G. M.; Valleau, J. P. J . Chem. Phys. 1980, 73(11), 5807. ( 5 ) Torrie, G. M.; Valleau, J. P. J. Phys. Chcm. 1982, 86, 325.
0 1992 American Chemical Society
Variational H N C for Electrical Double Layer
The Journal of Physical Chemistry, Vol. 96, No. 8,1992 3455
equations, which have no analytic solution and therefore must be solved numerically, usually by iteration. Obtaining numerical solutions is not easy, requiring what one author described as ‘heroic efforts” to ensure convergence.’ The situation is further complicated for the case of a fixed surface charge density, due to the fact that the ionic density profiles must exactly balance out the charge on the surface to satisfy the condition of electroneutrality, requiring that the density profdes be rescaled after each iteration. Alternatively, the equations can be solved for a fixed potential on the wall rather than for a fixed charge density. Additionally, the iterative approach requires a high degree of mixing of input and output to ensure convergence. Other computational methods of solution have also been de~cribed,~.’ eg., the finite element method, but these too require complex numencal methods and algorithms which make them less than convenient for use in applications. Variational solutions to integral equations have been described previously for bulk electrolytes with various approximate closure relations.*-I0 The variational method has the advantage of a relatively simple computational scheme involving only numerical integration and multivariate optimization routines, which are well described by many sources.II In addition, the variational solution is expressed in terms of analytic functions, rather than as an array of points, making it much easier to use in applications. The drawback to the variational approach is that the solutions are not exact. In this paper we will describe a variational method of solution to the HNC/MSA equations for ions near a charged planar interface and present results for various conditions. The calculations are for the restricted primitive model of an electrolyte at a charged planar surface, with a fixed surface charge density Q. The surface is a smooth plane of infinite extent, such that edge effects can be ignored, and is in contact with an infinite reservoir of electrolyte solution. The ions are described as equal sized hard spheres of diameter d, with point charges embedded in their centers. The solvent is taken to be a continuum of dielectric constant e. The effects of images and the molecular nature of the solvent are neglected.
Theory The Ornstein-Zernike equation can be applied to the problem of a planar interface using the limiting procedure of Henderson et aLi2 In their approach, the charged surface is described as a species whose density p 0, while the particle diameter d -. Henderson and Blum extended this method to the electrical double layer by considering a Coulomb system.13 Using the hypemetted chain approximation for the wall-ion closure and the bulk electrolyte MSA approximation for the ion-ion direct correlation functions results in the HNC/MSA equations for a planar double layeri4 In [ l + hi(x)] =
-
-
x
>0
(1)
h,(x) = -1 x 0
(11)
Thus the troublesome electroneutrality condition is automatically satisfied by the construction of the trial solution. The accuracy of the solution obtained with this trial function will depend on the number of variational parameters used in the calculation, as will the amount of calculation time required. The only difficulties encountered with this trial function were with systems of high surface charge and low concentration where the co-ion density is near to zero at small distances from the wall. In these cases the requirement that the ion density always remain positive makes it more difficult to optimize the variational parameters in the co-ion expansion. It is possible to avoid this problem by the addition of the zeroth-order Laguerre polynomial to the expansion (with the requirement that it also be added to the counterion trial function with the same coefficient to preserve electroneutrality). We did not, however, use that approach in the results presented here to ensure that all the calculations would be done in a consistent manner. To demonstrate this method of solution, we have carried out calculations on 1:1, 2:2, 1:2, and 1:1:2 electrolytes near a charged planar surface. Concentrations ranged from 0.010 to 1.0 M, with fixed surface charge densities from -0.16 to -0.27 C m-2. The conditions chosen for these calculations are those where the
Figure 3. Density profiles for a 1:l electrolyte. The curves have the same meaning as in Figure 2.
HNC/MSA results show large deviations from the MGC results, which are the most difficult cases to treat. For all these cases, analytic solutions to the Poisson-Boltzmann + Stem layer equation exist.I5-I7 The PB + S solutions are generally given in terms of the reduced distance KX and the fixed electrostatic potential at the Stern Layer rather than the fixed surface charge density. Using the fact that the derivative of the potential with respect to distance is related to the fixed charge density, one can express the Stern layer potential in terms of the fixed charge density. For a symmetric electrolyte the result is
$(:I (F) (-) =
arcsinh
(12)
For the cases where the electrolyte is not symmetrical, the Stern layer potential may be related to the fixed surface charge density by the solution to a cubic equation, an example of which is described in ref 17. All the results presented here were obtained using a trial function of the type in eq 9, using two variational parameters for each total correlation function hi’. The functional was minimized with respect to these parameters using a quasi-Newton method algorithm Fortran subroutine from the IMSL library.18 The functional (eq 8) evaluations were carried out with a Gauss quadrature numerical integration subroutine, also from the IMSL 1ibra1y.l~ An initial estimate of zero was used for each variational parameter when carrying out the optimization. All work was done on a Vax station 3100 M38. ( 1 5 ) Overbeek. J . Th. G.Colloid Science: Elsevier: New York, 1952; p 130. (16) Andrietti, F.; Peres, A.; Pezzotta. R. Biopkys. J . 1976, 16, 1121. (17) Abraham-Shrauner, B. J . Math. B i d . 1975, 2, 333; Erratum: ibid 1977, 4 , 201 ( 1 8) Gill, P. E.; Murray, W . Minimizarion Subject ro Simple Bounds on the Variables; NPL Report NAC 72; National Physics Laboratory, England. ( I 9) Piessens, R.; deDoncker-Kapenga, E.; Uberhuber, C. W.; Kahaner, D. K. Quadpack; Springer-Verlag: New York, 1983.
The Journal of Physical Chemistry, Vol. 96, No. 8,1992 3451
Variational H N C for Electrical Double Layer
1 : l : Z electrolyte = 0.10 M
c.,
c. !IC. I = 0.20 2 2 electrolyte
0.50 M -0.16 C m 2
‘t Flgwe 4. Electrostatic potential profiles for a 2:2 electrolyte. The curves have the same meaning as in Figure 2.
1:2 electrolyte c + =~ 0.05 M
50.8
-0.18 C m . *
g(x) 15.0
5.0
2.0
I .o X(A,
F i i 5. Density profiles for a 1:2 electrolyte. The curves have the same meaning as in Figure 2.
\
20.0
30.0
40.0
50.0
(A,
Figure 7. Density profiles for a l:1:2electrolyte. The curves have the same meaning as in Figure 2.
with MGC theory. Figure 5 demonstrates a case of fairly low concentration and charge density, for a 1:2 electrolyte, which still exhibits marked differences from MGC results. The effect of adding divalent counterions is shown in Figures 6 and 7 where it can be seen that only a small concentration of divalent counterions leads to significant differences between the two approaches. The agreement of the variational solutions with numerically obtained HNC/MSA solutions and Monte Carlo results is good. Note that the ion density profiles show the nonmonotonic behavior which cannot be found in MGC results. To check the accuracy of the method the potential at the plane of closest approach, also called the zeta potential, was calculated for the conditions shown in Figures 2 and 4 and compared with those found from the numerically obtained results of Henderson and Lozada-Cassou.20 The results agreed to within less than 2% for the monovalent electrolyte and to within 4% for the divalent case. These comparisons were done on systems where the MGC results differed by 40 and 170%, respectively.
Discussion I:l:2 electrolyte c.] = 0.10 M c.,/c. I = 0.05
\
-0.15 C m.’
10.0
: x 10.0
20.0
30.0
40.0
: x
50.0
(A)
Figure 6. Density profiles for a l:l:2 electrolyte. The curves have the same meaning as in Figure 2.
Results of our variational solutions are shown in Figures 2-7. The ion density profiles are plotted on a sinh-’ scale to allow both the co-ion and counterion to be shown in the same graph. The counterion distribution for very small distances from the wall was not shown for a few cases where the HNC/MSA results coincided with the MGC results. Figures 2 and 3 show reduced density profdes for 1:l electrolytes. The concentrated electrolyte in Figure 2 shows the significant deviations from the Poisson-Boltzmann approach. The conditions shown in Figure 3 are for a more dilute case which could have applications in colloidal physics. The electrostatic potential is plotted vs distance for a 0.50 M 2:2 electrolyte in Figure 4, showing that the potential can change sign in solutions with divalent ions, a result which cannot be predicted
The variational method is a useful tool in the solution of the integral equation theories of statistical mechanics. It allows for approximate solutions to be obtained in a straightforward manner without extensive numerical calculations. Additionally, it allows for flexibility in that the solutions can be obtained to any desired level of accuracy by the choice of the number and type of variational parameters. For example, another trial function we investigated, which is not shown here, can be constructed from the solution to the Poisson-Boltzmann equation by using the Debye-Huckel screening parameter, K , as a variational parameter. Calculations with this simple one-parameter trial function can be carried out quickly and easily and may be sufficient in many applications. Another advantage of these type of solutions is that tables of variational parameters can be constructed allowing for easier interpolation between different conditions than can be obtained from numerical solutions. A notable feature of this method is the ease with which a third type of ion was added to the calculation (see Figures 6 and 7) requiring the addition of less than six lines to the computer code written for the two-component case. To our knowledge these are the first calculations in the HNC/MSA for a mixture of monovalent and divalent counterions, which is a system of interest in many biological and colloidal systems. Acknowledgment. This work has been supported by the National Science Foundation under Grant NSF EAR 8910530. (20) Lozada-Cassou, Marcelo. Private communication.
