Simulating an electrochemical interface using charge dynamics

We present a simple classical method for treating charge mobility in metals adjacent to liquid solutions. The method, known as electrode charge dynamics, effectively bridges the computational gap between ab initio calculations on small metal clusters and large-scale simulations of metal surfaces...

Ausführliche Beschreibung

Gespeichert in:
Bibliographische Detailangaben
Veröffentlicht in:Condensed Matter Physics
Datum:2005
Hauptverfasser: Guymon, C.G., Rowley, R.L., Harb, J.N., Wheeler, D.R.
Format: Artikel
Sprache:English
Veröffentlicht: Інститут фізики конденсованих систем НАН України 2005
Online Zugang:https://nasplib.isofts.kiev.ua/handle/123456789/119605
Tags: Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
Назва журналу:Digital Library of Periodicals of National Academy of Sciences of Ukraine
Zitieren:Simulating an electrochemical interface using charge dynamics / C.G. Guymon, R.L. Rowley, J.N. Harb, D.R. Wheeler // Condensed Matter Physics. — 2005. — Т. 8, № 2(42). — С. 335–356. — Бібліогр.: 29 назв. — англ.

Institution

Digital Library of Periodicals of National Academy of Sciences of Ukraine
id nasplib_isofts_kiev_ua-123456789-119605
record_format dspace
spelling Guymon, C.G.
Rowley, R.L.
Harb, J.N.
Wheeler, D.R.
2017-06-07T15:38:51Z
2017-06-07T15:38:51Z
2005
Simulating an electrochemical interface using charge dynamics / C.G. Guymon, R.L. Rowley, J.N. Harb, D.R. Wheeler // Condensed Matter Physics. — 2005. — Т. 8, № 2(42). — С. 335–356. — Бібліогр.: 29 назв. — англ.
1607-324X
PACS: 61.20.Qg, 61.20.Ja
DOI:10.5488/CMP.8.2.335
https://nasplib.isofts.kiev.ua/handle/123456789/119605
We present a simple classical method for treating charge mobility in metals adjacent to liquid solutions. The method, known as electrode charge dynamics, effectively bridges the computational gap between ab initio calculations on small metal clusters and large-scale simulations of metal surfaces with arbitrary geometry. We have obtained model parameters for a copper (111) metal surface using high-level quantum-mechanical calculations on a 10-atom copper cluster. We validated the model against the classical image-charge result and ab initio results on an 18-atom copper cluster. The model is used in molecular dynamics simulations to predict the structure of the fluid interface for neat water and for aqueous NaCl solution. We find that water is organized into a two-dimensional ice-like layer on the surface and that both Na⁺ and Cl⁻ are strongly bound to the copper. When charging the metal electrode, most of the electrolyte response occurs in the diffuse part of the double layer
Ми представляємо простий класичний метод для трактування зарядової мобільності у металах, які межують з рідкими розчинами. Метод, відомий як електродна зарядова динаміка, ефективно заповнює прогалину між ab initio розрахунками на малих металічних кластерах і велико-масштабними симуляціями металічних поверхонь з довільною геометрією. Ми отримали модельні параметри для металічної (111) поверхні міді, використовуючи квантово-механічні розрахунки високого порядку на 10-атомному мідному кластері. Модель була перевірена шляхом порівняння з класичними результатами по методу відображень, а також з ab initio результатами для 18-атомного мідного кластера. Модель використана в молекулярно динамічних розрахунках для структури флюїдної поверхні чистої води і водного розчину NaCl. Ми побачили, що вода організовується у двовимірний льодоподібний шар на поверхні і що обидва, Na⁺ i Cl⁻, іони є сильно прив’язані до міді. Коли металічний електрод заряджати, то основна реакція електроліту проявляється у дифузивній частині подвійного шару.
We want to thank Dr. Tapani Pakkanen and his students at the University of Joensuu for their work in performing much of the ab initio calculations presented here. We’d also like to acknowledge financial support from the NSF, grant number CTS–0215786 as well as support from Brigham Young University.
en
Інститут фізики конденсованих систем НАН України
Condensed Matter Physics
Simulating an electrochemical interface using charge dynamics
Симуляція електрохімічної поверхні з використанням зарядової динаміки
Article
published earlier
institution Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection DSpace DC
title Simulating an electrochemical interface using charge dynamics
spellingShingle Simulating an electrochemical interface using charge dynamics
Guymon, C.G.
Rowley, R.L.
Harb, J.N.
Wheeler, D.R.
title_short Simulating an electrochemical interface using charge dynamics
title_full Simulating an electrochemical interface using charge dynamics
title_fullStr Simulating an electrochemical interface using charge dynamics
title_full_unstemmed Simulating an electrochemical interface using charge dynamics
title_sort simulating an electrochemical interface using charge dynamics
author Guymon, C.G.
Rowley, R.L.
Harb, J.N.
Wheeler, D.R.
author_facet Guymon, C.G.
Rowley, R.L.
Harb, J.N.
Wheeler, D.R.
publishDate 2005
language English
container_title Condensed Matter Physics
publisher Інститут фізики конденсованих систем НАН України
format Article
title_alt Симуляція електрохімічної поверхні з використанням зарядової динаміки
description We present a simple classical method for treating charge mobility in metals adjacent to liquid solutions. The method, known as electrode charge dynamics, effectively bridges the computational gap between ab initio calculations on small metal clusters and large-scale simulations of metal surfaces with arbitrary geometry. We have obtained model parameters for a copper (111) metal surface using high-level quantum-mechanical calculations on a 10-atom copper cluster. We validated the model against the classical image-charge result and ab initio results on an 18-atom copper cluster. The model is used in molecular dynamics simulations to predict the structure of the fluid interface for neat water and for aqueous NaCl solution. We find that water is organized into a two-dimensional ice-like layer on the surface and that both Na⁺ and Cl⁻ are strongly bound to the copper. When charging the metal electrode, most of the electrolyte response occurs in the diffuse part of the double layer Ми представляємо простий класичний метод для трактування зарядової мобільності у металах, які межують з рідкими розчинами. Метод, відомий як електродна зарядова динаміка, ефективно заповнює прогалину між ab initio розрахунками на малих металічних кластерах і велико-масштабними симуляціями металічних поверхонь з довільною геометрією. Ми отримали модельні параметри для металічної (111) поверхні міді, використовуючи квантово-механічні розрахунки високого порядку на 10-атомному мідному кластері. Модель була перевірена шляхом порівняння з класичними результатами по методу відображень, а також з ab initio результатами для 18-атомного мідного кластера. Модель використана в молекулярно динамічних розрахунках для структури флюїдної поверхні чистої води і водного розчину NaCl. Ми побачили, що вода організовується у двовимірний льодоподібний шар на поверхні і що обидва, Na⁺ i Cl⁻, іони є сильно прив’язані до міді. Коли металічний електрод заряджати, то основна реакція електроліту проявляється у дифузивній частині подвійного шару.
issn 1607-324X
url https://nasplib.isofts.kiev.ua/handle/123456789/119605
citation_txt Simulating an electrochemical interface using charge dynamics / C.G. Guymon, R.L. Rowley, J.N. Harb, D.R. Wheeler // Condensed Matter Physics. — 2005. — Т. 8, № 2(42). — С. 335–356. — Бібліогр.: 29 назв. — англ.
work_keys_str_mv AT guymoncg simulatinganelectrochemicalinterfaceusingchargedynamics
AT rowleyrl simulatinganelectrochemicalinterfaceusingchargedynamics
AT harbjn simulatinganelectrochemicalinterfaceusingchargedynamics
AT wheelerdr simulatinganelectrochemicalinterfaceusingchargedynamics
AT guymoncg simulâcíâelektrohímíčnoípoverhnízvikoristannâmzarâdovoídinamíki
AT rowleyrl simulâcíâelektrohímíčnoípoverhnízvikoristannâmzarâdovoídinamíki
AT harbjn simulâcíâelektrohímíčnoípoverhnízvikoristannâmzarâdovoídinamíki
AT wheelerdr simulâcíâelektrohímíčnoípoverhnízvikoristannâmzarâdovoídinamíki
first_indexed 2025-11-26T02:05:56Z
last_indexed 2025-11-26T02:05:56Z
_version_ 1850607847675527168
fulltext Condensed Matter Physics, 2005, Vol. 8, No. 2(42), pp. 335–356 Simulating an electrochemical interface using charge dynamics C.G.Guymon, R.L.Rowley, J.N.Harb, D.R.Wheeler Department of Chemical Engineering, BYU, Provo UT 84602 Received November 29, 2004 We present a simple classical method for treating charge mobility in met- als adjacent to liquid solutions. The method, known as electrode charge dynamics, effectively bridges the computational gap between ab initio cal- culations on small metal clusters and large-scale simulations of metal sur- faces with arbitrary geometry. We have obtained model parameters for a copper (111) metal surface using high-level quantum-mechanical calcu- lations on a 10-atom copper cluster. We validated the model against the classical image-charge result and ab initio results on an 18-atom copper cluster. The model is used in molecular dynamics simulations to predict the structure of the fluid interface for neat water and for aqueous NaCl so- lution. We find that water is organized into a two-dimensional ice-like layer on the surface and that both Na+ and Cl− are strongly bound to the cop- per. When charging the metal electrode, most of the electrolyte response occurs in the diffuse part of the double layer. Key words: simulation, double layer, molecular dynamics, ab initio, potentials, copper (111) surface, water PACS: 61.20.Qg, 61.20.Ja 1. Introduction Increased understanding of the interface between metals and liquid solutions (i.e. the electrochemical interface) is vital for the design and improvement of catalytic, corrosive, and electrochemical systems. The configuration of liquid species at the in- terface is intimately connected with mechanistic pathways of ion adsorption, metal dissolution, and metal deposition. Atomically detailed simulation techniques such as molecular dynamics have yielded insight into these kinds of problems. Simulati- ons have been used to predict the adsorbed water orientation, ionic distributions, ionic transport coefficients, atomically resolved electrostatic potential profiles, and interfacial differential capacitance [1–6]. Accurate molecular simulations of an electrochemical interface depend on the quality of intermolecular potentials between species. Prior work has resulted in ad- c© C.G.Guymon, R.L.Rowley, J.N.Harb, D.R.Wheeler 335 C.G.Guymon et al. equate interaction potentials between species in bulk liquid electrolytes, however the potentials suitable at a metallic interface are less well established. For an inter- molecular potential to be suitable for use in large-scale MD simulations, it should be computationally inexpensive. In this work we focus on obtaining and using accurate but inexpensive potentials between solution species and the metal. Two common ways of treating a charged species interacting with a polarizable metal are (1) the image-charge method and (2) the use of pairwise additive interatomic potentials obtained from ab initio calculations. Both of these approximations have been used to simulate large atomically discrete systems [7–10]. The image-charge method is based on an analytic solution to Poisson’s equation. When a charged particle approaches a conducting surface, it polarizes the surface, resulting in an attractive force. The induced potential field outside the metal is perfectly mimicked by placing an image point charge beneath the surface of the metal. For a planar interface the image charge is equal and opposite to the “real” charge and is equidistant to the interface. In a molecular simulation that uses im- age charges, all real and image charges interact with each other. Due to its classical continuum derivation, at long range the method correctly reproduces the Coulombic potential between a point charge and a semi-infinite metallic surface. However, when the charge is within about 3 Å of the surface the method overpredicts the strength of interaction and does not resolve important atomic-level details and spatial non- uniformities. For example, the image-charge method predicts the same interaction energy for positive or negative charges, and does not differentiate between surface ad- sorption sites. Moreover, the method requires an ambiguous selection of the location of the image mirror plane. Spohr and coworkers partially remedy these difficulties by incorporating an additional species-surface interatomic potential to account for the surface corrugation [7]. An alternative is to assume that the potential energy between a molecule and the metal surface is well represented by the interaction between the molecule and a small metal cluster (typically 10–20 metal atoms). The advantage of this method is that the molecule-cluster interaction energy can be obtained by quantum chemical calculations. In a refinement to this approach, the molecule-cluster interaction is decomposed into pairwise atomic interactions that can then be used for the molecule interacting with a much larger metal surface during simulation [3,8]. The primary drawback of this method is that it neglects multibody charge induction in the metal, which will be significant in liquid simulations due to the high density of charges next to the metal. This is in addition to concerns about the assumption that a small metal cluster accurately represents a surface [11]. In addition to the two above highlighted treatments of charge dynamics near a metal surface, Allen et al. [12] and Boda et al. [13] have developed and implemented induced charge computation methods. These recent methods allow for image-charge- type results with geometries other than a flat surface, and thus parallel some of the advantages we claim for this work. We present here a simple method, termed electrode charge dynamics (ECD), that reproduces the image-charge result at long range and matches quantum chemical 336 Simulating an electrochemical interface calculations at short range. ECD is superior to using only pairwise additive potentials or image charges or a combination thereof. In particular, ECD is easy to use for irregular surface geometries. In the sections that follow we first describe the method and then validate it by comparing ECD with other calculations. We conclude with MD simulation results for an aqueous NaCl electrolyte and neat water next to a planar (111) copper surface. 2. Methodology The electrode charge dynamics method is based on a simple model of charge density in the metal, shown in figure 1. Metallic conduction-band electrons are free to move in response to Coulombic forces originating within and without the metal, resulting in regions of increased and decreased electron density. We represent the conduction or valence electrons throughout the metal with a diffuse negative charge located at each metal atomic center of variable magnitude qv i . A fixed positive point charge qc i is co-located at each center to represent the nucleus and core electrons. Figure 1. Charge on each metal atom. Each diffuse charge is modeled with a Gaussian charge-density distribution, ρi(r) = qv i γ 3 i π −3/2 exp ( −γ2 i r 2 ) , (1) where r is the distance from the atom center and γi is an inverse-width parameter. We treat charge mobility in the metal by allowing each diffuse charge magnitude to fluctuate in response to its environment. More specifically, each qv i changes so as to minimize the electrostatic energy or equivalently to equilibrate the chemical potential for qv i at each atom. The fluctuation in qv i is subject to two constraints. The first constraint is that qv i 6 0, corresponding to a limit on the depletion of electron density from a given atom. The second constraint is a fixed total charge Qtot in the simulated metal: n ∑ i=1 (qv i + qc i ) = Qtot , (2) where the sum is over all metal atoms. 337 C.G.Guymon et al. At this point we note the similarity between our model and the charge equili- bration schemes developed by Rappé and Goddard [14] and Rick, Stuart, and Berne [15]. In these models as well as in our model there are two adjustable parameters per unique charge site. In this case, since the metal surface is composed of only one kind of atom, there are only two parameters (here qc and γ) that characterize the polari- zation behavior of the metal. The explicit division of the charge into core and valence charges is unique to this work and leads to different behavior for short-range inter- actions. We also note our model’s classical resemblance to the quantum-mechanical approximations used by Price and Halley [16]. To obtain the magnitude of the diffuse charge at each electrode atom, we mini- mize the Coulombic potential energy U of the metal subject to the constraints. This can be done analytically using Lagrange’s method of undetermined multipliers. The Lagrangian to be minimized is L = U − λ [ n ∑ i=1 (qv i + qc i ) − Qtot ] − n ∑ i=1 µiq v i , (3) where λ and µi are undetermined multipliers, and U = 1 2 n ∑ i=1 n ∑ j=1 qv i q v j Cij + n ∑ i=1 n ∑ j=1 qv i q c jC ∗ ij + 1 2 n ∑ i=1 n ∑ j=1, j 6=i qc i q c j 1 rij + n ∑ i=1 qv i φ set i + Uext . (4) Cij and C∗ ij are Coulomb overlap integrals given by Cij = erf(γijrij) rij , C∗ ij = erf(γirij) rij , (5) with γij = (γ−2 i + γ−2 j )−1/2. Uext is the Coulomb energy from interactions of the electrode charges with charges external to the electrode. φset i is a user specified offset potential or voltage that permits a voltage difference to be maintained between two portions of the electrode. According to the Kuhn-Tucker conditions for Lagrangian problems with inequali- ty constraints, the solution should satisfy ∂L ∂qv i = ∂U ∂qv i − λ − µi = 0, (6) with µi positive when qv i is not negative or µi zero when qv i is negative. In this scheme, µi acts as a switch to maintain the inequality constraint that qv i should be negative, and λ corresponds to the total-charge equality constraint (equation 2). This is a linear minimization problem and can be solved by matrix inversion at each time step throughout the simulation. If there are very many metal atoms, the computational cost of inversion can become excessive. (For a thorough discussion of solving large linear sets of equations see reference [17].) Therefore, we have chosen an alternative implementation. 338 Simulating an electrochemical interface To avoid the CPU-costly matrix inversion at each step, we treat the diffuse charges as dynamic degrees of freedom that respond to forces that decrease the value of L. This is very similar to the way in which the solution particles seek out the energy minima of phase space. We assign each diffuse charge a fictitious mass mq, as well as a charge “velocity” vi. Each qv i is forced down the gradient in the Lagrangian L while the “temperature” of the diffuse charges is kept near zero. In this way the diffuse charges in the electrode are efficiently held near the electrostatic energy minimum at a significant computational savings over matrix inversion. Note that the particle and charge degrees of freedom are coupled as the charge in the metal responds to the structure of the fluid. For this reason separate thermostats must be applied to the particle and charge degrees of freedom in order to maintain them at different temperatures. The diffuse charge and its velocity are propagated using the following equations of motion: q̇v i = vi − ξ, (7a) v̇i = Fi mq − ζvi , (7b) where over-dots indicate time derivatives. ξ, given by ξ = 1 nτ [ n ∑ j=1 (qv j + qc j) − Qtot ] , (8) is an additional parameter to correct long-term numerical drift in the total charge of the electrode. τ is the time constant equal to about 100 timesteps. Long-term drift is possible because it is the derivative of the total charge constraint that is incorporated into the form of the force Fi on each charge. Since we use the Ewald sum to account for long-range Coulombic interactions, it is important to control charges closely in order to maintain cell neutrality. ζ is an integral-feedback control variable that maintains the average charge temperature at a set point T̄q. Using a Nosé-Hoover temperature-control scheme [18], the equation of motion for ζ is ζ̇ = Tq − T̄q T̄qτ 2 , (9) where T̄q is the temperature setpoint – we find that 5 K works well. The instanta- neous charge temperature is Tq = mq nkB n ∑ j=1 v2 j , (10) where kB is Boltzmann’s constant. In fact, the choice of mass mq is tied to the choice of T̄q and the timestep, in order to maintain numerical stability. We use mq = kBT̄q(10 ps/|e|)2. 339 C.G.Guymon et al. The force on each valence charge is Fi = ( 1 n n ∑ j=1 φj ) − φi , (11) where φi, the chemical potential for valence charge i, is φi = ∂U ∂qv i − µi = n ∑ j=1 Cijq v j + n ∑ j=1 C∗ ijq c j + ∂Uext ∂qv i + φset i − µi . (12) Recall that µi is the Lagrange multiplier that acts to constrain qv i to be negative or zero. This inequality constraint acts like a hard wall for charge degrees of freedom. When using continuous dynamics, a softer repulsion is desirable and so we set µi = −kBT̄q Qµ exp ( qv i Qµ ) , (13) which acts as an exponential barrier with Qµ = 0.01 |e| being a stiffness parameter. To summarize, when charged species approach the electrode, the diffuse charges in the electrode adjust to minimize the energy, or to equilibrate the chemical po- tential φi for valence charge at each metal atom. To this point the discussion has focused only on charge-charge interactions that represent gross polarization of the metal conduction-band electrons. In addition to this, electron exchange and corre- lation effects, manifested at short range as a Pauli repulsion and at long range as a dispersion attraction, must be included to properly represent the solution-metal interaction. We therefore incorporate pairwise van der Waals potentials fitted from quantum-mechanical calculations. That procedure is discussed in the next section. Validation of the model is given in section 4 and model predictions are given in section 5. Nevertheless, for efficient use of space some of the figures include at the same time the fit, validation, and predictions of ECD. 3. ECD parameter fit to Cu 10 ab initio data We have applied the ECD method to a copper (111) electrode surface. In this section we report the van der Waals parameters for water, chloride, and sodium ion interacting with a copper surface, in order to supplement the Coulombic interactions from electrode charge dynamics. It should be clearly understood that in the fitting of the van der Waals potentials, Coulombic ECD interactions were not fit. That is, γ and qc are characteristic of the metal and not a function of the approaching species. γ and qc for Cu were chosen to be 0.816 Å−1 and 1 |e| respectively, based on the following rationales. We chose qc = 1 |e| as this results in qv i = −1 |e|, on average, corresponding to 1 conduction electron per copper atom. If the electrode, a face-centered-cubic copper crystal with lattice constant ao = 3.61496 Å, is modeled as a free-electron gas, the Fermi wave vector kF is equal to 1.36 Å−1. At the Fermi level the probability that an available 340 Simulating an electrochemical interface electron energy state is occupied is one-half. If we transform the ECD diffuse charge density (equation 1) to Fourier space and equate the probability that a given wave vector is occupied to the magnitude of the Fourier coefficients, we can relate γ to kF according to 1/2 = exp[−k2 F/(4γ2)]. The resulting value of γ appears to be reasonable: the overall valence charge density undulates smoothly over the electrode and yet large local charge variations are possible due to the fact that the diffuse charge on each copper is confined mostly to the volume within the Weigner-Seitz radius (1.41 Å) [19]. Ab initio scans (energy vs. distance) of molecular and ionic interactions with copper clusters having an exposed (111) plane have been performed to obtain accu- rate Cu-species interactions. The calculations were performed with GAUSSIAN98 [20] at the MP2 level of theory for a ten-copper-atom (Cu10) cluster interacting re- spectively with sodium cation, chloride anion, and water. Details can be found in reference [21] on the water-copper scans as well as the split-valence basis set used for copper for all the calculations here. The basis set 6− 31 + G* was used for both sodium and chloride. On-top, bridge, and hollow routes [21] were probed in the case of chloride, whereas only the on-top site was calculated with sodium. In all of the ab initio results presented here the counter-poise (CP) correction is not included because the basis-set-superposition error (BSSE) is large, particularly for chloride, and the CP method may overcorrect BSSE in the case of metal clusters. Three parameters were regressed for each site (H, O, Cl−, and Na+) interacting with copper atoms by fitting the modified-Morse potential, U(r) = −ε ( 1 − {1 − exp [−A(r − r∗)]}2 ) , to the difference between the ab initio and ECD Coulombic interactions. The SPC/E water geometry and charges [22] are used for ECD calculations and simulations. The resultant modified-Morse parameters and the atomic partial charges are given in table 1. Table 1. Interaction parameters with Cu metal. O H Cl− Na+ ε, kJ/mol 5.000 0.8430 29.42 10.98 A, Å−1 1.350 1.303 1.409 0.9872 r∗, Å 2.890 3.302 2.767 3.741 qsite, |e| –0.8476 0.4238 –1.000 1.000 Shown in figure 2 are the ab initio data points and the ECD results for Cl− approaching a Cu10 cluster over the on-top, bridge, and hollow sites. The sum of the ECD Coulombic and fitted van der Waals potentials agree very well with the ab initio results in both the repulsive and attractive regions (such an agreement is not guaranteed even with three adjustable parameters). The agreement is also good for the calculations of Na+ approaching the on-top site of the Cu10 cluster. 341 C.G.Guymon et al. Figure 2. ECD and ab initio potential energy as a function of distance from the Cu10 (111) plane for Cl−. -50 0 U , k J/ m ol 54321 z, Å -50 0 -50 0 Cu10 ECD Cu10 ab initio Cu18 ECD Cu18 ab initio Cu160 ECD Cu430 ECD (a) hydrogens up, top (b) hydrogens down, top (c) hydrogens parallel, hollow Figure 3. Interaction energy of water with copper clusters of differing size for three unique water approaches. Distances are measured from the oxygen normal to the copper centers in the exposed (111) plane. Figure 3 shows the sum of Coulombic and fitted van der Waals potentials versus the ab initio data for three orientations of water approaching the Cu10 cluster. Nine unique routes were calculated, but only three representative routes are shown. 342 Simulating an electrochemical interface Figure 3 (a), (b), and (c) represent the average, best, and worst fits, respectively, of our model to the water-Cu10 potential scans. (Also shown are validation results with Cu18 clusters and predicted results for larger clusters; these will be discussed shortly.) We suspect a source of the difficulty for ECD here is our use of the rigid (with respect to geometry and dipole moment) SPC/E model of water, for the sake of simplicity. While the fits are not as good with water as with the ions, we still conclude that ECD combined with a van der Waals potential satisfactorily reproduces the ab initio data. Of the nine routes (hydrogens up, down, and parallel approaching the on-top, bridge, and hollow sites), the hydrogens parallel over the on-top site are the most energetically favorable as predicted by the ab initio uncorrected MP2 results [21]. The ECD and van der Waals fit gives the most favorable route as the hydrogens parallel over the bridge site. Both of these results agree with recent DFT calculations on platinum [23] and silver [24] that show the most favorable water orientation to be the hydrogens parallel to the surface. Authors of both studies state that the high-energy p atomic orbital on the oxygen originating from the lone pairs interacts most strongly with the metal surface. The p orbital’s highest density is lateral to the plane containing the oxygen and two hydrogens. 4. ECD validation The ECD model for copper correlates reasonably well the quantum-mechanical interaction potentials of Na+, Cl−, and water with a small Cu10 cluster. In essence we use the ECD method and associated van der Waals potentials as an extrapolation scheme to obtain inexpensive interactions between species and large metal surfaces. We seek to validate and test the ECD copper-surface model by comparing its pre- dictions to image-charge results at long range and to larger cluster calculations. Ab initio data have also been obtained for a Cu18 cluster [25]. We compare the ab initio and ECD predicted energies of chloride, sodium, and water near a Cu18 cluster. We also show ab initio and ECD results of a hydrated sodium atom near the Cu10 and Cu18 (111) cluster surfaces. ECD reproduces the image-charge result at long range and successfully predicts the ab initio energies at the larger cluster size. 4.1. ECD and image charges As stated previously, the simple image-charge potential correctly reproduces the energy of a charge interacting with a metal surface at long range. It is not possible to use the ECD method for a true semi-infinite metal. Therefore we compare the ECD model result with the theoretical (image-charge) result for a point charge approaching a sphere of infinite dielectric constant. In the ECD calculation, a negative and a positive unit charge are separately brought toward a neutral spherical cluster of simulated copper atoms. The distance from the central copper in the spherical cluster to the outer-most copper was less 343 C.G.Guymon et al. than 7 Å with 135 coppers in the cluster. The nearest neighbor distances in the cluster correspond to the bulk fcc lattice, namely ao/ √ 2. On the other hand, the energy of interaction of a point charge q with a neutral sphere of infinite dielectric is U(r) = − q2a3 2r2(r2 − a2) , (14) where a is the sphere radius and r is the distance to the external charge from the center of the sphere. -100 -50 0 U , kJ /m o l 20151050 r - a, Å image pos. charge ECD neg. charge ECD Figure 4. Comparison of ECD with image-charge method for a spherical copper cluster of effective radius a = 7.27 Å. Inset shows negative charge at position r − a = 5.73 Å with cluster atoms shaded according to net charge: black is positive, white is negative, and gray is neutral. Figure 4 shows that the predicted ECD Coulombic energy agrees with the image- charge result at r−a distances greater than about 5 Å. Radius a was adjusted until the ECD and image-charge results agreed at the furthest data point from the sphere, resulting in a = 7.27 Å. This value is consistent with a cutoff radius that could generate the corresponding cluster of 135 discrete atoms. Shown in the figure is the fact that the ECD routine predicts unique energies for the positive and negative unit charges at close range. Inset in the figure is a representation of the copper cluster when a negative charge is at distance r− a = 5.73 Å. Each copper is shaded according to its net charge, qc + qv i , where black is positive, white negative, and gray is neutral. 4.2. Cluster-size results with ECD Figures 3 and 5 show (among other things) the ab initio results for a Cu18 cluster [25] with water, Cl−, and Na+. Due to the large computational expense, only a few points could be calculated. The same level of theory and basis sets as with the Cu10 cluster are used. Distances are measured normal from the surface Cu nuclei to the 344 Simulating an electrochemical interface -300 -200 -100 0 U , k J/ m ol 54321 z, Å -300 -200 -100 0 Cu10 ECD Cu10 ab initio Cu18 ECD Cu18 ab initio Cu160 ECD Cu430 ECD (a) Cl¯ (b) Na + Figure 5. Interaction energy of chloride and sodium with copper clusters of dif- ferent size. The ions approach the on-top site of the (111) plane in all cases. oxygen, chloride, or sodium. The corresponding ECD predictions compare well with the ab initio results in terms of relative displacement of the potential from the Cu10 results. In particular, there is excellent agreement between the Na+ and Cl−–Cu18 ab initio and ECD result. Note that this is a purely predictive result without any re-parameterization. We conclude that the ECD model correctly predicts cluster size dependence. 4.3. ECD and ab initio results for hydrated-sodium clusters near copper We have also compared ECD model predictions to ab initio results for hydrated Na+ near a Cu10 and Cu18 (111) surface, shown in figure 6. Calculations were per- formed for sodium only due to computational ease and because the BSSE is much smaller than with chloride. For the Cu10 cluster, sodium is hydrated with one to five waters. In the case of Cu18, one or two waters surround the sodium. A. Karttunen and T. Pakkanen calculated the optimum geometry of the system at the Hartree-Fock level and the final energy at the MP2 level of theory [25]. In the geometry optimization, all waters were allowed to move freely while Na+ motion was restricted to the direction normal to the surface. Na+ was placed directly over an on-top and a hollow site for the Cu10 and Cu18 clusters respectively. The optimized hydration energy of the isolated sodium-water clusters were also computed. The energies given in figure 6 are the energies of the hydrated ion near a copper cluster, relative to an isolated copper cluster and isolated sodium-water cluster. 345 C.G.Guymon et al. -250 -200 -150 -100 -50 0 U , k J/ m ol 543210 number of hydrating waters Cu10 ab initio Cu10 ECD Cu18 ab initio Cu18 ECD Figure 6. MP2 ab initio optimum and ECD energies as a function of the number of waters near a sodium ion. Geometries were optimized using Hartree-Fock level of theory in the ab initio case. The sodium is above the on-top and hollow site for Cu10 and Cu18 respectively. As shown in the figure, the same calculations were performed using the ECD method, but the oxygen and sodium positions were taken from the ab initio geometry optimization. We used SPC/E [22] water in the ECD calculations and located the SPC/E hydrogens as close as possible to the ab initio hydrogens. For the energy of Na+ interacting with the water we used the parameters of sodium cation found in reference [26]. Although not shown in figure 6, the sodium hydration energy of the ECD and the uncorrected ab initio calculations agree within 4% in all cases. While the energies in the presence of the copper surface for the two techniques compare less well, the ECD model correctly predicts the energetic trend with increasing water. The qualitative change in the ECD curve when 4 and 5 waters are present can be ascribed to the fact that 1 and 2 waters, respectively, are located between the sodium and the metal surface. The ab initio geometries show these waters to be distorted (HOH bisector angle is ∼ 105◦) from the geometry of the other hydrating waters, which more closely resemble SPC/E water. The fact that the ECD model does not reproduce this geometry distortion, in addition to discrepancies in the water-copper interaction are the sources of the inaccuracies in figure 6. 5. ECD predictions In the previous section we showed that ECD does reasonably well in predicting the energies for situations in which it was not parametrized. Here we discuss the ECD predictions of Na+, Cl−, and water on the (111) surface of copper clusters with 160 and 430 atoms-cluster sizes that are far from accessible to fully ab initio calculations. We also show molecular dynamics simulation results using ECD to simulate pure water and an electrolyte solution between two copper electrodes. 346 Simulating an electrochemical interface 5.1. ECD large-cluster predictions Also shown in figures 3 and 5 are the ECD results with clusters of size 160 and 430 atoms for Na+, Cl−, and water. These large clusters are circular-disk in shape and have two metal layers with the (111) plane exposed. At these sizes, the interac- tion energy at the energy minimum is about 60% Coulombic for the ions and about 30% Coulombic for water (hydrogens up, on-top), with the remainder of the energy due to van der Waals site-site interactions. Note that the ECD interaction energy for water converges at a cluster size around 160 atoms while a much larger cluster size is necessary for the ions. ECD predicts the most favorable water orientation to be hydrogens parallel for the Cu10 cluster, but for Cu430 there is slightly greater preference for hydrogens up with oxygen on the bridge site. This prediction is consis- tent with the conventional picture that water chemisorbs to catalytic metal surfaces preferentially through the oxygen. Interestingly, this shift in behavior with cluster size, as predicted by ECD, could help resolve discrepancies between experimental results for extended surfaces and small-cluster ab initio calculations. 5.2. MD simulation details and results The ECD model developed for the copper interface was employed in a series of preliminary molecular dynamics simulations to gain insight into the adsorbed water structure on the (111) copper surface as well as the surface sodium and chloride densities. Figure 7. 3D representation of simulation cell with metal atoms on both sides and a single water molecule, shown to indicate scale, located in the middle of the fluid region. 5.2.1. Simulation details We simulated a sodium-chloride aqueous solution between two copper electrodes at various voltages and two different ionic concentrations. The copper electrodes were each 5 layers thick, with the (111) plane exposed. SPC/E water is used as the solvent and the ion-ion and ion-water potentials are given in reference [26]. The cell dimensions are 97.7, 20.45, and 22.14 Å in the x, y, and z directions respectively, for simulations with no ions and 1200 water molecules. Simulations with ions were 347 C.G.Guymon et al. conducted in a cell of dimensions 101.5, 28.12, 30.99 Å with 80 sodium cations, 80 chloride anions, and 2340 water molecules. Metal atomic positions were fixed during the simulations, with 800 and 1540 total copper atoms for the smaller and larger cell, respectively. Shown in figure 7 is the larger cell geometry with a representative water molecule shown in the middle of the cell. In fact, because of periodic boundary conditions, the two electrodes are joined across the boundary and form one metal slab. The offset potential φset in equation 4 allows us to set the voltage of the two electrode halves at different relative values. In every simulation the total charge on the combined electrode is zero. The total charge in the fluid region is also zero. The charges in the electrodes adjust so as to maintain the voltage difference, as well as to respond to local field changes due to the liquid. Since they constitute constant-potential surfaces, the electrodes behave as a Faraday cage. That is, they shield the usual interaction in a 3-D Ewald sum between periodic images of the unit cell in the direction normal to the electrodes. Effectively, we obtain a 2-D Ewald sum at no additional cost. Here, we use the P3M method due to Hockney and Eastwood [27] to perform the Ewald reciprocal sum. Each simulation run was equilibrated for 50,000 time steps. Data were collected after equilibration for 200,000 or 400 ps for each run unless otherwise stated. With no ions present, 5 runs were simulated with 0, 1 and a 5 volt potential difference between the perfectly polarizable electrodes (i.e., no reactions are permitted) for a total of 15 runs. For simulations with ions, 6 runs were simulated with 0, 0.5 and a 1 volt potential difference between the electrodes for a total of 18 runs. The density of the water in the middle of the solution-filled region is 54 M for electrolyte simulations and 55 M for neat water. Ion concentrations in the middle of the cell are around 1.7 M. -80 -60 -40 -20 0 20 40 60 80 la ye r ch ar ge , µ C /c m 2 108642 layer ∆φ = 0 V ∆φ = 5 V Figure 8. Average charge on the copper layers for two different cases of ∆φ between electrode halves. Layers 1 and 10 are in contact with the surrounding neat liquid water. Error bars show the average of the standard deviation of the charges, with the wider error-bar caps for case ∆φ = 0V. 348 Simulating an electrochemical interface 5.2.2. Simulation results Shown in figure 8 is a graphical representation of the average charge on the electrode atoms as a function of the electrode layer when there is a 0 and 5 volt electrostatic potential difference between electrodes for water only in the liquid. Notice the charge dipole in the metal center that is required to generate the potential difference between the two electrode halves. This charge dipole found at the location of the step change in potential is in agreement with Poisson’s equation (∇2φ = −ρ/ε). When ions are present in the solution with a zero potential drop across the fluid, the charge on the layer is equivalent to the charge without ions but the standard deviation of that charge is twice as large. The large average standard deviations in copper charges for the surface layers indicate that these layers respond most strongly to the fluctuating fields due to the fluid, as we expect. 20 15 10 5 0 de ns ity , m ol /L 806040200 x, Å 40 20 250 200 150 100 50 0 Cl _ Na + oxygen hydrogen (a) (b) Figure 9. Average species density profiles for (a) water sites and (b) ions versus cell position for zero potential difference between the electrodes. Shown in figure 9 is the oxygen and hydrogen density profiles for water with ions present and with a uniform electrostatic potential across the electrodes. When a 1-volt potential difference is placed across the fluid, the resulting changes in the local water density are small and cannot be visibly distinguished. However, the addition of ions to the solution does generate more noticeable differences in the density of the water nearest the surface as the contact-adsorbed ions displace some of the adsorbed water molecules. As for the ion distributions, notable in figure 9 are strongly adsorbed chloride and sodium ions on the surface, in nearly equal numbers. There is a second peak for chloride next to the first, however, with no corresponding sodium peak. This second peak is responsible for the small amount of excess chloride on the surface of the charge-neutral electrode. As with the water density profile, the ion surface density changes little with the imposed voltage between electrodes; most of the change occurs in the diffuse part of the double layer. 349 C.G.Guymon et al. 7x10 -3 6 5 4 3 2 1 0 pr ob ab ili ty , a rb . u ni ts -1.0 -0.5 0.0 0.5 1.0 cosine θ without ions with ions Figure 10. Probability of the cosine of the angle between the surface normal and the surface adsorbed (within 3.5 Å of the surface) water dipole moment. Figure 10 gives the distributions of water dipole orientation in the surface- adsorbed layer (within 3.5 Å of the surface) as indicated by the cosine of the angle between the dipole vector and the surface normal. The electrode potential differ- ence is zero. The waters closest to the surface are strongly oriented with the dipole moment parallel to the surface, with a slight hydrogen-up tendency. Although not shown, this distribution is virtually unchanged under a 5 V range of potential differ- ence between the electrodes. To significantly orient the dipole of the water towards the surface would require a very strong electrode charge. As a consequence, the ECD model predicts a much larger barrier to flipping the surface-water dipole than does prior work with a nonpolarizable charged surface [1,5]. However, the presence of ions in solution is strong enough to create a noticeable disruption to the neat-water dipole distribution. The large surface density for water and narrow dipole distribution indicate that the water is highly structured near the surface. To gain some insight into this phe- nomenon, we rendered a snapshot of the first and second layers of water next to the copper surface within a neat water solution. Figure 11 shows a very ordered two-dimensional rhombus ice-like structure in the first layer. Although other inves- tigators[1,8,28] do not see such a pronounced ordering of the surface layer, X-ray spectroscopy and DFT results [29] suggest that the surface adsorbed water forms “flat ice”. In our model’s ice-like layer, the oxygen-oxygen nearest neighbor distance is about 2% smaller than the bulk liquid value of 2.79 Å. While there is a mismatch between the surface-water lattice and the underlying copper lattice, it appears that the waters somewhat prefer to occupy the bridge and hollow sites. This energetic preference is in agreement with the results obtained by Price and Halley [16]. Figure 12 is a corresponding snapshot with ions present. Only one layer of ad- sorbed water and ions is shown. Sodium cations are shown in black and chloride anions are shown in white. As expected, the ions disrupt the regular array of the water near the surface. This is consistent with the ion-induced changes in figure 10. Somewhat unexpectedly, both ions readily contact the adsorb. 350 Simulating an electrochemical interface Figure 11. Snapshot of the first two water layers on the copper surface for a simulation of neat water. The darker colored waters indicate the layer closest to the surface. Figure 12. Snapshot of the water orientation on the copper surface with ions present. Chloride is shown in white and sodium in black. Shown in figure 13 are the electrostatic potentials as a function of position, with and without ions present when the electrodes are the same potential. The convex 351 C.G.Guymon et al. shape of the electrostatic potential profile in the middle of the cell indicates that there is an excess of sodium ion away from the surface, or a net positive ionic charge in the middle of the fluid (∇2φ = −ρ/ε). Apparently the simulation cell is not large enough to contain the entire diffuse double layer. This is in stark contrast to earlier electrolyte simulations with static surface charges in that the distribution of ions and solvent generates a flat potential profile in the liquid within several Angstroms of the surface [5]. 10 8 6 4 2 0 el ec tr os ta tic p ot en tia l, V 806040200 x, Å NaCl soln. water only Figure 13. Electrostatic potential in the fluid for NaCl solution (solid line) and water only (dotted line). The two systems have different cell lengths as indicated at the right of the plot. The reader may notice that the average ion densities in figure 9 and the elec- trostatic potentials in figure 13 are not symmetric. This is a manifestation of slow relaxation of the double layer that requires longer simulation times than we used in order to generate the expected symmetry in the average properties. Again, the hi- ghly responsive charges in the electrode for the ECD model tend to stabilize charge separation in the double layer to a much greater degree than has been observed in prior simulation work. For example, with an electrolyte between electrodes at the same potential, the average charge-density difference between the two electrode halves averaged to +0.6 µC/cm2 over a 2.4-ns simulation. Since the diffuse region of the double layer is longer than half the length of the fluid-filled region, differential-capacitance calculations of the interface are not reliable. However, we did calculate a cell integral capacitance of 13 and 12 µF/cm2 for the solutions with and without ions, respectively. 352 Simulating an electrochemical interface 6. Conclusion We have developed a simple and robust method to treat the charge mobility in a metal for use in molecular dynamics simulations. Electrode charge dynamics acts as a bridge between small-scale ab initio metal-cluster calculations and large-scale molecular dynamics simulations of metal surfaces of arbitrary geometry. Our use of ECD to simulate a water-NaCl solution near copper electrodes has revealed be- havior with, in some instances, marked differences from prior simulation work. Our simulations show the presence of a dense 2D ice-like rhombus structure of water on the surface that is relatively impervious to perturbation by typical electrode poten- tials/charges. The model also predicts that sodium and chloride ions are strongly adsorbed to the copper surface at both positive and negative electrode charge, but that there is generally an excess of chloride associated with the surface. Based on the potential profile in the cell, it appears that the diffuse part of the double layer ex- tends beyond 40 Å. Further work is required to assess the veracity of the simulations with regard to double-layer properties and to refine the intermolecular potentials. Nevertheless, we believe the approach described here for treating metal polarizability will permit more accurate representations of the electrochemical interface. 7. Acknowledgements We want to thank Dr. Tapani Pakkanen and his students at the University of Joensuu for their work in performing much of the ab initio calculations presented here. We’d also like to acknowledge financial support from the NSF, grant number CTS–0215786 as well as support from Brigham Young University. References 1. Spohr E., Molecular simulation of the electrochemical double layer. Electrochimica Acta, 1999, 44, 1697–1705. 2. Philpott M.R., Glosli J.N. Molecular Dynamics Simulation of Interfacial Electro- chemical Processes: Electric Double Layer Screening, vol. 656(13) of ACS Symposium Series, p. 11–30. Solid-Liquid Electrochemical Interfaces, eds. G. Jerkiewicz, M.P. So- riaga, K. Uosaki, and A. Wieckowski. American Chemical Society, 1997. 3. Yeh In-Chul, Berkowitz M.L., Aqueous solution near charged Ag(111) surfaces: com- parison between a computer simulation and experiment. Chemical Physics Letters, 1999, 301, 81–86. 4. Boda D., Chan K.-Y., Henderson D., Monte Carlo simulation of an ion-dipole mixture as a model of an electrical double layer. Journal of Chemical Physics, 1998, 109, No. 17, 7362–7371. 5. Guymon C.G., Hunsaker M.L., Harb J.N., Henderson D., Rowley R., Effects of sol- vent model flexibility on aqueous electrolyte behavior between electrodes. Journal of Chemical Physics, 2003, 118, No. 22, 10195–10202. 353 C.G.Guymon et al. 6. Crozier P.S., Rowley R.L., Henderson D., Molecular dynamics calculations of the elec- trochemical properties of electrolyte systems between charged electrodes. Journal of Chemical Physics, 2000, 113, No. 20, 9202–9207. 7. Spohr E., Molecular dynamics simulations of water and ion dynamics in the electro- chemical double layer. Solid State Ionics, 2002, 150, 1–12. 8. Ignaczak A., Gomes J.A.N.F., Simulations of liquid water in contact with a Cu(100) surface. Journal of Molecular Structure (Theochem), 1999, 464, 227–238. 9. Evangelakis G.A., Vamvakopoulos E., Pantelios D., Papanicolaou N.I., Coverage de- pendent self-diffusion on Cu(111) by molecular dynamics. Surface Science, 1999, 425, L393–L399. 10. Ursenbach C.P., Calhoun A., Voth G.A., A first-principles simulation of the semicon- ductor/water interface. Journal of Chemical Physics, 1997, 106, No. 7, 2811–2818. 11. Masel R.I. Principles of Adsorption and Reaction on Solid Surfaces. John Wiley & Sons, Inc, 1996. 12. Allen R., Hansen J.-P., Melchionna S., Electrostatic potential inside ionic solutions confined by dielectrics: a variational approach. Physical Chemistry Chemical Physics, 2001, 3, 4177. 13. Boda D., Gillespie D., Nonner W., Henderson D., Eisenberg B., Computing induced charges in inhomogeneous dielectric media: Application in a Monte Carlo simulation of complex ionic systems. Phys. Rev. E, 2004, 69, 046702. 14. Rappé A.K., Goddard III W.A., Charge equilibration for molecular dynamics simula- tions. Journal of Physical Chemistry, 1991, 95, 3358–3363. 15. Rick S.W., Stuart S.J., Berne B.J., Dynamical fluctuating charge force fields: Appli- cation to liquid water. Journal of Chemical Physics, 1994, 101, No. 7, 6141–6156. 16. Price D.L., Halley J.W., Molecular dynamics, density functional theory of the metal- electrolyte interface. Journal of Chemical Physics, 1995, 102, No. 16, 6603–6612. 17. Press W.H., Flannery B.P., Teukolsky S.A., VetterlingJohn W.T. Numerical recipes: the art of scientific computing, vol. 1. Cambridge University Press, 1986. 18. Nose S., A molecular-dynamics method for simulation in the canonical ensemble. Molecular Physics, 1984, 52, No. 2, 255–268. 19. Borg R.J., Dienes G.J. The physical chemistry of solids. Academic Press, San Diego, 1992. 20. Frisch M.J., Trucks G.W., Schlegel H.B., Scuseria G.E., Robb M.A., Cheeseman J.R., Zakrzewski V.G., Montgomery J.A., Jr., Stratmann R.E., Burant J.C., Dapprich S., Millam J.M., Daniels A.D., Kudin K.N., Strain M.C., Farkas O., Tomasi J., Barone V., Cossi M., Cammi R., Mennucci B., Pomelli C., et al. Gaussian 98. Gaussian, Inc. Pittsburgh PA, 2001. Revision A.ll. 21. Ruuska H., Rowley R.L., Pakkanen T.A., MP2 study on water adsorption on cluster models of Cu(111). Journal of Physical Chemistry B, 2004, 108, No. 8, 2614–2619. 22. Berendsen H.J.C., Grigera J.R., Straatsma T.P., The missing term in effective pair potentials. Journal of Physical Chemistry, 1987, 91, 6269–6271. 23. Ohwaki T., Kamegai K., Yamashita K., Electric field effects on the adsorption, charge transfer and vibrational state at metal electrodes: A DFT study on H2O/Pt(111), H2O/Pt(100) and (H2O)2/Pt(111). The Chemical Society of Japan, 2001, 74, 1021– 1029. 24. Sanchez C.G., Molecular reorientation of water adsorbed on charged Ag(111) surfaces. Surface Science, 2003, 527, 1–11. 354 Simulating an electrochemical interface 25. Pakkanen T.A. Private communication, 2004. 26. Wheeler D.R. Molecular Simulations of Diffusion in Electrolytes. PhD thesis, Univer- sity of California, Berkeley, 2002. 27. Hockney R.W., Eastwood J.W. Computer Simulation Using Particles. Institute of Physics, Bristol, 1988. 28. Izvekov S., Voth G.A., Ab initio molecular dynamics simulation of the Ag(111)–water interface. Journal of Chemical Physics, 2001, 115, No. 15, 7196–7206. 29. Ludwig R., How does water bind to metal surfaces: hydrogen atoms up or hydrogen atoms down? Angewandte Chemie-International Edition, 2003, 42, No. 30, 3458–3460. 355 C.G.Guymon et al. Симуляція електрохімічної поверхні з використанням зарядової динаміки К.Дж.Гуймон, Р.Л.Роулей, Дж.Н.Гарб, Д.Р.Віхлер Відділення хімічної інженерії, Прово, США Отримано 29 листопада 2004 р. Ми представляємо простий класичний метод для трактування за- рядової мобільності у металах, які межують з рідкими розчинами. Метод, відомий як електродна зарядова динаміка, ефективно заповнює прогалину між ab initio розрахунками на малих металічних кластерах і велико-масштабними симуляціями металічних повер- хонь з довільною геометрією. Ми отримали модельні параметри для металічної (111) поверхні міді, використовуючи квантово-механічні розрахунки високого порядку на 10-атомному мідному кластері. Модель була перевірена шляхом порівняння з класичними резуль- татами по методу відображень, а також з ab initio результатами для 18-атомного мідного кластера. Модель використана в молекуляр- но динамічних розрахунках для структури флюїдної поверхні чистої води і водного розчину NaCl. Ми побачили, що вода організовується у двовимірний льодоподібний шар на поверхні і що обидва, Na + i Cl − , іони є сильно прив’язані до міді. Коли металічний електрод за- ряджати, то основна реакція електроліту проявляється у дифузивній частині подвійного шару. Ключові слова: симуляція, подвійний шар, молекулярна динаміка, ab initio потенціали, (111) поверхня міді, вода PACS: 61.20.Qg, 61.20.Ja 356