- Email: [email protected]
electrolyte interface
39
.J. Electroanal. Gem., 318 (1991) 39-52 Elsevier Sequoia S.A., Lausanne
JEC 01738
A Monte Carlo study of intermolecular electrostatic interactions at the electrode/ electrolyte interface Stanislaw Lamperski Department of Physical Chemistry, 60- 780 Poznan (Poland)
Faculty of Chemistry, Adam
Mickiewicz
University,
Grunwaldzka
6,
(Received 20 February 1991; in revised form 9 July 1991)
Abstract
A standard Metropolis Monte Carlo method has been adopted to simulate properties of solvent molecules occurring at the inner part of the charged electrode/electrolyte solution interface: the mean dipole moment, the differential dipole polarizability and the mean energy of intermolecular interactions. Molecules have been distributed over the electrode surface into a hexagonal or square lattice. Considerations have been restricted to the orientation effect. In the calculations, the energy of intermolecular interactions of the dipole-dipole type and interactions between the electrode surface charge and the molecular dipole moment have been taken into account. It has been observed that the increase in the energy of the intermolecular interactions caused by the increase in the permanent dipole moment of the solvent molecules, CL,leads to: (i) a decrease in the reduced mean dipole moment, (k)/p, at a given electrode surface charge, which is contrary to the results obtained from the Langevin function, (ii) an increase in the width of the dipole polarizability-surface charge curves (in some cases the transformation of a polarizability maximum into a plateau with small maxima on both sites takes place), and (iii) an almost steady dipole polarizability at the uncharged electrode. The mean energy of the intermolecular interactions is negative for small electrode charges and positive for higher values.
The energy of intermolecular electrostatic interactions between two polar solvent molecules is usually comparable to the energy of interactions between an electric dipole moment of the solvent molecule present at the inner part of the electrode/electrolyte solution interface and the electric field arising from the electrode charge. Therefore, a proper description of the intermolecular electrostatic interactions is of importance in the theoretical analysis of the interface properties. In molecular models, the intermolecular interactions of the dipole-dipole type are usually calculated from the mean molecular field approximation [l-5]. The mean molecular field has become popular because it helps to simplify the description of these interactions. However, this method allows only an approxi0022-0728/92/$05.00
0 1992 - Elsevier Sequoia S.A. All rights reserved
40
mate description of the intermolecular interactions to be obtained. Some limitations of the method have been pointed out by Marshall and Conway [6]. Here, let us note that by using the mean molecular field we cannot describe the influence of the orientation of one molecule on the orientation of its neighbours (short range correlation effect). Because of limitations of the mean molecular field, when studying the influence of the intermolecular dipole-dipole interactions on the properties of the interface, we decided to apply another treatment that involves the precisely calculated value of the intermolecular interaction energy. A computer simulation seems to be a proper method to achieve this effect. From the two main computer simulation techniques, Monte Carlo and molecular dynamics, we have chosen the Monte Carlo method as appropriate for the analysis of equilibrium states. The molecular dynamics also require the definition of additional model parameters like the molecular mass or moments of inertia. The Monte Carlo method has been used several times to study the properties of the electrode/electrolyte interface. Parsons and Reeves [7], Schmickler [8] and Aloisi et al. [91 investigated the structure of the interface at the uncharged electrode surface, whereas Guidelli et al. [lo] used the Monte Carlo method to evaluate their results obtained from a quasi-chemical approximation. Brasseur and Hurwitz used their own original simulation method to study the mercury/ propylene carbonate interface [ll]. A review of some other Monte Carlo and molecular dynamics simulations of water molecules near polar and non-polar surfaces and of an electrolyte near a non-polar surface can be found in a paper by Carnie [12]. Guidelli et al. showed that the electrical properties of the inner part of the electrode/electrolyte interface calculated for a square and hexagonal lattice of point dipoles with a continuous distribution of orientations are similar when the lattice constant is the same [lo]. Nevertheless, we have decided to analyse both the hexagonal and square lattice models. We assumed that our system consists of N2 solvent molecules distributed over the electrode surface into a regular rhombus or square lattice with dimension N X N and there are no vacant sites. Boundary conditions were assured in two ways: (i) by a periodic orientation of the dipoles for the short-range interactions, and (ii> by the mean molecular field for the long-range interactions. As we did not intend to focus on any special kind of molecule, we assumed the solvent molecules to be spherical. For the same reasons, we reduced the intermolecular interactions to interactions of the dipole-dipole type. For such an ideal system, it is convenient to suppose that the molecules do not exhibit any translation motions, but rotate only. When the molecular dipole moments are treated as point dipoles, the energy of the interactions between dipoles can be calculated from the well known equation [13]:
where & and jij are the vectors of the permanent electric dipole moment of the ith and jth molecules, cj is a vector connecting the centres of these molecules,
41
and E, is the relative electric permittivity of the inner layer if the orientation of the dipoles is fixed. We confined the application of eqn. (1) to the short-range interactions, i.e., between a given molecule and its nearest, and next nearest, neighbours localized in two rings surrounding the first (central) molecule. These rings are composed of 18 solvent molecules for the hexagonal and of 24 molecules for the square lattice. The influence of the orientation of the molecules occurring in the third ring on the energy of the dipole-dipole interactions with the central molecule is, on average, about 30 times weaker than the corresponding influence of the molecules located in the first ring. Therefore, we assumed that the long-range interactions can be described on the basis of the mean molecular field approximation [l-4]. According to this theory, the energy of the dipole can be calculated from the following equation: Q mfa= c/_L(jL) cos
~i/(47TE,E,d3)
(2)
where 4i is the angle between the vectors of the dipole moment of the ith molecule and the external electric field, (p > is the mean value of the projection of the dipole moment on the normal to the electrode surface, d is the molecular diameter and c is the effective coordination number. Topping [141 found that c is equal to 11.0342 for the hexagonal lattice and to 9.0336 for the square one. The contribution of the molecules occurring in the first two rings (short-range interactions) to the coordination number is equal to 7.9047 and 6.8065 for the hexagonal and square lattice, respectively. Thus, if we want to describe the long-range dipole interactions using eqn. (2), the coordination number, c, should be equal to 3.1295 for the hexagonal and 2.2271 for the square lattice. Furthermore, we assumed that the mean dipole moment, (p), occurring in eqn. (2) represents the average of the projections of N2 electric dipole moments on the normal to the electrode surface calculated for a given state of the molecules considered:
(3) Thus, the value of (CL) changes during the simulation process. The energy of the dipole interacting with the electric field that arises from the electrode surface charge, (T, is given by the formula: q.ud = --/_&La cos &/E&
(4)
In the present considerations we have taken no account of the interactions with dipole images in the electrode and electrolyte phases, non-electrostatic (residual) interactions and interactions with and between the induced dipole moments. To carry out the calculations, we adopted the standard Metropolis Monte Carlo algorithm 1151.At the beginning of the calculations the initial orientations for each of the N2 solvent molecules were defined by two angles: 4: and 13; (0’ is the azimuthal angle around the normal to the electrode surface; the superscript 1
42
refers to the first orientation
of the ith molecule),
cos 4: = 2( rnd - 0.5)
(5)
t$ =29r rnd
(6)
Here, md is a random number of value between 0 and 1. Succeeding orientations of the molecules were determined equations:
from the following
cos &F” = 2( rnd - 0.5)
(7)
0;”
(8)
= 0: + a,(rnd - 0.5)
Index n determines the nth orientation of the molecule. The calculations were executed exactly according to the Metropolis procedure [15] to obtain the mean value of the projection of the dipole moment on the normal to the electrode surface, (CL), the mean energy of the intermolecular dipole-dipole interactions, (U)dd, (related to one molecule), and the total mean energy of the molecule, (U). Before executing the calculations, the following quantities had to be evaluated: the value of the parameter a,, the number of configurations and of molecules required in the simulation process. To evaluate a,, we executed calculations for different values of this parameter. Some of the results are presented in Fig. 1. The outcomes are shown for a, = r (0) and a, = 27r ( n> (lines A). We see that for the higher values of a, stability of the results is obtained faster (a smaller number of configurations is required). The value of 23r for a, seems to be sensible to use in further calculations. The analysis of the number of configurations required to obtain credible results is also presented in Fig. 1. It appears from these plots that, starting from about 50000 configurations, the values of the tested quantities, i.e. (p)/p and (U), are nearly constant, independent of the electrode charge. Thus, we have taken that a number of configurations equal to 100000 should be sufficiently high to obtain sufficiently precise results. Next, we tested the influence of the number of molecules on the stability of the results. As follows from Fig. 2, both the mean electric dipole moment and the mean energy are nearly constant in the full range of the numbers of molecules analysed. We therefore decided to take into consideration 100 molecules in the simulation process. In compliance with eqn. (11, the energy of the intermolecular interactions depends on the permanent dipole moment of the solvent molecules, p, on the background electric permittivity of the inner layer, E,, and on the distance between the dipoles, rij, the value of which is connected with the diameter of the solvent molecules, d. From these parameters we have chosen the electric dipole moment, EL,as the parameter employed to model the energy of the intermolecular interactions. In the calculations we used the following values for p/C m: 5 x 10P3’ (in all figures the results for this value are denoted by (a)), 6 x 10m3’ (b), 7 x 10P3’ (c),
43
/kT 1
1
A --I C - -2 B
0.8
0.8
B --a 0.4
A-4 m.. ...... --5 ..........=...................................................... C
0.2
01
I
200 500
I
I
I
lE3 2E3 5E3
Number
I
I
I
lE4 2E4 5E4
I
1
I
lE5 2E5 6E6
I
lE8
‘-8
of configurations
-
“.*.. /kT
Fig. 1. The influence of the number of configurations on the stability of the results. The following values were used in the calculations: p = 6X lo- 3o C m, N2 = 100 molecules, d = 0.4nm, E, = 2, T = 298.15 K and a, = r (0) or 277 ( W). Calculations were executed for the hexagonal lattice at three different surface electrode charges: 0.030 (A), 0.062 (B) and 0.094 C mm2 (C). [( -1 (PcL)/PcL; (. ‘. ‘. .) (U)/kT.]
8 x 10d3’ (d), 9 X 10h3’ (e), 10 X 10P3’ (f), 11 X 10P3’ (g) and 12 X 10V3’ (h); also d = 0.4 nm, E, = 2 and T = 298.15 K. Calculations were executed at 40 different values of the electrode surface charge. The electrical properties of the electrode/electrolyte interface usually serve as a good starting point for further analysis. They pertain to the mean electric dipole moment. Therefore, we first present the dependences (p)/p =f(u> obtained from the Monte Carlo simulation. The results for hexagonal and square lattices are shown in Figs. 3 and 4, respectively. We see that in the predominant range of electrode charges an increase in p at u = const leads to a decrease in (p)/p. This outcome is contrary to the results following from the Langevin function [16]. As one might expect, the intermolecular interactions are responsible for this behaviour. This problem is discussed in the later part of this paper. Because of some fluctuations of the results presented in Figs. 3 and 4, we could not use them immediately for further calculations and that is why we decided to approximate them by an appropriate function. For perfect systems, like gases under low pressure, the dependence of (mu> on the electric field strength, E, is described by the Langevin function [16]: L(x)
= (,u)/p
= coth( x) - l/x
(9)
44 o 37
/kT
_1 15
- -1.16
0.36 - -1.17 _,...............-.' ~._........................-....,,,~ 0.35 -.......~~.~~".............~._.~....,....
0.34'
0.331 25
- -1.19 I 36
I 49
I 64
I 31
I 100
1 121
'-1.2 144
Number of molecules “‘C’.
0.8
0.6
0.4 -
0.2 a /
0 If 0
I
,
I
,
I
I
Cm
-2
I
0.02 0.04 0.06 0.08 0.1 0.12 0.14 Fig. 3. The dependence of the mean dipole moment on the electrode charge density obtained from the Monte Carlo method for the following values of the permanent dipole moment: 10e3’ p/C m: (a) 5, (b) 6, (c) 7, (d) 8, (e) 9, (f) 10, (g) 11 and (h) 12 (hexagonal lattice, d = 0.4 nm, E,= 2, T = 298.15 K, a, = 27r, number of configurations = 100000, N2 = 100).
From adip, it is easy to calculate the differential capacity of the inner layer. However, the values of the parameters p and d used in our calculations favour the appearance of the so-called Cooper-Harrison catastrophe [17] (negative capacity) in some regions of the electrode charge. Using values for the parameters Al.and d that allowed us to avoid the catastrophe, we could reduce the energy of the intermolecular interactions, and the influence of these interactions on the properties of the interface became hardly perceptible. Therefore, we restricted our further analyses to the dipolar polarizability. The dependences (Y+-(T calculated for different values of p are presented in Fig. 5 (hexagonal lattice) and in Fig. 6 (square lattice). A large number of curves exhibit a single maximum located at (T= 0. The maximum reflects the effect of reorientation of the solvent molecules. Increasing the permanent dipole moment leads to an increase in width of the curves. In some cases (curves e-h in Fig. 6) the maximum is transformed into a wide plateau with small maxima on both sites. At first, we analysed the behaviour of Q, at v = 0. As is known, the dipolar polarizability of perfect systems at E = 0 is proportional to the dipole moment squared: adip = p2/3kT
(14)
46
-2 u/Cm
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Fig. 4.The dependence of the mean dipole moment on the electrode charge density obtained from the Monte Carlo method for the square lattice; values and symbols as in Fig. 3.
~,*10
4o /Fm2
10
8
6
4
2
0 -0.15
-0.1
-0.05
0
0.05
0.1
0.15
Fig. 5. The differential dipolar polarizability as a function of the electrode surface charge for the hexagonal lattice; values and symbols as in Fig. 3.
47
a-d I
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
Fig. 6. The differential dipolar polarixability as a function of the electrode surface charge for the square lattice; values and symbols as in Fig. 3.
adip*lo
4o / Fm2
14
I
12 10 5 6 4 2
p2*lO ‘O/
(Cm)2
0 0 20 40 60 20 100 120 140 160 Fig. 7. The differential dipolar polarixability at (T= 0 as a function of the dipole moment squared calculated for perfect gases (. . . . . . ) and obtained from the Monte Carlo simulation for the hexagonal ( W) and square ( q ) lattices.
48
150 120 90 60 30 0 Fig. 8. The mean energy of intermolecular interactions at (+ = 0 as a function of the dipole moment squared obtained from the Monte Carlo simulation for the hexagonal (B) and square (0) lattices.
Therefore, in Fig. 7 the values of adip are presented as a function of p2. The dotted line corresponds to the perfect system, while the graphic characters ( W: hexagonal lattice, 0: square lattice) correspond to the results obtained from the Monte Carlo simulation. It is evident that the increase in the intermolecular interactions, caused by the increase in the dipole moment, hinders distinctly the increase in the dipolar polarizability. For higher dipole moments the polarizability reaches a nearly steady value, as is especially evident for the square lattice. It is symptomatic that the polarizability calculated for the square lattice is initially higher than that for the hexagonal one. To explain this phenomenon let us use the mean energy of the dipole-dipole interactions calculated for the uncharged electrode surface. The appropriate results are presented in Fig. 8 ( n : hexagonal lattice, 0: square lattice). For small dipole moments the interactions are stronger for molecules occupying a hexagonal lattice. This is connected with the effective coordination number, c, which is higher for the hexagonal lattice. The weaker intermolecular interactions lead to a higher dipolar polarizability. On the other hand, for higher values of the dipole moment (CL> 8 x 10W3’ C m) the intermolecular interactions are evidently stronger for molecules occupying the square lattice. Here, not the number of “effective” neighbours but the mutual orientation of molecules, leading to the formation of stable surface structures, plays a more important role. The increase in width of the polarizability curves and, in some cases, the transformation of the maximum into a wide plateau can also be related to the
49
-6’
I
I
I
I
I
I
I
0.14 0.12 0.1 0.08 0.00 0.02 0.04 Fig. 9. The dependence of the mean energy of intermolecular interactions on the electrode charge density obtained from the Monte Carlo simulation for the hexagonal lattice; values and symbols as in Fig. 3.
0
intermolecular interactions. The mechanism of this effect will become clearer if we analyse the results of the mean energy of intermolecular interactions, (U& in a wide range of electrode charges. The dependences of (U )dd on u for the hexagonal and square lattices are shown in Figs. 9 and 10, respectively. It can be seen that for small electrode charges the values of (U )dd are negative. Thus, we can conclude that in this range of electrode charges the solvent molecules form stable structures on the electrode surface. Such structures have been described in the literature before [7]. One of the most stable structures obtained during the simulation process for the square lattice is shown schematically in Fig. 11. The stability of these structures increases with an increase in the dipole moment of the solvent molecules. For high dipole moments, when the intermolecular interactions are strong, the weak external electric field arising from the small electrode charges is not able to disturb these structures completely but only to make a slight change. Therefore, the observed dipolar polarizability near CT= 0 reflects the ability of the molecular, structures to change. As the range of the electrode charge where the steady molecular structures occur increases with an increase in the dipole moment, we observe an increase in width of the polarizability curves or even a plateau. The appropriately strong electric field arising from the higher electrode surface charges is able to disturb these structures. The molecules are gradually arranged by the external electric field in such a way that their dipole moments are oriented parallel
50
-61
/kT
I
I
I
I
I
1
I
-0
0.04 0.02 0.02 0.06 0.1 0.12 0.14 Fig. 10. The dependence of the mean energy of intermolecular interactions on the electrode charge density obtained from the Monte Carlo simulation for the square lattice; values and symbols as in Fig. ?
to the vector of this field. The new arrangement of the dipoles leads to a decrease in the dipolar polarizability (the saturation effect). The mean energy of the intermolecular interactions increases and becomes positive. Some comments should be provided for the polarizability curves e-h, obtained for the square lattice. As strong intermolecular interactions are required to observe the polarizability plateau, for real systems with more or less spherical
Fig. 11. One of the more stable dipole structures for the square lattice at uncharged electrode surface (all dipoles are nearly parallel to the electrode surface).
51
molecules it is less probable that this effect will occur. However, for plane molecules with a high electric dipole moment like ethylene or propylene carbonate, the intermolecular interactions are strong. Such molecules are capable of forming especially stable molecular structures. The presence of such structures may be used to explain the electrical properties of the mercury/propylene carbonate interface, as we did earlier [lS]. Additionally, let us note that, contrary to our results, the mean energy of the intermolecular interactions calculated from the mean molecular field [l-.5] is not negative in the whole range of the electrode charges. This comparison shows the limitations of the mean field when it is used to describe intermolecular interactions. The results presented in this paper evidently show the existence of stable structures of solvent molecules present on the electrode surface at small electric charges if the intermolecular interactions of the dipole type are sufficiently strong. The presence of such structures in the liquid phase was proved experimentally earlier [193 and confirmed by other investigations [20]. Such structures, composed of three solvent molecules adsorbed on the electrode surface where each molecule can take three orientations, have been regarded recently by Fawcett [21]. However, a precise description of the intermolecular dipole-dipole interactions in molecular models of the interface requires further investigations. We hope that the results presented in this paper will be useful in estimating the value of new models. ACKNOWLEDGEMENTS
I would like to thank Professor E. Dutkiewicz for helpful discussions. This work was carried out within the framework of the KBN Grant No. 207119101. REFERENCES 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
J.O’M. Bockris, M.A.V. Devanathan and K. Miiller, Proc. R. Sot. (London), Ser. A, 274 (1963) 55. S. Levine, B.M. Bell and A.L. Smith, J. Chem. Phys., 73 (1969) 3534. I.L. Cooper and J.A. Harrison, J. Electroanal. Chem., 60 (1975) 85. W.R. Fawcett, J. Chem. Phys., 82 (1978) 1385. E. Dutkiewicz and S. Lamperski, J. Electroanal. Chem., 289 (1990) 1. S.L. Marshall and B.E. Conway, J. Electroanal. Chem., 227 (1987) 245. R. Parsons and R.M. Reeves, J. Electroanal. Chem., 123 (1981) 141. W. Schmickler, J. Electroanal. Chem., 157 (1983) 1. G. Aloisi, R. Guidelli, R.A. Jackson, SM. Clark and P. Barnes, J. Electroanal. Chem., 206 (1986) 131. R. Guidelli, G. Aloisi, M. Carla and M.R. Moncelli, J. Electroanal. Chem., 197 (1986) 143. R. Brasseur and H.D. Hunvitz, J. Electroanal. Chem., 148 (1983) 249. S.L. Camie, Ber. Bunsenges. Phys. Chem., 91 (1987) 262. For example: C.F.J. Bottcher, Theory of Electric Polarization, 2nd ed., Elsevier, Amsterdam, 1973, p. 118. J. Topping, Proc. R. Sot. (London), Ser. A, 114 (1927) 67. N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys., 21 (1953) 1087.
52 16 T.L. Hill, An Introduction to Statistical Thermodynamics, Dover Publications, New York, 1986, Ch. 17 18 19 20 21
I’.;. Cooper and J.A. Harrison, J. Electroanal. Chem., 66 (1975) 85. E. Dutkiewicz and S. Lamperski, J. Electroanal. Chem., 247 (1988) 57. A. Piekara, Proc. R. Sot. (London), Ser. A, 172 (1939) 360. M. Dutkiewicz, Chem. Phys. Lett., 121 (1985) 73. W.R. Fawcett, J. Chem. Phys., 93 (1990) 6813.
.J. Electroanal. Gem., 318 (1991) 39-52 Elsevier Sequoia S.A., Lausanne
JEC 01738
A Monte Carlo study of intermolecular electrostatic interactions at the electrode/ electrolyte interface Stanislaw Lamperski Department of Physical Chemistry, 60- 780 Poznan (Poland)
Faculty of Chemistry, Adam
Mickiewicz
University,
Grunwaldzka
6,
(Received 20 February 1991; in revised form 9 July 1991)
Abstract
A standard Metropolis Monte Carlo method has been adopted to simulate properties of solvent molecules occurring at the inner part of the charged electrode/electrolyte solution interface: the mean dipole moment, the differential dipole polarizability and the mean energy of intermolecular interactions. Molecules have been distributed over the electrode surface into a hexagonal or square lattice. Considerations have been restricted to the orientation effect. In the calculations, the energy of intermolecular interactions of the dipole-dipole type and interactions between the electrode surface charge and the molecular dipole moment have been taken into account. It has been observed that the increase in the energy of the intermolecular interactions caused by the increase in the permanent dipole moment of the solvent molecules, CL,leads to: (i) a decrease in the reduced mean dipole moment, (k)/p, at a given electrode surface charge, which is contrary to the results obtained from the Langevin function, (ii) an increase in the width of the dipole polarizability-surface charge curves (in some cases the transformation of a polarizability maximum into a plateau with small maxima on both sites takes place), and (iii) an almost steady dipole polarizability at the uncharged electrode. The mean energy of the intermolecular interactions is negative for small electrode charges and positive for higher values.
The energy of intermolecular electrostatic interactions between two polar solvent molecules is usually comparable to the energy of interactions between an electric dipole moment of the solvent molecule present at the inner part of the electrode/electrolyte solution interface and the electric field arising from the electrode charge. Therefore, a proper description of the intermolecular electrostatic interactions is of importance in the theoretical analysis of the interface properties. In molecular models, the intermolecular interactions of the dipole-dipole type are usually calculated from the mean molecular field approximation [l-5]. The mean molecular field has become popular because it helps to simplify the description of these interactions. However, this method allows only an approxi0022-0728/92/$05.00
0 1992 - Elsevier Sequoia S.A. All rights reserved
40
mate description of the intermolecular interactions to be obtained. Some limitations of the method have been pointed out by Marshall and Conway [6]. Here, let us note that by using the mean molecular field we cannot describe the influence of the orientation of one molecule on the orientation of its neighbours (short range correlation effect). Because of limitations of the mean molecular field, when studying the influence of the intermolecular dipole-dipole interactions on the properties of the interface, we decided to apply another treatment that involves the precisely calculated value of the intermolecular interaction energy. A computer simulation seems to be a proper method to achieve this effect. From the two main computer simulation techniques, Monte Carlo and molecular dynamics, we have chosen the Monte Carlo method as appropriate for the analysis of equilibrium states. The molecular dynamics also require the definition of additional model parameters like the molecular mass or moments of inertia. The Monte Carlo method has been used several times to study the properties of the electrode/electrolyte interface. Parsons and Reeves [7], Schmickler [8] and Aloisi et al. [91 investigated the structure of the interface at the uncharged electrode surface, whereas Guidelli et al. [lo] used the Monte Carlo method to evaluate their results obtained from a quasi-chemical approximation. Brasseur and Hurwitz used their own original simulation method to study the mercury/ propylene carbonate interface [ll]. A review of some other Monte Carlo and molecular dynamics simulations of water molecules near polar and non-polar surfaces and of an electrolyte near a non-polar surface can be found in a paper by Carnie [12]. Guidelli et al. showed that the electrical properties of the inner part of the electrode/electrolyte interface calculated for a square and hexagonal lattice of point dipoles with a continuous distribution of orientations are similar when the lattice constant is the same [lo]. Nevertheless, we have decided to analyse both the hexagonal and square lattice models. We assumed that our system consists of N2 solvent molecules distributed over the electrode surface into a regular rhombus or square lattice with dimension N X N and there are no vacant sites. Boundary conditions were assured in two ways: (i) by a periodic orientation of the dipoles for the short-range interactions, and (ii> by the mean molecular field for the long-range interactions. As we did not intend to focus on any special kind of molecule, we assumed the solvent molecules to be spherical. For the same reasons, we reduced the intermolecular interactions to interactions of the dipole-dipole type. For such an ideal system, it is convenient to suppose that the molecules do not exhibit any translation motions, but rotate only. When the molecular dipole moments are treated as point dipoles, the energy of the interactions between dipoles can be calculated from the well known equation [13]:
where & and jij are the vectors of the permanent electric dipole moment of the ith and jth molecules, cj is a vector connecting the centres of these molecules,
41
and E, is the relative electric permittivity of the inner layer if the orientation of the dipoles is fixed. We confined the application of eqn. (1) to the short-range interactions, i.e., between a given molecule and its nearest, and next nearest, neighbours localized in two rings surrounding the first (central) molecule. These rings are composed of 18 solvent molecules for the hexagonal and of 24 molecules for the square lattice. The influence of the orientation of the molecules occurring in the third ring on the energy of the dipole-dipole interactions with the central molecule is, on average, about 30 times weaker than the corresponding influence of the molecules located in the first ring. Therefore, we assumed that the long-range interactions can be described on the basis of the mean molecular field approximation [l-4]. According to this theory, the energy of the dipole can be calculated from the following equation: Q mfa= c/_L(jL) cos
~i/(47TE,E,d3)
(2)
where 4i is the angle between the vectors of the dipole moment of the ith molecule and the external electric field, (p > is the mean value of the projection of the dipole moment on the normal to the electrode surface, d is the molecular diameter and c is the effective coordination number. Topping [141 found that c is equal to 11.0342 for the hexagonal lattice and to 9.0336 for the square one. The contribution of the molecules occurring in the first two rings (short-range interactions) to the coordination number is equal to 7.9047 and 6.8065 for the hexagonal and square lattice, respectively. Thus, if we want to describe the long-range dipole interactions using eqn. (2), the coordination number, c, should be equal to 3.1295 for the hexagonal and 2.2271 for the square lattice. Furthermore, we assumed that the mean dipole moment, (p), occurring in eqn. (2) represents the average of the projections of N2 electric dipole moments on the normal to the electrode surface calculated for a given state of the molecules considered:
(3) Thus, the value of (CL) changes during the simulation process. The energy of the dipole interacting with the electric field that arises from the electrode surface charge, (T, is given by the formula: q.ud = --/_&La cos &/E&
(4)
In the present considerations we have taken no account of the interactions with dipole images in the electrode and electrolyte phases, non-electrostatic (residual) interactions and interactions with and between the induced dipole moments. To carry out the calculations, we adopted the standard Metropolis Monte Carlo algorithm 1151.At the beginning of the calculations the initial orientations for each of the N2 solvent molecules were defined by two angles: 4: and 13; (0’ is the azimuthal angle around the normal to the electrode surface; the superscript 1
42
refers to the first orientation
of the ith molecule),
cos 4: = 2( rnd - 0.5)
(5)
t$ =29r rnd
(6)
Here, md is a random number of value between 0 and 1. Succeeding orientations of the molecules were determined equations:
from the following
cos &F” = 2( rnd - 0.5)
(7)
0;”
(8)
= 0: + a,(rnd - 0.5)
Index n determines the nth orientation of the molecule. The calculations were executed exactly according to the Metropolis procedure [15] to obtain the mean value of the projection of the dipole moment on the normal to the electrode surface, (CL), the mean energy of the intermolecular dipole-dipole interactions, (U)dd, (related to one molecule), and the total mean energy of the molecule, (U). Before executing the calculations, the following quantities had to be evaluated: the value of the parameter a,, the number of configurations and of molecules required in the simulation process. To evaluate a,, we executed calculations for different values of this parameter. Some of the results are presented in Fig. 1. The outcomes are shown for a, = r (0) and a, = 27r ( n> (lines A). We see that for the higher values of a, stability of the results is obtained faster (a smaller number of configurations is required). The value of 23r for a, seems to be sensible to use in further calculations. The analysis of the number of configurations required to obtain credible results is also presented in Fig. 1. It appears from these plots that, starting from about 50000 configurations, the values of the tested quantities, i.e. (p)/p and (U), are nearly constant, independent of the electrode charge. Thus, we have taken that a number of configurations equal to 100000 should be sufficiently high to obtain sufficiently precise results. Next, we tested the influence of the number of molecules on the stability of the results. As follows from Fig. 2, both the mean electric dipole moment and the mean energy are nearly constant in the full range of the numbers of molecules analysed. We therefore decided to take into consideration 100 molecules in the simulation process. In compliance with eqn. (11, the energy of the intermolecular interactions depends on the permanent dipole moment of the solvent molecules, p, on the background electric permittivity of the inner layer, E,, and on the distance between the dipoles, rij, the value of which is connected with the diameter of the solvent molecules, d. From these parameters we have chosen the electric dipole moment, EL,as the parameter employed to model the energy of the intermolecular interactions. In the calculations we used the following values for p/C m: 5 x 10P3’ (in all figures the results for this value are denoted by (a)), 6 x 10m3’ (b), 7 x 10P3’ (c),
43
/kT 1
1
A --I C - -2 B
0.8
0.8
B --a 0.4
A-4 m.. ...... --5 ..........=...................................................... C
0.2
01
I
200 500
I
I
I
lE3 2E3 5E3
Number
I
I
I
lE4 2E4 5E4
I
1
I
lE5 2E5 6E6
I
lE8
‘-8
of configurations
-
“.*.. /kT
Fig. 1. The influence of the number of configurations on the stability of the results. The following values were used in the calculations: p = 6X lo- 3o C m, N2 = 100 molecules, d = 0.4nm, E, = 2, T = 298.15 K and a, = r (0) or 277 ( W). Calculations were executed for the hexagonal lattice at three different surface electrode charges: 0.030 (A), 0.062 (B) and 0.094 C mm2 (C). [( -1 (PcL)/PcL; (. ‘. ‘. .) (U)/kT.]
8 x 10d3’ (d), 9 X 10h3’ (e), 10 X 10P3’ (f), 11 X 10P3’ (g) and 12 X 10V3’ (h); also d = 0.4 nm, E, = 2 and T = 298.15 K. Calculations were executed at 40 different values of the electrode surface charge. The electrical properties of the electrode/electrolyte interface usually serve as a good starting point for further analysis. They pertain to the mean electric dipole moment. Therefore, we first present the dependences (p)/p =f(u> obtained from the Monte Carlo simulation. The results for hexagonal and square lattices are shown in Figs. 3 and 4, respectively. We see that in the predominant range of electrode charges an increase in p at u = const leads to a decrease in (p)/p. This outcome is contrary to the results following from the Langevin function [16]. As one might expect, the intermolecular interactions are responsible for this behaviour. This problem is discussed in the later part of this paper. Because of some fluctuations of the results presented in Figs. 3 and 4, we could not use them immediately for further calculations and that is why we decided to approximate them by an appropriate function. For perfect systems, like gases under low pressure, the dependence of (mu> on the electric field strength, E, is described by the Langevin function [16]: L(x)
= (,u)/p
= coth( x) - l/x
(9)
44 o 37
/kT
_1 15
- -1.16
0.36 - -1.17 _,...............-.' ~._........................-....,,,~ 0.35 -.......~~.~~".............~._.~....,....
0.34'
0.331 25
- -1.19 I 36
I 49
I 64
I 31
I 100
1 121
'-1.2 144
Number of molecules “‘C’.
/p -
/kT
Fig. 2. The influence of the number of molecules on the stability of the results. The following values were used in the calculations: p = 6 X lo- 3o C m , d = 0.4 nm, E, = 2, T = 298.15 K, n, = 27~, number of configurations = 200000. Calculations were executed at c = 0.030 C m-*. [(. . . .) (p)/k; (1
(W/%T.l where x =pE/kT
(10)
To approximate our results, we used the following formula: (P&,=P[L(u~)
+bo
exp( -ca2)]
(11) where the second term and the parameter a in the Langevin function reflect the perturbation caused by the intermolecular interactions. The values of the coefficients a, b and c were obtained from the non-linear least-squares fitting method. Some problems with fitting the approximation function (eqn. 11) to the Monte Carlo results appeared in the case of high electric dipole moments of molecules occupying the square lattice (Fig. 4, curves e-h). Here, we used a polynomial approximation. Now the analytical form for the dependence (p) = f(a) was known, we could calculate the differential dipolar polarizability, Lydip,of the molecules present in the inner layer. This value was obtained by differentiating the approximating function of the mean dipole moment with respect to the external electric field (dielectric displacement), E: adip = a( p )app/aE
(12)
where E = a/q, as the field E can be regarded as an intensive thermodynamic
parameter
(13) [16].
45
1
0.8
0.6
0.4 -
0.2 a /
0 If 0
I
,
I
,
I
I
Cm
-2
I
0.02 0.04 0.06 0.08 0.1 0.12 0.14 Fig. 3. The dependence of the mean dipole moment on the electrode charge density obtained from the Monte Carlo method for the following values of the permanent dipole moment: 10e3’ p/C m: (a) 5, (b) 6, (c) 7, (d) 8, (e) 9, (f) 10, (g) 11 and (h) 12 (hexagonal lattice, d = 0.4 nm, E,= 2, T = 298.15 K, a, = 27r, number of configurations = 100000, N2 = 100).
From adip, it is easy to calculate the differential capacity of the inner layer. However, the values of the parameters p and d used in our calculations favour the appearance of the so-called Cooper-Harrison catastrophe [17] (negative capacity) in some regions of the electrode charge. Using values for the parameters Al.and d that allowed us to avoid the catastrophe, we could reduce the energy of the intermolecular interactions, and the influence of these interactions on the properties of the interface became hardly perceptible. Therefore, we restricted our further analyses to the dipolar polarizability. The dependences (Y+-(T calculated for different values of p are presented in Fig. 5 (hexagonal lattice) and in Fig. 6 (square lattice). A large number of curves exhibit a single maximum located at (T= 0. The maximum reflects the effect of reorientation of the solvent molecules. Increasing the permanent dipole moment leads to an increase in width of the curves. In some cases (curves e-h in Fig. 6) the maximum is transformed into a wide plateau with small maxima on both sites. At first, we analysed the behaviour of Q, at v = 0. As is known, the dipolar polarizability of perfect systems at E = 0 is proportional to the dipole moment squared: adip = p2/3kT
(14)
46
-2 u/Cm
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Fig. 4.The dependence of the mean dipole moment on the electrode charge density obtained from the Monte Carlo method for the square lattice; values and symbols as in Fig. 3.
~,*10
4o /Fm2
10
8
6
4
2
0 -0.15
-0.1
-0.05
0
0.05
0.1
0.15
Fig. 5. The differential dipolar polarizability as a function of the electrode surface charge for the hexagonal lattice; values and symbols as in Fig. 3.
47
a-d I
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
Fig. 6. The differential dipolar polarixability as a function of the electrode surface charge for the square lattice; values and symbols as in Fig. 3.
adip*lo
4o / Fm2
14
I
12 10 5 6 4 2
p2*lO ‘O/
(Cm)2
0 0 20 40 60 20 100 120 140 160 Fig. 7. The differential dipolar polarixability at (T= 0 as a function of the dipole moment squared calculated for perfect gases (. . . . . . ) and obtained from the Monte Carlo simulation for the hexagonal ( W) and square ( q ) lattices.
48
150 120 90 60 30 0 Fig. 8. The mean energy of intermolecular interactions at (+ = 0 as a function of the dipole moment squared obtained from the Monte Carlo simulation for the hexagonal (B) and square (0) lattices.
Therefore, in Fig. 7 the values of adip are presented as a function of p2. The dotted line corresponds to the perfect system, while the graphic characters ( W: hexagonal lattice, 0: square lattice) correspond to the results obtained from the Monte Carlo simulation. It is evident that the increase in the intermolecular interactions, caused by the increase in the dipole moment, hinders distinctly the increase in the dipolar polarizability. For higher dipole moments the polarizability reaches a nearly steady value, as is especially evident for the square lattice. It is symptomatic that the polarizability calculated for the square lattice is initially higher than that for the hexagonal one. To explain this phenomenon let us use the mean energy of the dipole-dipole interactions calculated for the uncharged electrode surface. The appropriate results are presented in Fig. 8 ( n : hexagonal lattice, 0: square lattice). For small dipole moments the interactions are stronger for molecules occupying a hexagonal lattice. This is connected with the effective coordination number, c, which is higher for the hexagonal lattice. The weaker intermolecular interactions lead to a higher dipolar polarizability. On the other hand, for higher values of the dipole moment (CL> 8 x 10W3’ C m) the intermolecular interactions are evidently stronger for molecules occupying the square lattice. Here, not the number of “effective” neighbours but the mutual orientation of molecules, leading to the formation of stable surface structures, plays a more important role. The increase in width of the polarizability curves and, in some cases, the transformation of the maximum into a wide plateau can also be related to the
49
-6’
I
I
I
I
I
I
I
0.14 0.12 0.1 0.08 0.00 0.02 0.04 Fig. 9. The dependence of the mean energy of intermolecular interactions on the electrode charge density obtained from the Monte Carlo simulation for the hexagonal lattice; values and symbols as in Fig. 3.
0
intermolecular interactions. The mechanism of this effect will become clearer if we analyse the results of the mean energy of intermolecular interactions, (U& in a wide range of electrode charges. The dependences of (U )dd on u for the hexagonal and square lattices are shown in Figs. 9 and 10, respectively. It can be seen that for small electrode charges the values of (U )dd are negative. Thus, we can conclude that in this range of electrode charges the solvent molecules form stable structures on the electrode surface. Such structures have been described in the literature before [7]. One of the most stable structures obtained during the simulation process for the square lattice is shown schematically in Fig. 11. The stability of these structures increases with an increase in the dipole moment of the solvent molecules. For high dipole moments, when the intermolecular interactions are strong, the weak external electric field arising from the small electrode charges is not able to disturb these structures completely but only to make a slight change. Therefore, the observed dipolar polarizability near CT= 0 reflects the ability of the molecular, structures to change. As the range of the electrode charge where the steady molecular structures occur increases with an increase in the dipole moment, we observe an increase in width of the polarizability curves or even a plateau. The appropriately strong electric field arising from the higher electrode surface charges is able to disturb these structures. The molecules are gradually arranged by the external electric field in such a way that their dipole moments are oriented parallel
50
-61
/kT
I
I
I
I
I
1
I
-0
0.04 0.02 0.02 0.06 0.1 0.12 0.14 Fig. 10. The dependence of the mean energy of intermolecular interactions on the electrode charge density obtained from the Monte Carlo simulation for the square lattice; values and symbols as in Fig. ?
to the vector of this field. The new arrangement of the dipoles leads to a decrease in the dipolar polarizability (the saturation effect). The mean energy of the intermolecular interactions increases and becomes positive. Some comments should be provided for the polarizability curves e-h, obtained for the square lattice. As strong intermolecular interactions are required to observe the polarizability plateau, for real systems with more or less spherical
Fig. 11. One of the more stable dipole structures for the square lattice at uncharged electrode surface (all dipoles are nearly parallel to the electrode surface).
51
molecules it is less probable that this effect will occur. However, for plane molecules with a high electric dipole moment like ethylene or propylene carbonate, the intermolecular interactions are strong. Such molecules are capable of forming especially stable molecular structures. The presence of such structures may be used to explain the electrical properties of the mercury/propylene carbonate interface, as we did earlier [lS]. Additionally, let us note that, contrary to our results, the mean energy of the intermolecular interactions calculated from the mean molecular field [l-.5] is not negative in the whole range of the electrode charges. This comparison shows the limitations of the mean field when it is used to describe intermolecular interactions. The results presented in this paper evidently show the existence of stable structures of solvent molecules present on the electrode surface at small electric charges if the intermolecular interactions of the dipole type are sufficiently strong. The presence of such structures in the liquid phase was proved experimentally earlier [193 and confirmed by other investigations [20]. Such structures, composed of three solvent molecules adsorbed on the electrode surface where each molecule can take three orientations, have been regarded recently by Fawcett [21]. However, a precise description of the intermolecular dipole-dipole interactions in molecular models of the interface requires further investigations. We hope that the results presented in this paper will be useful in estimating the value of new models. ACKNOWLEDGEMENTS
I would like to thank Professor E. Dutkiewicz for helpful discussions. This work was carried out within the framework of the KBN Grant No. 207119101. REFERENCES 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
J.O’M. Bockris, M.A.V. Devanathan and K. Miiller, Proc. R. Sot. (London), Ser. A, 274 (1963) 55. S. Levine, B.M. Bell and A.L. Smith, J. Chem. Phys., 73 (1969) 3534. I.L. Cooper and J.A. Harrison, J. Electroanal. Chem., 60 (1975) 85. W.R. Fawcett, J. Chem. Phys., 82 (1978) 1385. E. Dutkiewicz and S. Lamperski, J. Electroanal. Chem., 289 (1990) 1. S.L. Marshall and B.E. Conway, J. Electroanal. Chem., 227 (1987) 245. R. Parsons and R.M. Reeves, J. Electroanal. Chem., 123 (1981) 141. W. Schmickler, J. Electroanal. Chem., 157 (1983) 1. G. Aloisi, R. Guidelli, R.A. Jackson, SM. Clark and P. Barnes, J. Electroanal. Chem., 206 (1986) 131. R. Guidelli, G. Aloisi, M. Carla and M.R. Moncelli, J. Electroanal. Chem., 197 (1986) 143. R. Brasseur and H.D. Hunvitz, J. Electroanal. Chem., 148 (1983) 249. S.L. Camie, Ber. Bunsenges. Phys. Chem., 91 (1987) 262. For example: C.F.J. Bottcher, Theory of Electric Polarization, 2nd ed., Elsevier, Amsterdam, 1973, p. 118. J. Topping, Proc. R. Sot. (London), Ser. A, 114 (1927) 67. N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys., 21 (1953) 1087.
52 16 T.L. Hill, An Introduction to Statistical Thermodynamics, Dover Publications, New York, 1986, Ch. 17 18 19 20 21
I’.;. Cooper and J.A. Harrison, J. Electroanal. Chem., 66 (1975) 85. E. Dutkiewicz and S. Lamperski, J. Electroanal. Chem., 247 (1988) 57. A. Piekara, Proc. R. Sot. (London), Ser. A, 172 (1939) 360. M. Dutkiewicz, Chem. Phys. Lett., 121 (1985) 73. W.R. Fawcett, J. Chem. Phys., 93 (1990) 6813.