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 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.


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 interface 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.Simulations have been used to predict the adsorbed water orientation, ionic distributions, ionic transport coefficients, atomically resolved electrostatic potential profiles, and interfacial differential capacitance [1][2][3][4][5][6].
Accurate molecular simulations of an electrochemical interface depend on the quality of intermolecular potentials between species.Prior work has resulted in ad-equate interaction potentials between species in bulk liquid electrolytes, however the potentials suitable at a metallic interface are less well established.For an intermolecular 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][8][9][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 image 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 nonuniformities.For example, the image-charge method predicts the same interaction energy for positive or negative charges, and does not differentiate between surface adsorption 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-chargetype 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 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.

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 q v i .A fixed positive point charge q c i is co-located at each center to represent the nucleus and core electrons.Each diffuse charge is modeled with a Gaussian charge-density distribution, 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 q v i changes so as to minimize the electrostatic energy or equivalently to equilibrate the chemical potential for q v i at each atom.The fluctuation in q v i is subject to two constraints.The first constraint is that q v i 0, corresponding to a limit on the depletion of electron density from a given atom.The second constraint is a fixed total charge Q tot in the simulated metal: where the sum is over all metal atoms.
At this point we note the similarity between our model and the charge equilibration 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 q c and γ) that characterize the polarization 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 interactions.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 minimize 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 where λ and µ i are undetermined multipliers, and C ij and C * ij are Coulomb overlap integrals given by with . U ext 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 inequality constraints, the solution should satisfy with µ i positive when q v i is not negative or µ i zero when q v i is negative.In this scheme, µ i acts as a switch to maintain the inequality constraint that q v 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.
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 m q , as well as a charge "velocity" v i .Each q v 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: where over-dots indicate time derivatives.ξ, given by 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 F i 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 Tq .Using a Nosé-Hoover temperature-control scheme [18], the equation of motion for ζ is where Tq is the temperature setpoint -we find that 5 K works well.The instantaneous charge temperature is where k B is Boltzmann's constant.In fact, the choice of mass m q is tied to the choice of Tq and the timestep, in order to maintain numerical stability.We use The force on each valence charge is where φ i , the chemical potential for valence charge i, is Recall that µ i is the Lagrange multiplier that acts to constrain q v 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 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 potential φ 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 correlation 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.

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 q c are characteristic of the metal and not a function of the approaching species.γ and q c for Cu were chosen to be 0.816 Å−1 and 1 |e| respectively, based on the following rationales.We chose q c = 1 |e| as this results in q v 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 a o = 3.61496 Å, is modeled as a free-electron gas, the Fermi wave vector k F is equal to 1.36 Å−1 .At the Fermi level the probability that an available 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 k F according to 1/2 = exp[−k 2 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 accurate Cu-species interactions.The calculations were performed with GAUSSIAN98 [20] at the MP2 level of theory for a ten-copper-atom (Cu 10 ) cluster interacting respectively 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, 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.  2 are the ab initio data points and the ECD results for Cl − approaching a Cu 10 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 Cu 10 cluster.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.

ECD validation
The ECD model for copper correlates reasonably well the quantum-mechanical interaction potentials of Na + , Cl − , and water with a small Cu 10 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 predictions to image-charge results at long range and to larger cluster calculations.
Ab initio data have also been obtained for a Cu 18 cluster [25].We compare the ab initio and ECD predicted energies of chloride, sodium, and water near a Cu 18 cluster.We also show ab initio and ECD results of a hydrated sodium atom near the Cu 10 and Cu 18 (111) cluster surfaces.ECD reproduces the image-charge result at long range and successfully predicts the ab initio energies at the larger cluster size.

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 than 7 Å with 135 coppers in the cluster.The nearest neighbor distances in the cluster correspond to the bulk fcc lattice, namely a o / √ 2. On the other hand, the energy of interaction of a point charge q with a neutral sphere of infinite dielectric is where a is the sphere radius and r is the distance to the external charge from the center of the sphere.Figure 4 shows that the predicted ECD Coulombic energy agrees with the imagecharge 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, q c + q v i , where black is positive, white negative, and gray is neutral.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 Cu 10 results.In particular, there is excellent agreement between the Na + and Cl − -Cu 18 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.

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 Cu 10 and Cu 18 (111) surface, shown in figure 6. Calculations were performed for sodium only due to computational ease and because the BSSE is much smaller than with chloride.
For the Cu 10 cluster, sodium is hydrated with one to five waters.In the case of Cu 18 , 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 Cu 10 and Cu 18 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.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.

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.

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 interaction 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 Cu 10 cluster, but for Cu 430 there is slightly greater preference for hydrogens up with oxygen on the bridge site.This prediction is consistent 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.

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.

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 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.

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.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.Figure 10 gives the distributions of water dipole orientation in the surfaceadsorbed 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 difference 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 difference between the 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 phenomenon, 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 investigators [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 adsorbed 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.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 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].The reader may notice that the average ion densities in figure 9 and the electrostatic 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 highly 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/cm 2 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/cm 2 for the solutions with and without ions, respectively.

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 behavior 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 potentials/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 extends 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.

Figure 1 .
Figure 1.Charge on each metal atom.

Figure 2 .Figure 3 .
Figure 2. ECD and ab initio potential energy as a function of distance from the Cu 10 (111) plane for Cl − .

Figure 3
Figure3shows the sum of Coulombic and fitted van der Waals potentials versus the ab initio data for three orientations of water approaching the Cu 10 cluster.Nine unique routes were calculated, but only three representative routes are shown.

Figure 3 (
Figure3(a), (b), and (c) represent the average, best, and worst fits, respectively, of our model to the water-Cu 10 potential scans.(Also shown are validation results with Cu 18 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.

Figure 4 .
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.

Figures 3 Figure 5 .
Figures3 and 5show (among other things) the ab initio results for a Cu 18 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 Cu 10 cluster are used.Distances are measured normal from the surface Cu nuclei to the

Figure 6 .
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 Cu 10 and Cu 18 respectively.

Figure 7 .
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.

Figure 8 .
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 ∆φ = 0 V.

Figure 9 .
Figure 9. Average species density profiles for (a) water sites and (b) ions versus cell position for zero potential difference between the electrodes.

Figure 10 .
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 11 .
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 .
Figure 12.Snapshot of the water orientation on the copper surface with ions present.Chloride is shown in white and sodium in black.

Figure 13 .
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.

Table 1 .
Interaction parameters with Cu metal.