Applied Potentials in Variable Charge Reactive Force Fields for

Subscriber access provided by University of Florida | Smathers Libraries

Article

Applied Potentials in Variable Charge Reactive Force Fields for Electrochemical Systems Tao Liang, Andrew C Antony, Sneha A. Akhade, Michael J. Janik, and Susan B. Sinnott J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b06064 • Publication Date (Web): 19 Dec 2017 Downloaded from http://pubs.acs.org on December 21, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Applied Potentials in Variable Charge Reactive Force Fields for Electrochemical Systems Tao Liang1,*, Andrew C. Antony1,2, Sneha A. Akhade3,4, Michael J. Janik3, and Susan B. Sinnott1,* 1

Department \of Materials Science and Engineering, The Pennsylvania State University, University Park, PA, 16802, USA

2

Department of Materials Science and Engineering, The University of Florida, Gainesville, FL, 32611, USA

3

Department of Chemical Engineering, The Pennsylvania State University, University Park, PA, 16802, USA

4

Institute for Integrated Catalysis, Pacific Northwest National Laboratory, Richland, WA, 99352 USA

Abstract An atomic description of water dynamics and electrochemical properties at electrode-electrolyte interfaces is presented using molecular dynamics with the third generation of the charge optimized many body (COMB3) potential framework. Externally applied potentials in electrochemical applications were simulated by offsetting electronegativity on electrode atoms. This approach is incorporated into the variable charge scheme within COMB3 and is used to investigate electrochemical systems consisting of two Cu electrodes and a water electrolyte with varying concentrations of hydroxyls (OH-) and protons (H+). The interactions between the electronegativity offset method and the charge equilibration method in a variable charge scheme are analyzed. In addition, a charge equilibration method for electrochemical applications is proposed where the externally applied potentials are treated by the electronegativity offset on the electrodes thus enforcing charge neutrality on the electrolyte. This method is able to qualitatively capture the relevant electrochemistry and predict consistently correct voltages with pre-calibration.

*

corresponding authors: [email protected] and [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction Atomic-scale modeling of metal-liquid interfaces is of interest in numerous applications, including those that involve catalytic, corrosive, and electrochemical systems. Of the available atomistic modeling techniques, molecular dynamics (MD) simulations with empirical interatomic force fields, or potentials, are attractive for their ability to provide atomic-level, mechanistic information about dynamic processes within a system at a length scale that is beyond the typical accessibility of quantum mechanical techniques. A number of nonreactive, fixed charge water force fields, including TIP and the SPC series [1], have been successfully developed and used to predict the properties of ions in aqueous solution [2-4], water adsorption on metal surfaces [5-7] and metal-liquid interfaces under externally applied voltages [8-11]. More recently, the bond order concept and self-consistent charge equilibration scheme have been introduced to interatomic force fields [12, 13]. An interatomic force field with the bond order concept is often referred to as a reactive force field, or potential, as it allows bond formation or dissociation during simulations. The self-consistent charge equilibration scheme, on the other hand, is referred as a variable charge model since it dynamically determines atomic charge according to an atom’s local electrostatic environment. Several force fields that incorporate the bond order concept and/or a variable charge equilibration scheme have been developed either exclusively for water [14-16] or for a wide of range of materials. Some of them have been used to explore the surface chemistry at metal-liquid interfaces [17-19]. Among the reactive variable charge potentials, ReaxFF [20, 21] and chargeoptimized many body (COMB) [13, 22], are widely used for multicomponent systems. The influence of externally applied potentials on the electrode-electrolyte interfaces has been previously studied using MD simulations. Several studies [23-26] have examined the properties of the electrolyte between charged electrodes, in which a constant electrode charge is assumed. While giving new insights at atomic-level, this constant electrode charge method has limitations to represent the induced charge on electrodes. Rather a constant charge induction, the induced charge on real electrodes is dynamically changing to maintain a constant potential at electrodes. Siepmann and Sprik [27] developed a model to simulate the constant potential at metallic electrodes. In this method, the atomic charge on each electrode atoms is expressed as a Gaussian distribution function. The electrostatic potential acting on any atom, which is equal to the derivative of the induced electrostatic energy with respect to the atomic charge, is a constant. The induced charge density distribution on electrodes is satisfied the Poisson’s equation and thus the overall electrostatic potential drop at electrodes is equal to the pre-set value V. By parameterizing the width of the Gaussian charge distribution of different metals, this method is able to reproduce the induced electrostatic potential on different metallic electrodes in electrochemical applications. Based on the work of Siepmann and Sprik, several theoretical methods [9, 28-30] have been developed to simulate the metallic electrodes in classic dynamics. For example, Reed et al. [28] and Limmer et al.[29] developed an algorithm to simulate the metallic electrodes, in which the induced charge on each atom is adjusted on-the-fly to maintain the constant potential drop on electrodes. Onofrio and Strachan introduced the electric resistivity in the electrochemical dynamics with implicit degrees of freedom (EChemDID) [30] method, which can be used to estimate electronic current in electrochemical processes. Guymon et al. [9]


Page 2 of 17

Page 3 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


developed the electrode charge dynamics (ECD) method, in which the electrode atom is assumed to consist of a nucleus with a unit point charge and valence electrons with the Gaussian distribution around its nucleus. The externally applied potential in ECD method is simulated as an add-on electrostatic force field on the electrodes in a very similar manner to Siepmann and Sprik’s method. Yeh et al. [10] used the ECD method to investigate the electrode-electrolyte interface associated with water at Pt (111) electrodes, where the liquid water is described by nonreactive SPC/E [31] and reactive modified central force (mCF) [32] fixed charge force fields. The results indicate that the interfacial configuration as well as water dynamics has a strong dependence on the applied potentials. However, because both SPC/E and mCF are fixed charge force fields, they are unable to model local charge transfer at the electrode-electrolyte interface. In addition, parameterization was required since no Pt-H2O interactions were included in the SPC/E and mCF force fields. To more realistically simulate electrochemical surface reactions that include proton transfer with solution, a reactive potential with a variable charge scheme is desirable. In a force field with the variable charge scheme, the partial derivative of the electrostatic energy over the atomic charge of the ith atom is defined as the electronegativity of the atom [33, 34]. Therefore, the induced electrostatic potential on each electrode atom can be regarded as the electronegativity difference of the atom, which is acting as an extra force on the atomic charge. If the induced electronegativity of any electrode atom is a constant, the concept of a constant electrostatic drop at electrodes can be easily achieved by offsetting the electronegativity at the electrodes in the variable charge reactive force fields [35, 36]. For example, Onofrio et al [35] investigated the electrometallization of Cu on amorphous SiO2 using the ReaxFF reactive, dynamic charge force fields. To realize the applied potentials on electrodes, they modified the electronegativity parameter within the force field for the electrode atoms. Here, we present the electronegativity offset method to simulate the effects of explicit applied potentials on charge transfer in combination with the third-generation charge optimized manybody (COMB3) potential [22, 37], in which the atomic charge is determined by dynamic charge equilibration (QEq) scheme that was proposed by Rick et al. [38]. No change to existing COMB3 parameters is required and there is no need for additional parameterization of the COMB3 potential. This electronegativity offset method is then used to investigate water dynamics and electrochemical reactions within Cu/water/Cu electrochemical systems. The results are compared to the findings of previous models from the literature. Due to the limitation of the variable scheme in such potentials, the results indicate that it is best to apply the charge equilibration (QEq) separately on the electrodes and electrolyte and use the electronegativity offset method to account for the applied potentials in electrochemical systems. 2. Computational Details 2.1 Methodology As shown in Eq. (1), the total energy (Utot) of the system in COMB3 consists of the self-energy term (Uself), the Coulomb interactions (UCoul), short-range bond energy (Ubond), van der Waals interactions (UvdW) and several other ad hoc terms (Uothers).



Utot = Uself + UCoul + UvdW + Ubond + Uothers

(1)

The first two terms in Eq. (1) are the electrostatic energy in COMB3. In particular, the self energy of the system equals to the summation of the self-energy of an ion ( ), which is approximated as second order of Taylor expansion series with respected to its charge (qi) [39]:

= ∑

= ∑( + )

.

(2)

In Eq. (2), the parameters and are the electronegativity and the chemical hardness of an element, which are fitted to the electron affinity and the ionization energy of a single isolated atom. COMB3 makes use of the charge density distribution functions proposed by Streitz and Mintmire [40], in which the charge density of an atom is a function of its charge, the spatial location (r) and atomic position (ri): ⍴ () = δ(| − |) + ( − ) (| − |), (| − |) = exp #−2 | − |%

.

(3) (4)

Here, Zi is an effective point core charge treated as a fitted parameter, δ is the Dirac delta function, f(r) is a function that captures the radial decay of the electron density of the S-type orbital, and ξi is an orbital exponent that controls the length scale associated with this decay. Thus, the Coulomb term (UCoul) includes charge-charge and nuclear charge-charge interactions between each ion pair. The Wolf summation [41] with a cutoff radius of 11 Å is employed to compute the Coulomb interactions in COMB3. The details of COMB3 potentials are reviewed in Ref. [22]. In the variable charge scheme, the electronegativity (χi) of atom i is the partial derivative of total energy (Utot) with respect to its charge (qi). 1 ∂U tot and χ = ∑ χ i χi = (5) ∂qi N i In Eq. (5), ̅ is the mean value of electronegativities of the N-atom system. The atomic charge of each atom in COMB3 is determined by dynamic QEq scheme, in which the equation of motion for the charge of the ith atom is expressed as:

M Q q&&i = χ − χ i + υq&i ,

(6)

where )* is a uniform fictitious mass of charge and υ is a damping factor on the velocity of the atomic charge. The quantity ̅ − is the force acting on the charge of the ith atom. The electronegativity equalization (EE) condition is achieved by iteratively solving Eq. (6) for all atoms until the electronegativity of every atom is equal to ̅ . To ensure charge conservation, COMB3 computes the summed value of atomic charges or net charge (Q) of the system at any moment and compares it with the value at time zero. If the difference between these two


Page 4 of 17

Page 5 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


numbers (∆Q) is larger than 0.001 e, COMB3 evenly shifts the atomic charge on each atom by ∆Q/N to ensure charge conservation, where N is the number of atoms of the system. Both ReaxFF and COMB3 lack electronic structure degree of freedom; therefore, neither can simulate electrical conductivity or rate of charge transfer within a system. Both assume that charges are equilibrated before ionic motion takes place, which implies that the material under QEq is an instantaneous conductor. Before introducing externally applied voltages, we first performed MD simulation for 2000 ps for the electrochemical system in Fig. 1, with no external applied voltage. The dynamic QEq scheme is applied to the whole system using the scheme expressed in Eq. 6. The resulting water dynamics are provided in the Supplementary Information. The current variable charge scheme in COMB3 leads to a significant amount of net charge on water electrolyte (-0.039 e per water molecule in Fig. 1S), which is unrealistic for the system under consideration. As pointed out by Chen and Martinez [42], this unrealistic charge transfer is primarily due to Mulliken electronegativity, χ0 in Eq. (2), which is described as a distance-independent parameter in QEq model. This distance-independent feature of Mulliken electronegativity has led to this unrealistic charge transfer as discussed in the Supplementary Information. It is therefore desirable to use separate QEq schemes in hetero-junctions consisting of materials with substantial differences in electronegativities with only local charge transfer enabled. This separate QEq scheme underestimates the local charge transfer at the electrode and electrolyte interfaces. The focus of this work is on how to apply the electronegativity offset method on electrodes within the variable charge scheme to account for the externally applied potentials in electrochemical systems. Compared with the constant potential model [27, 28] for metallic electrodes with Eq. (5), the induced electrostatic potential on electrodes has the same quantity with the electronegativity drop in QEq scheme. If the induced potential experienced by atomic charge on each electrode atom is a constant v0, the equations of motion of charges of the electrolyte and electrodes can be expressed in Eq. (7). The charge neutrality is enforced on the electrodes and electrolyte, respectively. : )* 0 = ̅1 − + 2 3 and ∑ 1 = 0 . ,8: )* 0 = ̅9: + ;< v − + 2 3 ;= >;< ;=

-@: )* 0 = ̅9: − v − + 2 3 ;= >;< , and ∑ 9: = 0 +

(7)

Here, ̅1 and ̅9: are the mean value of electronegativities of the electrolyte and electrodes (cathode and anode). 1 and 9: are the atomic charge of the electrolyte and electrodes (cathode and anode), respectively. A9 and A: are the number of atoms at cathode and anode. If A9 = A: , the electronegativity offset on each electrode atom is ± v0/2. In the case where A9 ≠ A: , the coefficients ahead of v0 in Eq. (7) are applied to conserve the total electronegativity of the system.

The magnitude of v0 is depending on the decay factor for the charge density function, i.e. ξi in Eq. (4) and Gaussian width parameter in Ref. [27, 28], and the dimensions of the electrodes



for a given system. Similar to the previous work [9, 27, 28], v0 is a preset number in this work to realize the overall induced electrode potential is the applied potential V. In particular, the value of v0 is set to be 0.792/nL per volt, where nL is the number of copper layers of each electrode (Fig. 1). The calibration of v0 in this paper is presented in Supplementary Information.

2.2 Simulation Setup The electrochemical cell in the simulations consists of two parallel 6-layer Cu (111) slabs (each layer has 168 atoms and each slab has 1008 atoms) as the electrodes and 2240 water molecules as the electrolyte (Fig. 1). The dimensions of the electrochemical system are 90.1 × 31.0 × 30.7 Å in length (parallel to the surface normal of both electrodes), width and height, respectively. The density of the liquid water is 1.0 g cm-3, which is the value predicted by COMB3 for bulk liquid water [17]. A vacuum of width 20 Å is added along the surface normal to avoid interactions between electrodes and periodic boundary conditions are applied in all three directions.

Figure 1. System configuration of Cu/water/Cu electrochemical cell. Yellow atoms are Cu, red atoms are O and white atoms are H. The letters C, E and A represent the cathode, electrolyte and anode, respectively. The letters b and s represents the bulk and surface, respectively. In order to investigate the influence of applied external potentials on the transport of hydroxyl/proton ions, randomly distributed hydroxyl ions (OH-) and protons (H+) pairs are introduced to the water. The number of introduced OH-/H+ ions is estimated from the capacitance of electrochemical double-layer capacitors (EDLCs). For convenience, the induced charge on the electrode is selected to be 10 e⸱V-1. Since the surface area of the electrodes is 951.7 Å2, the selected value of 10 e⸱V-1 corresponds to a capacitance value of 16.9 µF⸱cm-2 for EDLCs, which is within the scope of estimated values (5 to 40 µF⸱cm-2) [43, 44]. If the countercharges in EDLCs that balance the net charge on the electrodes come from OH-/H+ ions, we have thus introduced four, seven and thirteen OH-/H+ ions for the applied potentials of 0.625 (±3.13 e), 1.25 (±6.25 e) and 2.50 V (±12.50 e), respectively. These concentrations of ions are higher than average values in liquid water, which are around 1.0×10-7 mol⸱L-1. For example, four pairs of ions lead to an ionic concentration of about 6.4×10-2 mol⸱L-1. This selected concentration is the smallest concentration that allows for reasonable sampling of ion behavior,


Page 6 of 17

Page 7 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and helps accelerate the generation of ions that would otherwise proceed slowly even with applied potentials in a finite system. Two systems (with or without OH-/H+ ions) were subjected to MD simulations, where the QEq schemes were applied to electrodes and electrolyte separately as indicated in Eq. (7). Since nL equals to 6 for these two systems of interest, the offset electronegativities (v0 in Eq. (7)) on each electrode atom are 0.165, 0.33, and 0.495 V to represent the applied potentials of 0.625, 1.25 and 2.5 V between cathode and anode, respectively. In all the MD simulations, the interatomic forces are determined using COMB3 potentials [37] as implemented within the LAMMPS software [45]. Throughout the simulations, the Cs and As layers in Fig. 1 are kept fixed. The rest of the system is considered within the NVT ensemble, in which the thermostat used is the Nose-Hoover algorithm [46, 47] on Cu atoms and Langevin [48] algorithm on water molecules to maintain the temperature of the system at 300 K. The damping parameter in both thermostat algorithms is set to 1 ps and the time step is set to 0.00025 ps or 0.25 fs. Each MD simulation for all cases is allowed to run for 2000 ps, which is about double the typical relaxation time [10] for a dipole moment in bulk water. 3. Results and Discussion After evolving the electrochemical system in Fig. 1 within an NVT ensemble for 2000 ps, we analyzed the atomic charge distribution, density profile of water, water dipole moment, and calculated potential. The data presented for these properties of water dynamics are the time averaged values of the last 20 ps of the MD simulation.



Page 8 of 17

Figure 2. Average charge per atom of different regions in Figure 1 for a) systems without OH-/H+ ions; b) for systems with OH-/H+ ions; and c) total induced charge on the cathode. To better quantify the charge distribution, the system is divided into seven regions as indicated in Fig. 1, where Eb represents the bulk water region. The thickness of the interfacial electrolyte (EC and EA in Fig. 1) is set to 10 Å from the cathode and anode surfaces, where EC and EA combined takes up about 30% of the total electrolyte. For the case of no applied external potential (Fig. 2a), the average atomic charge is about 0.007 e for surface Cu atoms (Cs and As), -0.003 e for bulk Cu atoms (Cb and Ab), 0.003 e for Cu atoms at cathode interfaces (CE), and 0.007 e for Cu atoms at anode interface (AE). Since the local electrostatic environment of Cu atoms is different, Cu atoms at surfaces and interfaces have a small amount of positive charge. As charge neutrality is enforced, Cu atoms at bulk are thus slightly negatively charged. As the applied potential increases to 2.5 V, the average atomic charge increases to 0.010 e for cathode Cu atoms (CE) and decreases to 0.002 e for anode Cu atoms (AE) for the systems with or without ions. Based on Fig. 2a and 2b, the total induced charge (Qind) on each electrode by applied potentials can be computed using the following equation ^ _` = ∑ a − ∑ ,

(8)

where qiV and qi0 are the atomic charge of ith atom with the applied potential V and without the applied potential. Since the countercharges are expected on anode, we only plot the induced charges on atoms at cathode in surface in Fig. 2c, in which black diamonds, blue squares and red triangles represent the induced charge for Cs (168 atoms), Cb (672 atoms) and CE (168 atoms), respectively. The systems without ions are plotted in solid lines and with ions are plotted in dotted lines in Fig. 2c. There is no significant difference in induced charges on Cu atoms for both systems with or without ions. For both systems, the induced charges on the cathode increase linearly with the applied potential for Cs, Cb and CE atoms. As shown in Fig. 2c, the induced


Page 9 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


charges by 2.5 V of applied potential are 0.004, 0.001 and 0.008 e per atom for Cs, Cb and CE atoms, respectively. The summed induced charges on these three types of atoms are about 0.6 e, 0.6 e and 1.4 e, respectively. Compared with the charge fluctuation model in Ref. [28, 29], in which the induced charges are almost 100% accumulated on CE, the induced charge on CE only takes 54% of the total induced charge in this work. This underestimation is primarily due to the separated QEq scheme in Eq. (7), which prevents the local charge transfer at interfaces. The induced charge per CE atom in this work has the similar value with the work of Ref. [28, 29]. It is worthwhile to point out that the induced charges on anode atoms are not fully symmetric with cathode atoms. In particular, the induced charge per AE atom is -0.005 e at V = 2.5 V. Since atomic charge of AE atoms at V=0 is about 0.007 e, the atomic charge of AE atoms at V=2.5 V is still positively charged (0.002 e in Fig. 2a and 2b) for both cases.

Figure 3. Water density profile as a function of applied potentials. The legend is shown at the bottom of the figure. The density profiles of water as a function of the applied potentials for the two types of systems are provided in Fig. 3, where the horizontal line at 1.00 g⸱cm-3 indicates the COMB3 predicted density of bulk liquid water at 300 K. Since the density profile of water is dependent on water structure near the electrode surfaces. The combined density profile and water structure near the cathode surface for the system without ions at V=0 is shown in Fig. 4, where the O and H atoms between first and second peaks on the density profile are highlighted in blue and gray, respectively.



Figure 4. The water structure near cathode surface with the density profile (black line) for the system without ions at V=0. Blue atoms are O and gray atoms are H. For all the cases, the density profiles in Fig. 3 have the similar pattern in which the first peak (4.2 to 5.1 g⸱cm-3) occurs at 2.4 Å away from the electrode surfaces and the second peak (1.6 to 1.9 g⸱cm-3) at 5.2 Å away from the electrode surfaces. The separation between the first and second peaks is about 2.8 Å, which is the O-O distance in water. As indicated in Fig. 4, the negatively charged O atoms orientate towards the positively charged cathode surface, which leads to the much higher water density (4.5 g⸱cm-3 in Fig. 3) at the cathode/water interface. Since the tetrahedrally coordinated hydrogen bonds are favorable in liquid water, the density of water in between the two peaks (blue and gray highlighted atoms in Fig. 4) decreases to 0.5 g⸱cm-3, which maximize the four hydrogen bond coordinated H2O molecules near the interface. With increased potentials, the H atoms are expected to orientate towards the negatively charged Cu atoms at the anode, which leads to the smaller and left-shifted intensification as indicated in the bottom right image in Fig. 3. Within this study, the averaged atomic charge on anode atoms is small and remains positive under an applied potential of 2.5 V (Fig. 2). As a result, the applied potentials have no significant effects on the density profiles.


Page 10 of 17

Page 11 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The average water dipole along the Z direction (µ(Z)). Upper figures are the systems without ions and the lower figures are the systems with ions. The images on the left are for water molecules near the cathode and the images on the right are for water molecules near the anode. The legends are the same with those in Figure 3. Fig. 5 illustrates the Z component of the water dipole (µ(Z)), where Z is the direction from the cathode to the anode. The averaged magnitude of dipole moments in liquid water from COMB3 [17] is about 0.45 e⸱Å with the direction of positive charge center aimed at the O atom in a water molecule. If the water dipole is normal to the electrode surface with the O atom towards Cu, µ(Z) is about -0.45 e⸱Å at the cathode and 0.45 e⸱Å at the anode. For a randomly orientated water molecule, the time average of dipole moments of water molecules is zero. As applied voltage increases, µ(Z) decreases near the cathode which indicates that more O-down water molecules orientate towards Cu atoms. Near the anode area, the values of µ(Z) decrease but remain positive as the applied voltage increases, which suggests that the main orientation of the water molecules is still O-down towards anode surface. For all the cases, the minimum µ(Z) near the cathode or the maximum µ(Z) near the anode occurs at about 4.0 Å from electrode surfaces, which represents the higher ordering of the water molecules between the two peaks in the density profiles (the blue and gray highlighted water molecules in Fig. 4). Due to the intensification and orientation of water molecules at the electrode interfaces, the H atoms from the first water layer form a positively charged layer. To maximum the hydrogen bonds, the oxygen atoms that are highlighted in blue color in Fig. 4 orientate towards these H atoms from the first water layer. As a result, the value of µ(Z) reaches the minimum at the cathode and the maximum at the anode.



To visualize the dynamic of the water ions (OH- and H+), three snapshots of the systems with an applied external potential of 0.625 V (4 ion pairs), 1.25 V (7 ion pairs) and 2.50 V (13 ion pairs) are shown in Fig. 6a to 6c, respectively. The averaged displacement along the Z direction (D(z)) of OH- and H+ ions for three simulated voltages are provided in Fig. S2.

Figure 6. Snapshots of the system after 2000 ps with applied potentials of a) 0.625 V (4 pairs of OH-/H+ ions), b) 1.25 V (7 pairs of OH-/H+ ions), and c) 2.5 V (13 pairs of OH-/H+ ions). The atoms are colored by type, light blue is O and yellow is H in the OH- ions and red atoms are H+ ions. The water molecules are not shown for clarity. The super script * indicates the water ions are absorbed on an electrode surface. The OH- and H+ ions form H2O and H2 molecules over the course of the MD simulations as shown in Figs. 6a to 6c. One, three and six H2O molecules formed from OH- and H+ ions under external applied potentials of 0.625, 1.25 and 2.50 V, respectively. One H2 molecule is formed in the cases of 1.25 and 2.5 V applied potentials. Recombination of OH- and H+ ions is physically reasonable since the concentration of water ions is high and the formation of water and hydrogen gas is energetically favorable in COMB3. Due to the positively charged anode atoms (0.002 e in Fig. 2) at 2.5 V, two OH- ions are absorbed on the anode surface rather than move towards the cathode surface.


Page 12 of 17

Page 13 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Using the charge density profile, we can calculate the potential, V, via Poisson’s equation:

b c = −

⍴(d) eef

,

(9)

where ⍴(z) is the charge density profile along the length direction, which is obtained from summed atomic charges divided by the unit volume. The potential profile or electronegativity drop with respect to the anode is presented in Fig.7, where the vertical lines indicate the cathode and anode surfaces.

Figure 7. The potential profile of the systems a) without water ions and b) with water ions. The vertical lines indicate the positions of cathode and anode surfaces. The legends are the same as in Fig. 3. In the case of the systems without water ions (Fig. 7a), the differences of V between the cathode and anode surface, i.e. ∆V = VC - VA, are about 0.07, 0.63, 1.45 and 2.50 V for the preset voltages of 0.0 (black solid line), 0.625 (purple dotted dash line), 1.25 (blue dash line) and 2.5 (red dotted line), respectively. In the case of the systems with water ions (Fig. 7b), the differences of V between the cathode and anode surfaces are about 0.74, 1.25 and 2.63 V corresponding to V = 0.625, 1.25 and 2.5 voltages, respectively. With calibration of v0, the calculated potentials are consistent with the pre-applied potentials. In summary, the electronegativity offset method depicted in Eq. 7 is able to simulate the externally applied potential in the variable charge force fields for electrochemical application. The work presented in this paper suggests that using the electronegativity offset method in COMB3 is able to qualitatively capture the influence of applied potentials on electrochemical reactions and water dynamics in electrochemical applications. It is worthwhile to point out that all the methods including ECD [9] and EChemDiD [30] need to pre-define the electrodes and electrolyte before the simulation. This prohibits the methods to simulate the phenomenon such as sacrificial anode.

4. Conclusions The developments outlined within this paper are applying the classic constant electrode model within the QEq charge equilibration scheme in an empirical variable charge force field to



simulate externally applied potentials in electrochemical applications. The implementation of this method inside QEq scheme is simple and straightforward: offsetting electronegativity on electrodes to account for electrode potential induced by the externally applied potential. To prevent from unrealistic net charge accumulation on electrolyte, which is one of limitation of variable charge schemes in current empirical force fields, the charge neutrality is therefore enforced on electrolyte in this study. We apply the electronegativity offset method to MD simulations of electrochemical cells consisting of Cu (electrodes) and water (electrolyte) using COMB3 force field to study the effects of applied potentials on water dynamics. For the systems with or without OH- and H+ ions, the induced atomic charge of cathode Cu atoms increases about 0.003 eV-1. Such charge transfer effects determine electrochemical characteristics such as overall electrolyte density, dipole moment and displacement of water ions. In contrast, in the case of the system with OH- and H+ ions, we predict several chemical reactions including formation of water molecule, hydrogen gas and OH- and H+ co-absorption on the electrode surfaces. The calculated potential profiles between cathode and anode surface for the systems with or without OH- and H+ ions are in line with the preset applied potentials. With calibration of offsetting electronegativity on each electrode atom, the electronegativity offset method can simulate electrochemical phenomena.

Acknowledgements We gratefully acknowledge the support of the National Science Foundation CBET-1264173. 1. 2.

3.

4.

5. 6.

7. 8.

Vega, C. and J.L.F. Abascal, Simulating water with rigid non-polarizable models: a general perspective. Physical Chemistry Chemical Physics, 2011. 13(44): p. 19663-19688. Patra, M. and M. Karttunen, Systematic comparison of force fields for microscopic simulations of NaCl in aqueous solutions: Diffusion, free energy of hydration, and structural properties. Journal of Computational Chemistry, 2004. 25(5): p. 678-689. Joung, I.S. and T.E. Cheatham, Molecular Dynamics Simulations of the Dynamic and Energetic Properties of Alkali and Halide Ions Using Water-Model-Specific Ion Parameters. Journal of Physical Chemistry B, 2009. 113(40): p. 13279-13290. Li, P.F., et al., Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent. Journal of Chemical Theory and Computation, 2013. 9(6): p. 2733-2748. Spohr, E. and K. Heinzinger, Molecular Dynamics Simulation of a Water Metal Interface. Chemical Physics Letters, 1986. 123(3): p. 218-221. Yeh, I.C. and M.L. Berkowitz, Structure and dynamics of water at water vertical bar Pt interface as seen by molecular dynamics computer simulation. Journal of Electroanalytical Chemistry, 1998. 450(2): p. 313-325. Panczyk, T., et al., Dynamics of water adsorption on Pt{110}-(1x2): A molecular dynamics study. Journal of Chemical Physics, 2009. 131(6): p. 11. Yeh, I.-C. and M.L. Berkowitz, Aqueous solution near charged Ag(111) surfaces: comparison between a computer simulation and experiment. Chemical Physics Letters, 1999. 301(1–2): p. 81-86.


Page 14 of 17

Page 15 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


9. 10.

11. 12. 13. 14. 15.

16. 17. 18.

19. 20. 21. 22.

23. 24.

25. 26. 27. 28. 29.

Guymon, C.G., et al., Simulating an electrochemical interface using charge dynamics. Condensed Matter Physics, 2005. 8(2): p. 335-356. Yeh, K.Y., M.J. Janik, and J.K. Maranas, Molecular dynamics simulations of an electrified water/Pt(111) interface using point charge dissociative water. Electrochimica Acta, 2013. 101: p. 308-325. Hens, A., G. Biswas, and S. De, Evaporation of water droplets on Pt-surface in presence of external electric field-A molecular dynamics study. Journal of Chemical Physics, 2015. 143(9). Phillpot, S.R. and S.B. Sinnott, Simulating multifunctional structures. Science, 2009. 325(5948): p. 1634-1635. Liang, T., et al., Reactive Potentials for Advanced Atomistic Simulations, in Annual Review of Materials Research, Vol 43, D.R. Clarke, Editor. 2013, Annual Reviews: Palo Alto. p. 109-129. Jiang, H., et al., Hydrogen-Bonding Polarizable Intermolecular Potential Model for Water. The Journal of Physical Chemistry B, 2016. 120(48): p. 12358-12370. Asthana, A. and D.R. Wheeler, A polarizable reactive force field for water to enable molecular dynamics simulations of proton transport. The Journal of Chemical Physics, 2013. 138(17): p. 174502. Hofmann, D.W.M., L. Kuleshova, and B. D’Aguanno, A new reactive potential for the molecular dynamics simulation of liquid water. Chemical Physics Letters, 2007. 448(1–3): p. 138-143. Antony, A.C., et al., Effect of Surface Chemistry on Water Interaction with Cu(111). Langmuir, 2016. 32(32): p. 8061-8070. van Duin, A.C.T., et al., Development and Validation of a ReaxFF Reactive Force Field for Cu Cation/Water Interactions and Copper Metal/Metal Oxide/Metal Hydroxide Condensed Phases. Journal of Physical Chemistry A, 2010. 114(35): p. 9507-9514. Bedrov, D., J. Vatamanu, and Z.Z. Hu, Ionic liquids at charged surfaces: Insight from molecular simulations. Journal of Non-Crystalline Solids, 2015. 407: p. 339-348. van Duin, A.C.T., et al., ReaxFF: A reactive force field for hydrocarbons. Journal of Physical Chemistry A, 2001. 105(41): p. 9396-9409. Shin, Y.K., et al., Variable charge many-body interatomic potentials. Mrs Bulletin, 2012. 37(5): p. 504-512. Liang, T., et al., Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials. Materials Science & Engineering R-Reports, 2013. 74(9): p. 255-279. Spohr, E., Molecular simulation of the electrochemical double layer. Electrochimica Acta, 1999. 44(11): p. 1697-1705. Crozier, P.S., R.L. Rowley, and D. Henderson, Molecular dynamics calculations of the electrochemical properties of electrolyte systems between charged electrodes. The Journal of Chemical Physics, 2000. 113(20): p. 9202-9207. Guymon, C.G., et al., Effects of solvent model flexibility on aqueous electrolyte behavior between electrodes. The Journal of Chemical Physics, 2003. 118(22): p. 10195-10202. Wiebe, J. and E. Spohr, Double layer effects in a model of proton discharge on charged electrodes. Beilstein Journal of Nanotechnology, 2014. 5: p. 973-982. Siepmann, J.I. and M. Sprik, Influence of surface topology and electrostatic potential on water/electrode systems. The Journal of chemical physics, 1995. 102(1): p. 511-524. Reed, S.K., O.J. Lanning, and P.A. Madden, Electrochemical interface between an ionic liquid and a model metallic electrode. The Journal of Chemical Physics, 2007. 126(8): p. 084704. Limmer, D.T., et al., Charge Fluctuations in Nanoscale Capacitors. Physical Review Letters, 2013. 111(10): p. 106102.



30. 31. 32. 33. 34. 35. 36.

37.

38. 39. 40. 41. 42.

43.

44. 45. 46. 47. 48.

Onofrio, N. and A. Strachan, Voltage equilibration for reactive atomistic simulations of electrochemical processes. Journal of Chemical Physics, 2015. 143(5): p. 8. Berendsen, H.J.C., J.R. Grigera, and T.P. Straatsma, The missing term in effective pair potentials. The Journal of Physical Chemistry, 1987. 91(24): p. 6269-6271. Stillinger, F.H. and A. Rahman, Revised central force potentials for water. The Journal of Chemical Physics, 1978. 68(2): p. 666-670. Parr, R.G., et al., Electronegativity: The density functional viewpoint. The Journal of Chemical Physics, 1978. 68(8): p. 3801-3807. Parr, R.G. and R.G. Pearson, Absolute hardness: companion parameter to absolute electronegativity. Journal of the American Chemical Society, 1983. 105(26): p. 7512-7516. Onofrio, N., D. Guzman, and A. Strachan, Atomic origin of ultrafast resistance switching in nanoscale electrometallization cells. Nature Materials, 2015. 14(4): p. 440-446. Antony, A.C., et al. Simulating an Applied Voltage in Molecular Dynamics Using Charge Optimized Many Body (COMB3) Potentials. in ECS Transactions. 2015. The Electrochemical Society. Liang, T., et al., Molecular dynamics simulations of CO2 reduction on Cu(111) and Cu/ZnO(10 (1)over-bar 0) using charge optimized many body potentials. Catalysis Communications, 2014. 52: p. 84-87. Rick, S.W., S.J. Stuart, and B.J. Berne, Dynamical fluctuating charge force-fields - application to liquid water. Journal of Chemical Physics, 1994. 101(7): p. 6141-6156. Mortier, W.J., K. Vangenechten, and J. Gasteiger, Electronegativity equalization - application and parametrization. Journal of the American Chemical Society, 1985. 107(4): p. 829-835. Streitz, F.H. and J.W. Mintmire, Electrostatic Potentials for Metal-Oxide Surfaces and Interfaces. Physical Review B, 1994. 50(16): p. 11996-12003. Wolf, D., et al., Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r(-1) summation. Journal of Chemical Physics, 1999. 110(17): p. 8254-8282. Chen, J. and T.J. Martínez, QTPIE: Charge transfer with polarization current equalization. A fluctuating charge model with correct asymptotics. Chemical Physics Letters, 2007. 438(4–6): p. 315-320. Bockris, J.O.M., M.A.V. Devanathan, and K. Muller, On the Structure of Charged Interfaces. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 1963. 274(1356): p. 55-79. Simon, P. and Y. Gogotsi, Materials for electrochemical capacitors. Nat Mater, 2008. 7(11): p. 845-854. Plimpton, S.J. and A.P. Thompson, Computational aspects of many-body potentials. Mrs Bulletin, 2012. 37(5): p. 513-521. Nose, S., A Unified Formulation of the Constant Temperature Molecular-Dynamics Methods. Journal of Chemical Physics, 1984. 81(1): p. 511-519. Hoover, W.G., Canonical dynamics - equilibrium phase-space distributions. Physical Review A, 1985. 31(3): p. 1695-1697. Adelman, S.A. and J.D. Doll, Generalized Langevin equation approach for atom-solid-surface scattering - general formulation for classical scattering off harmonic solids. Journal Chemical Physics, 1976. 64(6): p. 2375-2388.


Page 16 of 17

Page 17 of 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic for this paper


Applied Potentials in Variable Charge Reactive Force Fields for

Recommend Documents