Molecular dynamics calculations of the electrochemical properties of electrolyte systems between charged electrodes

Molecular dynamics calculations of the electrochemical properties of electrolyte systems between charged electrodes
We investigate the interfacial electrochemical properties of an aqueous electrolyte solution with discrete water molecules in slab geometry between charged atomistic electrodes. Long-range intermolecular Coulombic interactions are calculated using the particle particle particle mesh method with a modification to account for the slab geometry. Density distribution profiles and potential drops across the double layer are given for 0, 0.25, and 1 M aqueous electrolyte solutions each electrode surface charges. Results are compared qualitatively with experimental x-ray scattering findings, other computer simulation results, and traditional electrochemistry theory. The interfacial fluid structure characteristics are generally in good qualitative agreement with the conclusions obtained in some integral equation theories and in the experimental x-ray study. The potential in the simulations shows an oscillatory behavior near the electrode, which theories that do not include the molecular nature of water cannot reproduce for the given conditions. Surprisingly, the results also show that the water structure near the electrode is dominated by the charge on the electrode and is fairly insensitive to the ion concentrations. Except at large electrode charge, the potential drop across the double layer does not depend significantly upon the concentration of the ions.
Molecular dynamics calculations of the electrochemical properties of electrolyte systems between charged electrodes Paul S. Croziera) and Richard L. Rowleyb) Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602- 4100 Douglas Hendersonc) Department of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602- 5700 Received 27 June 2000; accepted 1 September 2000 ! We investigate the interfacial electrochemical properties of an aqueous electrolyte solution with discrete water molecules in slab geometry between charged atomistic electrodes. Long- range intermolecular Coulombic interactions are calculated using the particle particle particle mesh method with a modification to account for the slab geometry. Density distribution profiles and potential drops across the double layer are given for 0, 0.25, and 1 M aqueous electrolyte solutions each at 0, 6 0.1, 6 0.2, and 6 0.3 C/ m2 electrode surface charges. Results are compared qualitatively with experimental x- ray scattering findings, other computer simulation results, and traditional electrochemistry theory. The interfacial fluid structure characteristics are generally in good qualitative agreement with the conclusions obtained in some integral equation theories and in the experimental x- ray study. The potential in the simulations shows an oscillatory behavior near the electrode, which theories that do not include the molecular nature of water cannot reproduce for the given conditions. Surprisingly, the results also show that the water structure near the electrode is dominated by the charge on the electrode and is fairly insensitive to the ion concentrations. Except at large electrode charge, the potential drop across the double layer does not depend significantly upon the concentration of the ions. 2000 American Institute of Physics. @ S0021- 9606 00 ! 51344- 0 # I. INTRODUCTION The ability to model correctly electrochemical processes at interfaces is a key to our understanding of numerous very important chemical and biological processes. Especially of interest are aqueous solutions in contact with charged sur-faces, which can be found in such diverse fields as molecular biology, electrochemistry, and heterogeneous catalysis. Although unable to model quantum mechanical effects, such as electron transfer, classical molecular dynamics is the most valuable tool available for the elucidation of important interfacial processes. Such simulations provide realistic rep-resentations of interfacial electrochemistry processes, limited primarily by molecular model accuracy and computational demands. In order to save time, many interfacial electrochemistry molecular dynamics researchers neglect long- range Coulom-bic interactions while admitting the vital role that these play. 1 5 Other researchers have employed costly Ewald sum-mation approximations for long- range interactions. 6 9 Re-cently, mesh routines have been developed, 10 15 which present a more efficient way to perform these calculations. Of the mesh routines, the particle particle particle mesh P3M ! routine developed by Hockney and Eastwood10 has been shown to be the most efficient for a given accuracy level. 14,15 In this work, we modify P3M for application to slab geometry, as suggested by Pollock and Glosli16 and similar to the 3D Ewald summation modification used by Yeh and Berkowitz. 8 We show that this implementation of the slab- corrected P3M method, or P3MC as we shall call it, gives results in very good agreement with regular Ewald summation methods, yet at a fraction of the cost. Once we have established the validity of the P3MC methodology, we present simulation results obtained using this method dealing with the effect of electrode charge and ion concentration on the potential drop across the double layer of an aqueous electrolyte solution at a solid atomistic surface. We also note the similarity of our results with the experimental force measurements of Israelachvili et al. 17 and the x- ray scattering findings of Toney et al. 18 In addition, we compare our findings with the predictions of traditional in-terfacial electrochemistry theories, as well as with the find-ings from other simulation work. II. P3MC IMPLEMENTATION Many methods for approximating long- range Coulombic interactions have been devised, including the charged sheets method, 19,20 the modified reaction field, 21 24 Ewald sums, 6 9,25 27 the fast multipole method, 16,28 30 and mesh methods. 10 15 Recently, Crozier et al. 9 found that the charged sheets method is inadequate for correctly modeling strong long- range Coulombic interactions for systems with discrete solvent molecules. Ewald summations appear to be adequate, yet computationally expensive, and multipole methods are even more costly for systems of a reasonable size. 16 Mesh methods, particularly P3M, present the most flexible and reliable alternative. 14 16 Deserno and Holm give a ! Electronic mail: crozierp@ et. byu. edu b ! Electronic mail: rowley@ byu. edu c ! Electronic mail: doug@ huey. byu. edu JOURNAL OF CHEMICAL PHYSICS VOLUME 113, NUMBER 20 22 NOVEMBER 2000 0021- 9606/ 2000/ 113( 20)/ 9202/ 6/$ 17.00 9202 2000 American Institute of Physics Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp an excellent analysis of the various mesh methods, 14 along with useful recommendations for P3M implementation. 15 We have followed their recommendations and adopted their no-tation in the following. We use Fourier space, or ik, differentiation to calculate the intermolecular forces from the electrostatic energies, and a seventh- order assignment scheme ( P 5 7) to assign the charges to the mesh. Calculation of the optimal influence function, G opt , of Hockney and Eastwood10 is done at the beginning of each run and therefore presents no additional run- time overhead. The cutoff for the Lennard- Jones LJ ! interactions, and the real- space Coulombic interactions is set at rmax 5 10 , and the mesh size is set at 16 3 16 3 64 in the x, y, and z directions, respectively. With a corresponding fixed box size of 25.5 3 25.5 3 85.1232 , and the total number of ions and water molecules fixed at 1000 the opti-mum value of a is estimated to be 0.3007. We estimate the average error in our force calculations to be on the order of 10 2 5 dimensionless force units, according to the estimation method and dimensionless force units of Deserno and Holm. 14 Although not as precise as an Ewald sum, the amount of error introduced by this implementation of P3M is very reasonable. Extreme accuracy is not needed, especially when considering other possible sources of error in the MD simulation, such as the use of a randomly fluctuating ther-mostat, or a discrete time step e. g., 2.5 fs ! . In order to account for the slab geometry, we include an empty space, equal to the size of the fluid- occupied space, between repeating slabs of fluid, and a correction term to damp out interslab interactions. With the implementation of this correction, the 3D, periodically repeating simulation closely approximates the simulation of an isolated 2D slab of finite thickness, as shown by Yeh and Berkowitz. 8 To confirm the reliability of our P3MC implementation, we compare calculated forces with those calculated using the corrected 3D Ewald EW3DC ! method. 8 In particular, we calculate the force between two oppositely charged particles along with their periodic images in slab geometry. For easy comparison, we consider the same case as presented in Refs. 7 9 ( L 5 Lx 5 Ly 5 18 and Lz 5 90 ! . As seen in Figs. 1 and 2, the P3MC and EW3DC calculations are visually in-distinguishable. For Fig. 1, the difference between the two calculations was less than 0.013% in all cases, with an aver-age difference of 0.002%. Even though the deviation was up to 1% in the case of Fig. 2, percent deviation is less mean-ingful here, since the absolute value of the force approaches zero. The absolute difference between the two calculations on the scale of Fig. 2 was always less than 10 2 4, which represents a very small amount of error. We also performed a test case simulation of the smaller system size as described in Ref. 9, and we found the ion and water density profiles to be in good agreement with the EW3DC results presented there. With the chosen mesh parameters, our P3M algorithm was approximately four times faster than the corresponding Ewald algorithm. It should be noted that the flexibility of P3M allows for more exact force calculations if desired. The use of a tighter mesh would clearly be more precise, but more CPU intensive. III. SIMULATION DETAILS Simulations were performed at three ion concentrations, 0, 0.25, and 1 M, each at four electrode surface charges, 0, 6 0.1, 6 0.2, and 6 0.3 C/ m2. Previous studies that included discrete solvent molecules considered only higher ionic con-centrations. In all cases, an equal number of monovalent an-ions and cations of equal size were used, and opposite walls were given equal and opposite charges for an electrostati-cally neutral system. The well- known SPC/ E model of water was used, 31 and the mobile ions were given LJ parameters identical to the SPC/ E oxygen oxygen parameters. Initially, the water molecules and ions were randomly distributed FIG. 1. Comparison of the force acting between two oppositely charged point charges in a two- dimensional periodic system as calculated using EW3DC line ! , and P3MC diamonds ! . The scale is chosen for comparison with Fig. 8 of Ref. 7, Fig. 4 of Ref. 8, and Fig. 3 of Ref. 9. FIG. 2. Comparison of the y- component of force parallel to the slab sur-face ! acting between two oppositely charged point charges as calculated using EW3DC line ! and P3MC diamonds ! . The scale is chosen for com-parison with Fig. 4 of Ref. 9. J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Electrolyte systems between charged electrodes 9203 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp throughout the occupied half of the simulation cell between the two oppositely charged electrode surfaces, which con-sisted of single layers of a square lattice of fixed- in- space wall ions centered at 1.5 intervals. Each wall ion was assigned LJ parameters of s LJ 5 1.5 and e / k 5 50 K, along with a partial charge corresponding to the desired electrode charge density. Canonical ensemble NVT ! simulations were performed with the temperature set at 25 C and the volume chosen to give a bulk density of approximately 60 mol/ L. A somewhat high density was chosen in order to ensure that bulk liquid water would be present in all of the cases without the need to adjust the volume for the individual runs. The simulations were allowed to equilibrate for 50 ps prior to the density profile accumulation period of 200 ps 80 000 time steps of 2.5 fs each ! . Twenty repetitions, each with a unique starting configuration, were performed at each state point of interest in order to obtain adequate ensemble averaging and smooth ion density distribution profiles. Each repetition required ap-proximately 16 h of CPU time on a 32 node SGI Origin 2000 supercomputer. IV. RESULTS AND DISCUSSION We present in Figs. 3 5 water and ion z- direction per-pendicular to the electrode ! density distribution profiles for the 12 cases studied at various surface charges, s . No ions are present in Fig. 3, and the water density profiles are omit-ted from Figs. 4 and 5 for clarity. Water density profiles for Figs. 4 and 5 are essentially identical to the water density profiles of Fig. 3 with the corresponding electrode surface charge. In all cases, the positively charged electrode is on the left, and the negatively charged electrode is on the right. Flat, neutral density profiles seen in the central regions of all of the plots indicate the presence of completely shielded bulk electrolyte fluid. The presence of the bulk fluid shows that the positive and negative electrodes are isolated, and enables us to examine the entire shielding layer of fluid at either electrode. 2,28 The oxygen and hydrogen profiles are quite smooth in the central region, but the ionic profiles are less smooth due to the greater statistical uncertainty. Other au-thors have also noted this problem. Even in the absence of ions or electrode surface charge, we see a pronounced difference between bulk and interfacial water characteristics see Fig. 3, top panel ! . Other studies confirm this finding, 1,24,32 34 and point out that the difference is even more pronounced as the electrode is charged. 2,18,32 34 Especially encouraging is the agreement between these mo-lecular dynamics MD ! results and force measurements of Israelachvili et al. 17 and the x- ray scattering findings of Toney et al. 18 They found, contrary to commonly used theo-ries, such as the Gouy Chapman GC ! theory, 35,36 that the water is ordered in layers extending several molecular diam-eters from the electrode. They also observed strong dipole orientation near the charged surfaces. The oxygen and hydro-gen water profiles presented here agree, at least qualitatively, on these points. Earlier theoretical studies of electrochemical interfaces using simpler molecular solvents also show this effect. 32 34 Despite the generally good agreement between these MD results and the experimental findings of Toney et al., 18 we also note a disagreement concerning the density of water next to the charged surface. While the experimental findings indicate that the average density of water across the first peak and valley next to the surface is much higher than the bulk water density, our results indicate no such effect. This is in agreement with other MD simulation results37,38 even though different electrode models were used. We refer the reader to the article by Yeh and Berkowitz38 for a discussion of this discrepancy between the simulation and experimental results. FIG. 3. Oxygen solid line ! and hydrogen broken line ! average density distribution as a function of distance from the center of the simulation cell. The positive and negative electrodes are centered at 2 21.28 and 21.28 , respectively. Bulk anion and cation concentrations are 0 M. 9204 J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Crozier, Rowley, and Henderson Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp Without ions present, one might expect a pure solvent to behave similarly at equally and oppositely charged elec-trodes, but the lack of symmetry of the water molecules about all axes of rotation introduces asymmetry in the den-sity profiles, as can be seen by comparing the left and right hand sides of the plots of Fig. 3. Surprisingly, the addition of ions does little to change the water density distributions. Except in the extreme case of 6 0.3 C/ m2, the ions make very small contributions to the charge shielding of the electrode when compared to the con-tribution of the water molecules. Potential plots at varying ion concentration of the same moderate electrode charge show insignificant differences see Fig. 6 ! . Only for the physically extreme case of 6 0.3 C/ m2 is a significant devia-tion with change in ion concentration observed. It is well-known from experiment39 that the potential and capacitance show strong ion concentration dependence as well as elec-trode charge dependence at low ionic concentrations. This phenomenon was not observed in our simulations. It may be necessary to go to concentrations well below 0.25 M to ob- FIG. 4. Cation solid line ! and anion broken line ! average density distribu-tion as a function of distance from the center of the simulation cell. The positive and negative electrodes are centered at 2 21.28 and 21.28 , re-spectively. Bulk anion and cation concentrations are 0.25 M. FIG. 5. Cation solid line ! and anion broken line ! average density distribu-tion as a function of distance from the center of the simulation cell. The positive and negative electrodes are centered at 2 21.28 and 21.28 , re-spectively. Bulk anion and cation concentrations are 1 M. J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Electrolyte systems between charged electrodes 9205 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp serve this phenomenon. This, regrettably is exceedingly de-manding of computer resources. It is also well known that an ion s size has a strong influence on its solubility, and hence its adsorptivity onto the surface of an electrode. The results shown here are all for equal- sized ions; the effect of ion size on solvation, surface adsorption, and the ability of the ion to affect the double layer capacitance is the subject of a future study. Here, we are interested in the effects of charge and ion concentration. Even though the size of the cations and anions was chosen to be equal for the purposes of this study, the asymmetric charge distribution of the water molecules causes asymmet-ric anion and cation behavior. The cations are clearly more easily adsorbed onto the electrode surface, which also sug-gests better anion solvation, and slower anion transport. At zero charge we observe no wall contact adsorption of either species, and at 6 0.3 C/ m2 both species are strongly ad-sorbed to the electrode surfaces. It is expected that this ob-servation may change when the cations and anions have dif-ferent diameters. The effect of ion size difference will be considered in a subsequent study. Also evident at large surface charges is the well- known fact that the concentration of ions at a charged electrode can be much larger than the bulk ion concentration. 18,19,32 36,40 The bottom panes of Figs. 4 and 5, which correspond to an electrode charge of 6 0.3 C/ m2, show peaks of electrode-adsorbed ions at concentrations of more than two orders- of-magnitude larger than the corresponding bulk ion concentra-tions. Figure 7 shows the potential drop across the double layer as a function of distance from the electrode surface. Theo-ries, such as the GC theory, that do not consider the molecu-lar nature of the solvent, are unable to predict the oscillatory pattern seen here. More sophisticated theories do show this effect. 32 34 Even in the case of no charge and no ions, this oscillatory pattern is clearly evident as seen in Fig. 7. A comparison between the bottom two panes of Fig. 7 is also noteworthy, since it demonstrates how the asymmetric nature of the solvent molecules dramatically affects the potential drop across the double layer. Even though the electrode sur-face charge is equal and opposite, the plots are radically different. Although the potential drop across the double layer seems to be only mildly affected by ion concentration, it is strongly affected by electrode charge see Fig. 6 ! . FIG. 6. Potential drop across the double layer from the middle of the slab to the electrode surface ! as a function of surface charge density and bulk ion concentration. The n , 3 , and h represent 0, 0.25, and 1 M bulk anion and cation concentrations, respectively. FIG. 7. Potential drop from the neutral bulk fluid to the electrode surface as a function of distance from the electrode surface. Bulk anion and cation concentrations are 0 M. 9206 J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Crozier, Rowley, and Henderson Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp V. CONCLUSIONS The P3MC method is a well- defined and efficient way to model coulombic interactions for large slab geometry elec-trolytic systems. Its flexibility and adaptability make it an attractive way to properly account for long- range intermo-lecular interactions at the accuracy level and speed that the user desires. This method has enabled us to accurately per-form these MD calculations within a reasonable amount of computer time. Unlike many previous double layer simulations, the mo-lecular nature of the solvent was not neglected here, and it was found that the water molecules were responsible for most of the potential drop across the double layer and that ion concentration had little effect. Due to the charge distribution and shape of the water molecule, very different behavior was observed at equally and oppositely charged electrodes. With the models used here, the cations more easily adsorbed to the electrode sur-face than the anions, which indicates better solvation of the anions by the water molecules. Given sufficient electrode charge, contact adsorption of either species occurred. In fact, at very high charge, the ion concentration profile shows a peak at the electrodes in excess of 100 times the value of the bulk ion concentration. Overall, results of this MD study were in good qualitative agreement with experimental find-ings, and those found in earlier, less sophisticated theoretical results. ACKNOWLEDGMENTS This work was supported in part by the National Science Foundation Grant No. CHE98- 13729 ! and by the donors of the Petroleum Research Fund, administered by the American Chemical Society Grant No. ACS- PRF 31573- AC9 ! . 1 C. Y. Lee, J. A. McCammon, and P. J. Rossky, J. Chem. Phys. 80, 4448 1984 ! . 2 S. H. Lee and P. J. Rossky, J. Chem. Phys. 100, 3334 1994 ! . 3 J. N. Glosli and M. R. Philpott, J. Chem. Phys. 96, 6962 1992 ! . 4 D. A. Rose and I. Benjamin, J. Chem. Phys. 95, 6856 1991 ! . 5 D. A. Rose and I. Benjamin, J. Chem. Phys. 98, 2283 1993 ! . 6 B. Eck and E. Spohr, Electrochim. Acta 42, 2779 1997 ! . 7 E. Spohr, J. Chem. Phys. 107, 6342 1997 ! . 8 I.- C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 1999 ! . 9 P. S. Crozier, R. L. Rowley, E. Spohr, and D. Henderson, J. Chem. Phys. 112, 9253 2000 ! . 10 R. W. Hockney and J. W. Eastwood, Computer Simulation Using Par-ticles Institute of Physics, Bristol, 1988 ! . 11 T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 1993 ! . 12 B. A. Brock, G. T. Tironi, and W. F. van Gunsteren, J. Chem. Phys. 103, 3014 1995 ! . 13 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 1995 ! . 14M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 1998 ! . 15M. Deserno and C. Holm, J. Chem. Phys. 109, 7694 1998 ! . 16 E. L. Pollock and J. Glosli, Comput. Phys. Commun. 95, 93 1996 ! . 17 J. N. Israelachvili, Intermolecular and Surface Forces, 2nd ed. Aca-demic, London, 1992 ! . 18M. F. Toney, J. N. Howard, J. Richer, G. L. Borges, J. G. Gordon, O. R. Melroy, D. G. Wiesler, D. Yee, and L. B. Sorensen, Surf. Sci. 335, 326 1995 ! . 19 G. M. Torrie and J. P. Valleau, J. Chem. Phys. 73, 5807 1980 ! . 20 D. Boda, K.- Y. Chan, and D. Henderson, J. Chem. Phys. 109, 7362 1998 ! . 21 J. A. Barker and R. O. Watts, Chem. Phys. Lett. 3, 144 1969 ! . 22 R. O. Watts, Mol. Phys. 28, 1069 1974 ! . 23 N. I. Christou, J. S. Whitehouse, D. Nicholson, and N. G. Parsonage, Mol. Phys. 55, 397 1985 ! . 24 S.- B. Zhu and G. W. Robinson, J. Chem. Phys. 94, 1403 1991 ! . 25 D. E. Parry, Surf. Sci. 49, 433 1975 ! ; 54, 195 1976 ! . 26 D. M. Heyes, M. Barber, and J. H. R. Clarke, J. Chem. Soc., Faraday Trans. 2 73, 1485 1977 ! . 27 S. W. de Leeuw and J. W. Perram, Mol. Phys. 37, 1313 1979 ! . 28 L. Greengard and V. Rokhlin, J. Comput. Phys. 60, 187 1985 ! . 29M. R. Philpott and J. N. Glosli, J. Electrochem. Soc. 142, L25 1995 ! . 30 F. Figueirido, R. M. Levy, R. Zhou, and B. J. Berne, J. Chem. Phys. 106, 9835 1997 ! . 31 H. J. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 1987 ! . 32 L. Blum and D. Henderson, J. Chem. Phys. 74, 1902 1981 ! . 33 S. L. Carnie and D. Y. C. Chan, J. Chem. Phys. 73, 2949 1980 ! . 34 D. Henderson, F. F. Abraham, and J. A. Barker, Mol. Phys. 31, 1291 1976 ! . 35 G. Gouy, J. Phys. Paris ! 9, 457 1910 ! . 36 D. L. Chapman, Philos. Mag. 25, 475 1913 ! . 37 K. J. Schweighofer, X. Xia, and M. L. Berkowitz, Langmuir 12, 3747 1996 ! . 38 I.- C. Yeh and M. L. Berkowitz, Chem. Phys. Lett. 301, 81 1999 ! . 39 D. C. Grahame, J. Chem. Phys. 21, 1054 1953 ! . 40 I. Danielewicz- Ferchmin, A. R. Ferchmin, and A. Szlaferek, Chem. Phys. Lett. 288, 197 1998 ! . J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Electrolyte systems between charged electrodes 9207 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp

Click tabs to swap between content that is broken into logical sections.

We investigate the interfacial electrochemical properties of an aqueous electrolyte solution with discrete water molecules in slab geometry between charged atomistic electrodes. Long-range intermolecular Coulombic interactions are calculated using the particle particle particle mesh method with a modification to account for the slab geometry. Density distribution profiles and potential drops across the double layer are given for 0, 0.25, and 1 M aqueous electrolyte solutions each electrode surface charges. Results are compared qualitatively with experimental x-ray scattering findings, other computer simulation results, and traditional electrochemistry theory. The interfacial fluid structure characteristics are generally in good qualitative agreement with the conclusions obtained in some integral equation theories and in the experimental x-ray study. The potential in the simulations shows an oscillatory behavior near the electrode, which theories that do not include the molecular nature of water cannot reproduce for the given conditions. Surprisingly, the results also show that the water structure near the electrode is dominated by the charge on the electrode and is fairly insensitive to the ion concentrations. Except at large electrode charge, the potential drop across the double layer does not depend significantly upon the concentration of the ions.

(c) 2000 American Institute of Physics. This article may be downloaded for personal use only. Any other use requires prior permission of the author and the American Institute of Physics. The following article appeared in The Journal of Chemical Physics and may be found at http://link.aip.org/link/?JCPSA6/118/5474/1;

Molecular dynamics calculations of the electrochemical properties of electrolyte systems between charged electrodes
We investigate the interfacial electrochemical properties of an aqueous electrolyte solution with discrete water molecules in slab geometry between charged atomistic electrodes. Long-range intermolecular Coulombic interactions are calculated using the particle particle particle mesh method with a modification to account for the slab geometry. Density distribution profiles and potential drops across the double layer are given for 0, 0.25, and 1 M aqueous electrolyte solutions each electrode surface charges. Results are compared qualitatively with experimental x-ray scattering findings, other computer simulation results, and traditional electrochemistry theory. The interfacial fluid structure characteristics are generally in good qualitative agreement with the conclusions obtained in some integral equation theories and in the experimental x-ray study. The potential in the simulations shows an oscillatory behavior near the electrode, which theories that do not include the molecular nature of water cannot reproduce for the given conditions. Surprisingly, the results also show that the water structure near the electrode is dominated by the charge on the electrode and is fairly insensitive to the ion concentrations. Except at large electrode charge, the potential drop across the double layer does not depend significantly upon the concentration of the ions.
Molecular dynamics calculations of the electrochemical properties of electrolyte systems between charged electrodes Paul S. Croziera) and Richard L. Rowleyb) Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602- 4100 Douglas Hendersonc) Department of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602- 5700 Received 27 June 2000; accepted 1 September 2000 ! We investigate the interfacial electrochemical properties of an aqueous electrolyte solution with discrete water molecules in slab geometry between charged atomistic electrodes. Long- range intermolecular Coulombic interactions are calculated using the particle particle particle mesh method with a modification to account for the slab geometry. Density distribution profiles and potential drops across the double layer are given for 0, 0.25, and 1 M aqueous electrolyte solutions each at 0, 6 0.1, 6 0.2, and 6 0.3 C/ m2 electrode surface charges. Results are compared qualitatively with experimental x- ray scattering findings, other computer simulation results, and traditional electrochemistry theory. The interfacial fluid structure characteristics are generally in good qualitative agreement with the conclusions obtained in some integral equation theories and in the experimental x- ray study. The potential in the simulations shows an oscillatory behavior near the electrode, which theories that do not include the molecular nature of water cannot reproduce for the given conditions. Surprisingly, the results also show that the water structure near the electrode is dominated by the charge on the electrode and is fairly insensitive to the ion concentrations. Except at large electrode charge, the potential drop across the double layer does not depend significantly upon the concentration of the ions. 2000 American Institute of Physics. @ S0021- 9606 00 ! 51344- 0 # I. INTRODUCTION The ability to model correctly electrochemical processes at interfaces is a key to our understanding of numerous very important chemical and biological processes. Especially of interest are aqueous solutions in contact with charged sur-faces, which can be found in such diverse fields as molecular biology, electrochemistry, and heterogeneous catalysis. Although unable to model quantum mechanical effects, such as electron transfer, classical molecular dynamics is the most valuable tool available for the elucidation of important interfacial processes. Such simulations provide realistic rep-resentations of interfacial electrochemistry processes, limited primarily by molecular model accuracy and computational demands. In order to save time, many interfacial electrochemistry molecular dynamics researchers neglect long- range Coulom-bic interactions while admitting the vital role that these play. 1 5 Other researchers have employed costly Ewald sum-mation approximations for long- range interactions. 6 9 Re-cently, mesh routines have been developed, 10 15 which present a more efficient way to perform these calculations. Of the mesh routines, the particle particle particle mesh P3M ! routine developed by Hockney and Eastwood10 has been shown to be the most efficient for a given accuracy level. 14,15 In this work, we modify P3M for application to slab geometry, as suggested by Pollock and Glosli16 and similar to the 3D Ewald summation modification used by Yeh and Berkowitz. 8 We show that this implementation of the slab- corrected P3M method, or P3MC as we shall call it, gives results in very good agreement with regular Ewald summation methods, yet at a fraction of the cost. Once we have established the validity of the P3MC methodology, we present simulation results obtained using this method dealing with the effect of electrode charge and ion concentration on the potential drop across the double layer of an aqueous electrolyte solution at a solid atomistic surface. We also note the similarity of our results with the experimental force measurements of Israelachvili et al. 17 and the x- ray scattering findings of Toney et al. 18 In addition, we compare our findings with the predictions of traditional in-terfacial electrochemistry theories, as well as with the find-ings from other simulation work. II. P3MC IMPLEMENTATION Many methods for approximating long- range Coulombic interactions have been devised, including the charged sheets method, 19,20 the modified reaction field, 21 24 Ewald sums, 6 9,25 27 the fast multipole method, 16,28 30 and mesh methods. 10 15 Recently, Crozier et al. 9 found that the charged sheets method is inadequate for correctly modeling strong long- range Coulombic interactions for systems with discrete solvent molecules. Ewald summations appear to be adequate, yet computationally expensive, and multipole methods are even more costly for systems of a reasonable size. 16 Mesh methods, particularly P3M, present the most flexible and reliable alternative. 14 16 Deserno and Holm give a ! Electronic mail: crozierp@ et. byu. edu b ! Electronic mail: rowley@ byu. edu c ! Electronic mail: doug@ huey. byu. edu JOURNAL OF CHEMICAL PHYSICS VOLUME 113, NUMBER 20 22 NOVEMBER 2000 0021- 9606/ 2000/ 113( 20)/ 9202/ 6/$ 17.00 9202 2000 American Institute of Physics Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp an excellent analysis of the various mesh methods, 14 along with useful recommendations for P3M implementation. 15 We have followed their recommendations and adopted their no-tation in the following. We use Fourier space, or ik, differentiation to calculate the intermolecular forces from the electrostatic energies, and a seventh- order assignment scheme ( P 5 7) to assign the charges to the mesh. Calculation of the optimal influence function, G opt , of Hockney and Eastwood10 is done at the beginning of each run and therefore presents no additional run- time overhead. The cutoff for the Lennard- Jones LJ ! interactions, and the real- space Coulombic interactions is set at rmax 5 10 , and the mesh size is set at 16 3 16 3 64 in the x, y, and z directions, respectively. With a corresponding fixed box size of 25.5 3 25.5 3 85.1232 , and the total number of ions and water molecules fixed at 1000 the opti-mum value of a is estimated to be 0.3007. We estimate the average error in our force calculations to be on the order of 10 2 5 dimensionless force units, according to the estimation method and dimensionless force units of Deserno and Holm. 14 Although not as precise as an Ewald sum, the amount of error introduced by this implementation of P3M is very reasonable. Extreme accuracy is not needed, especially when considering other possible sources of error in the MD simulation, such as the use of a randomly fluctuating ther-mostat, or a discrete time step e. g., 2.5 fs ! . In order to account for the slab geometry, we include an empty space, equal to the size of the fluid- occupied space, between repeating slabs of fluid, and a correction term to damp out interslab interactions. With the implementation of this correction, the 3D, periodically repeating simulation closely approximates the simulation of an isolated 2D slab of finite thickness, as shown by Yeh and Berkowitz. 8 To confirm the reliability of our P3MC implementation, we compare calculated forces with those calculated using the corrected 3D Ewald EW3DC ! method. 8 In particular, we calculate the force between two oppositely charged particles along with their periodic images in slab geometry. For easy comparison, we consider the same case as presented in Refs. 7 9 ( L 5 Lx 5 Ly 5 18 and Lz 5 90 ! . As seen in Figs. 1 and 2, the P3MC and EW3DC calculations are visually in-distinguishable. For Fig. 1, the difference between the two calculations was less than 0.013% in all cases, with an aver-age difference of 0.002%. Even though the deviation was up to 1% in the case of Fig. 2, percent deviation is less mean-ingful here, since the absolute value of the force approaches zero. The absolute difference between the two calculations on the scale of Fig. 2 was always less than 10 2 4, which represents a very small amount of error. We also performed a test case simulation of the smaller system size as described in Ref. 9, and we found the ion and water density profiles to be in good agreement with the EW3DC results presented there. With the chosen mesh parameters, our P3M algorithm was approximately four times faster than the corresponding Ewald algorithm. It should be noted that the flexibility of P3M allows for more exact force calculations if desired. The use of a tighter mesh would clearly be more precise, but more CPU intensive. III. SIMULATION DETAILS Simulations were performed at three ion concentrations, 0, 0.25, and 1 M, each at four electrode surface charges, 0, 6 0.1, 6 0.2, and 6 0.3 C/ m2. Previous studies that included discrete solvent molecules considered only higher ionic con-centrations. In all cases, an equal number of monovalent an-ions and cations of equal size were used, and opposite walls were given equal and opposite charges for an electrostati-cally neutral system. The well- known SPC/ E model of water was used, 31 and the mobile ions were given LJ parameters identical to the SPC/ E oxygen oxygen parameters. Initially, the water molecules and ions were randomly distributed FIG. 1. Comparison of the force acting between two oppositely charged point charges in a two- dimensional periodic system as calculated using EW3DC line ! , and P3MC diamonds ! . The scale is chosen for comparison with Fig. 8 of Ref. 7, Fig. 4 of Ref. 8, and Fig. 3 of Ref. 9. FIG. 2. Comparison of the y- component of force parallel to the slab sur-face ! acting between two oppositely charged point charges as calculated using EW3DC line ! and P3MC diamonds ! . The scale is chosen for com-parison with Fig. 4 of Ref. 9. J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Electrolyte systems between charged electrodes 9203 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp throughout the occupied half of the simulation cell between the two oppositely charged electrode surfaces, which con-sisted of single layers of a square lattice of fixed- in- space wall ions centered at 1.5 intervals. Each wall ion was assigned LJ parameters of s LJ 5 1.5 and e / k 5 50 K, along with a partial charge corresponding to the desired electrode charge density. Canonical ensemble NVT ! simulations were performed with the temperature set at 25 C and the volume chosen to give a bulk density of approximately 60 mol/ L. A somewhat high density was chosen in order to ensure that bulk liquid water would be present in all of the cases without the need to adjust the volume for the individual runs. The simulations were allowed to equilibrate for 50 ps prior to the density profile accumulation period of 200 ps 80 000 time steps of 2.5 fs each ! . Twenty repetitions, each with a unique starting configuration, were performed at each state point of interest in order to obtain adequate ensemble averaging and smooth ion density distribution profiles. Each repetition required ap-proximately 16 h of CPU time on a 32 node SGI Origin 2000 supercomputer. IV. RESULTS AND DISCUSSION We present in Figs. 3 5 water and ion z- direction per-pendicular to the electrode ! density distribution profiles for the 12 cases studied at various surface charges, s . No ions are present in Fig. 3, and the water density profiles are omit-ted from Figs. 4 and 5 for clarity. Water density profiles for Figs. 4 and 5 are essentially identical to the water density profiles of Fig. 3 with the corresponding electrode surface charge. In all cases, the positively charged electrode is on the left, and the negatively charged electrode is on the right. Flat, neutral density profiles seen in the central regions of all of the plots indicate the presence of completely shielded bulk electrolyte fluid. The presence of the bulk fluid shows that the positive and negative electrodes are isolated, and enables us to examine the entire shielding layer of fluid at either electrode. 2,28 The oxygen and hydrogen profiles are quite smooth in the central region, but the ionic profiles are less smooth due to the greater statistical uncertainty. Other au-thors have also noted this problem. Even in the absence of ions or electrode surface charge, we see a pronounced difference between bulk and interfacial water characteristics see Fig. 3, top panel ! . Other studies confirm this finding, 1,24,32 34 and point out that the difference is even more pronounced as the electrode is charged. 2,18,32 34 Especially encouraging is the agreement between these mo-lecular dynamics MD ! results and force measurements of Israelachvili et al. 17 and the x- ray scattering findings of Toney et al. 18 They found, contrary to commonly used theo-ries, such as the Gouy Chapman GC ! theory, 35,36 that the water is ordered in layers extending several molecular diam-eters from the electrode. They also observed strong dipole orientation near the charged surfaces. The oxygen and hydro-gen water profiles presented here agree, at least qualitatively, on these points. Earlier theoretical studies of electrochemical interfaces using simpler molecular solvents also show this effect. 32 34 Despite the generally good agreement between these MD results and the experimental findings of Toney et al., 18 we also note a disagreement concerning the density of water next to the charged surface. While the experimental findings indicate that the average density of water across the first peak and valley next to the surface is much higher than the bulk water density, our results indicate no such effect. This is in agreement with other MD simulation results37,38 even though different electrode models were used. We refer the reader to the article by Yeh and Berkowitz38 for a discussion of this discrepancy between the simulation and experimental results. FIG. 3. Oxygen solid line ! and hydrogen broken line ! average density distribution as a function of distance from the center of the simulation cell. The positive and negative electrodes are centered at 2 21.28 and 21.28 , respectively. Bulk anion and cation concentrations are 0 M. 9204 J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Crozier, Rowley, and Henderson Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp Without ions present, one might expect a pure solvent to behave similarly at equally and oppositely charged elec-trodes, but the lack of symmetry of the water molecules about all axes of rotation introduces asymmetry in the den-sity profiles, as can be seen by comparing the left and right hand sides of the plots of Fig. 3. Surprisingly, the addition of ions does little to change the water density distributions. Except in the extreme case of 6 0.3 C/ m2, the ions make very small contributions to the charge shielding of the electrode when compared to the con-tribution of the water molecules. Potential plots at varying ion concentration of the same moderate electrode charge show insignificant differences see Fig. 6 ! . Only for the physically extreme case of 6 0.3 C/ m2 is a significant devia-tion with change in ion concentration observed. It is well-known from experiment39 that the potential and capacitance show strong ion concentration dependence as well as elec-trode charge dependence at low ionic concentrations. This phenomenon was not observed in our simulations. It may be necessary to go to concentrations well below 0.25 M to ob- FIG. 4. Cation solid line ! and anion broken line ! average density distribu-tion as a function of distance from the center of the simulation cell. The positive and negative electrodes are centered at 2 21.28 and 21.28 , re-spectively. Bulk anion and cation concentrations are 0.25 M. FIG. 5. Cation solid line ! and anion broken line ! average density distribu-tion as a function of distance from the center of the simulation cell. The positive and negative electrodes are centered at 2 21.28 and 21.28 , re-spectively. Bulk anion and cation concentrations are 1 M. J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Electrolyte systems between charged electrodes 9205 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp serve this phenomenon. This, regrettably is exceedingly de-manding of computer resources. It is also well known that an ion s size has a strong influence on its solubility, and hence its adsorptivity onto the surface of an electrode. The results shown here are all for equal- sized ions; the effect of ion size on solvation, surface adsorption, and the ability of the ion to affect the double layer capacitance is the subject of a future study. Here, we are interested in the effects of charge and ion concentration. Even though the size of the cations and anions was chosen to be equal for the purposes of this study, the asymmetric charge distribution of the water molecules causes asymmet-ric anion and cation behavior. The cations are clearly more easily adsorbed onto the electrode surface, which also sug-gests better anion solvation, and slower anion transport. At zero charge we observe no wall contact adsorption of either species, and at 6 0.3 C/ m2 both species are strongly ad-sorbed to the electrode surfaces. It is expected that this ob-servation may change when the cations and anions have dif-ferent diameters. The effect of ion size difference will be considered in a subsequent study. Also evident at large surface charges is the well- known fact that the concentration of ions at a charged electrode can be much larger than the bulk ion concentration. 18,19,32 36,40 The bottom panes of Figs. 4 and 5, which correspond to an electrode charge of 6 0.3 C/ m2, show peaks of electrode-adsorbed ions at concentrations of more than two orders- of-magnitude larger than the corresponding bulk ion concentra-tions. Figure 7 shows the potential drop across the double layer as a function of distance from the electrode surface. Theo-ries, such as the GC theory, that do not consider the molecu-lar nature of the solvent, are unable to predict the oscillatory pattern seen here. More sophisticated theories do show this effect. 32 34 Even in the case of no charge and no ions, this oscillatory pattern is clearly evident as seen in Fig. 7. A comparison between the bottom two panes of Fig. 7 is also noteworthy, since it demonstrates how the asymmetric nature of the solvent molecules dramatically affects the potential drop across the double layer. Even though the electrode sur-face charge is equal and opposite, the plots are radically different. Although the potential drop across the double layer seems to be only mildly affected by ion concentration, it is strongly affected by electrode charge see Fig. 6 ! . FIG. 6. Potential drop across the double layer from the middle of the slab to the electrode surface ! as a function of surface charge density and bulk ion concentration. The n , 3 , and h represent 0, 0.25, and 1 M bulk anion and cation concentrations, respectively. FIG. 7. Potential drop from the neutral bulk fluid to the electrode surface as a function of distance from the electrode surface. Bulk anion and cation concentrations are 0 M. 9206 J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Crozier, Rowley, and Henderson Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp V. CONCLUSIONS The P3MC method is a well- defined and efficient way to model coulombic interactions for large slab geometry elec-trolytic systems. Its flexibility and adaptability make it an attractive way to properly account for long- range intermo-lecular interactions at the accuracy level and speed that the user desires. This method has enabled us to accurately per-form these MD calculations within a reasonable amount of computer time. Unlike many previous double layer simulations, the mo-lecular nature of the solvent was not neglected here, and it was found that the water molecules were responsible for most of the potential drop across the double layer and that ion concentration had little effect. Due to the charge distribution and shape of the water molecule, very different behavior was observed at equally and oppositely charged electrodes. With the models used here, the cations more easily adsorbed to the electrode sur-face than the anions, which indicates better solvation of the anions by the water molecules. Given sufficient electrode charge, contact adsorption of either species occurred. In fact, at very high charge, the ion concentration profile shows a peak at the electrodes in excess of 100 times the value of the bulk ion concentration. Overall, results of this MD study were in good qualitative agreement with experimental find-ings, and those found in earlier, less sophisticated theoretical results. ACKNOWLEDGMENTS This work was supported in part by the National Science Foundation Grant No. CHE98- 13729 ! and by the donors of the Petroleum Research Fund, administered by the American Chemical Society Grant No. ACS- PRF 31573- AC9 ! . 1 C. Y. Lee, J. A. McCammon, and P. J. Rossky, J. Chem. Phys. 80, 4448 1984 ! . 2 S. H. Lee and P. J. Rossky, J. Chem. Phys. 100, 3334 1994 ! . 3 J. N. Glosli and M. R. Philpott, J. Chem. Phys. 96, 6962 1992 ! . 4 D. A. Rose and I. Benjamin, J. Chem. Phys. 95, 6856 1991 ! . 5 D. A. Rose and I. Benjamin, J. Chem. Phys. 98, 2283 1993 ! . 6 B. Eck and E. Spohr, Electrochim. Acta 42, 2779 1997 ! . 7 E. Spohr, J. Chem. Phys. 107, 6342 1997 ! . 8 I.- C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 1999 ! . 9 P. S. Crozier, R. L. Rowley, E. Spohr, and D. Henderson, J. Chem. Phys. 112, 9253 2000 ! . 10 R. W. Hockney and J. W. Eastwood, Computer Simulation Using Par-ticles Institute of Physics, Bristol, 1988 ! . 11 T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 1993 ! . 12 B. A. Brock, G. T. Tironi, and W. F. van Gunsteren, J. Chem. Phys. 103, 3014 1995 ! . 13 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 1995 ! . 14M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 1998 ! . 15M. Deserno and C. Holm, J. Chem. Phys. 109, 7694 1998 ! . 16 E. L. Pollock and J. Glosli, Comput. Phys. Commun. 95, 93 1996 ! . 17 J. N. Israelachvili, Intermolecular and Surface Forces, 2nd ed. Aca-demic, London, 1992 ! . 18M. F. Toney, J. N. Howard, J. Richer, G. L. Borges, J. G. Gordon, O. R. Melroy, D. G. Wiesler, D. Yee, and L. B. Sorensen, Surf. Sci. 335, 326 1995 ! . 19 G. M. Torrie and J. P. Valleau, J. Chem. Phys. 73, 5807 1980 ! . 20 D. Boda, K.- Y. Chan, and D. Henderson, J. Chem. Phys. 109, 7362 1998 ! . 21 J. A. Barker and R. O. Watts, Chem. Phys. Lett. 3, 144 1969 ! . 22 R. O. Watts, Mol. Phys. 28, 1069 1974 ! . 23 N. I. Christou, J. S. Whitehouse, D. Nicholson, and N. G. Parsonage, Mol. Phys. 55, 397 1985 ! . 24 S.- B. Zhu and G. W. Robinson, J. Chem. Phys. 94, 1403 1991 ! . 25 D. E. Parry, Surf. Sci. 49, 433 1975 ! ; 54, 195 1976 ! . 26 D. M. Heyes, M. Barber, and J. H. R. Clarke, J. Chem. Soc., Faraday Trans. 2 73, 1485 1977 ! . 27 S. W. de Leeuw and J. W. Perram, Mol. Phys. 37, 1313 1979 ! . 28 L. Greengard and V. Rokhlin, J. Comput. Phys. 60, 187 1985 ! . 29M. R. Philpott and J. N. Glosli, J. Electrochem. Soc. 142, L25 1995 ! . 30 F. Figueirido, R. M. Levy, R. Zhou, and B. J. Berne, J. Chem. Phys. 106, 9835 1997 ! . 31 H. J. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 1987 ! . 32 L. Blum and D. Henderson, J. Chem. Phys. 74, 1902 1981 ! . 33 S. L. Carnie and D. Y. C. Chan, J. Chem. Phys. 73, 2949 1980 ! . 34 D. Henderson, F. F. Abraham, and J. A. Barker, Mol. Phys. 31, 1291 1976 ! . 35 G. Gouy, J. Phys. Paris ! 9, 457 1910 ! . 36 D. L. Chapman, Philos. Mag. 25, 475 1913 ! . 37 K. J. Schweighofer, X. Xia, and M. L. Berkowitz, Langmuir 12, 3747 1996 ! . 38 I.- C. Yeh and M. L. Berkowitz, Chem. Phys. Lett. 301, 81 1999 ! . 39 D. C. Grahame, J. Chem. Phys. 21, 1054 1953 ! . 40 I. Danielewicz- Ferchmin, A. R. Ferchmin, and A. Szlaferek, Chem. Phys. Lett. 288, 197 1998 ! . J. Chem. Phys., Vol. 113, No. 20, 22 November 2000 Electrolyte systems between charged electrodes 9207 Downloaded 16 Sep 2009 to 128.187.0.164. Redistribution subject to AIP license or copyright; see http:// jcp. aip. org/ jcp/ copyright. jsp