Reactive Molecular Dynamics of the Initial Oxidation Stages of Ni

Article pubs.acs.org/JPCA

Reactive Molecular Dynamics of the Initial Oxidation Stages of Ni(111) in Pure Water: Effect of an Applied Electric Field O. Assowe,*,† O. Politano,†,⊥ V. Vignal,† P. Arnoux,‡ B. Diawara,§ O. Verners,∥ and A. C. T. van Duin∥ †

ICB, Université de Bourgogne, UMR 6303 CNRS, 9 Avenue A. Savary, Dijon, France CEA, DEN, DPC, SCCME, Laboratoire d’Etude de la Corrosion Aqueuse, F-91191 Gif-sur-Yvette, France § LPCS, ENSCP-CNRS (UMR 7045), 11 rue Pierre et Marie Curie, 75231 Paris Cedex 05, France ∥ Department of Mechanical and Nuclear Engineering, Pennsylvania State University, University Park, Pennsylvania 16802, United States ‡

ABSTRACT: Corrosion processes occurring in aqueous solutions are critically dependent upon the interaction between the metal electrode and the solvent. In this work, the interaction of a nickel substrate with water molecules has been investigated using reactive force field (ReaxFF) molecular dynamics simulations. This approach was originally developed by van Duin and co-workers to study hydrocarbon chemistry and the catalytic properties of organic compounds. To our knowledge, this method has not previously been used to study the corrosion of nickel. In this work, we studied the interaction of 480 molecules of water (ρ = 0.99 g·cm−3) with Ni(111) surfaces at 300 K. The results showed that a water “bilayer” was adsorbed on the nickel surface. In the absence of an applied electric field, no dissociation of water was observed. However, the nickel atoms at the surface were charged positively, whereas the first water layer was charged negatively, indicating the formation of an electric double layer. To study the corrosion of nickel in pure water, we introduced an external electric field between the metal and the solution. The electric field intensity varied between 10 and 20 MeV/cm. The presence of this electric field led to oxidation of the metal surface. The structural and morphological differences associated with the growth of this oxide film in the presence of the electric field were evaluated. The simulated atomic trajectories were used to analyze the atomic displacement during the reactive process. The growth of the oxide scale on the nickel surface was primarily due to the movement of anions toward the interior of the metal substrate and the migration of nickel toward the free surface. We found that increasing the electric field intensity sped up the corrosion of nickel. The results also showed that the oxide film thickness increased linearly with increasing electric field intensity.

1. INTRODUCTION The interaction of water with solid surfaces is fundamental to research on such processes as corrosion1 and electrochemical reactions.2,3 Water splitting and many heterogeneous catalytic reactions are based on interactions between H2O and metal electrodes.4,5 Recently, there has been considerable interest in the study of metal−electrolyte solutions, especially in terms of metal−water interface structures, as this promises to be a step toward understanding the nature of many electrochemical and surface phenomena. The structures of such interfaces have been investigated by a variety of experimental and theoretical methods.6,7 Some basic industrial catalytic processes, such as reactions in fuel cells and other electrocatalytic devices, involve the adsorption and dissociation of H2O molecules at metal electrode surfaces.8 One important characteristic of such behavior in an electrochemical environment is the variation of the electrode potential, which modifies the surface charge densities and the electric field strength.9,10 Nickel and nickel-based alloys are vitally important to modern industry because of their ability to withstand a wide variety of severe operating conditions, including corrosive © 2012 American Chemical Society

environments, high temperatures, high stresses, and combinations of all of these factors.11 In the nuclear industry, nickelbased alloys are employed in gas turbine engines and a whole other spectrum of relevant applications, such as high-pressure turbine disks, control rod drive mechanisms, valve stems, springs, ducting,12 and so on. Ni-based alloys possess excellent corrosion resistance,13 which is the most sought-after property in the nuclear industry. This excellent corrosion resistance is largely due to the presence of dense ultrathin oxide/hydroxide films that protect the alloy surfaces in aqueous electrolytes or moist air. Such films act as an effective barrier that separates the substrate from the corrosive environment, thereby protecting the substrate from further corrosion.14,15 Therefore, understanding the corrosion behavior of nickel and the growth of these passive films is of vital importance to both scientific research and engineering. Although experimental techniques that can reliably identify the detailed structure and composition Received: July 12, 2012 Revised: October 15, 2012 Published: October 23, 2012 11796

dx.doi.org/10.1021/jp306932a | J. Phys. Chem. A 2012, 116, 11796−11805

The Journal of Physical Chemistry A

Article

of passive films have been widely studied in the literature, it is nonetheless difficult to observe the earliest stages of oxidation as these techniques do not provide all of the necessary information, primarily because local measurements are difficult at the atomic level. Fortunately, modeling may be used as a complementary technique to study the early stages of oxidation. For first-principles investigations of metal anodic oxidation, tuning of the electrochemical potential remains a challenge. The double-reference method developed by Filhol and Neurock2 set the electrochemical potential by adding (or removing) a noninteger number of electrons in the cell. It has been used successfully for the dissociation of water on a Ni(111) surface.3 However, the use of this method to simulate the whole dynamical process of the growth of passive films on metals is not feasible because of computer resource limitations. Recent developments in reactive force fields have made it possible to study the reactivity of materials via a molecular dynamics (MD) approach. In the field of corrosion, the Streitz−Mintmire potential16 can be used to study the thermal oxidation of metallic surfaces or metal−impurity interactions.17−20 In the early 2000s, the study of hydrocarbon chemistry and catalytic properties of organic compounds motivated the development of a new reactive force field21 (ReaxFF). Since this pioneering work, many articles have been published in the literature concerning ceramics,22 metals and metal oxides,23,24 and metal−hydrocarbon interactions.25 Recent progress has also made it possible to investigate the reactivity of metals with water. This method has been used to study the corrosion of such metals as iron,26 copper,27 and aluminum28 in aqueous solution or in the presence of halide ions. To our knowledge, however, it has not been employed to study the corrosion of nickel exposed to water, providing the motivation for the present investigation. For this study, MD simulations using a ReaxFF (reactive force field) were carried out with the aim of investigating the interaction of water molecules with nickel at 300 K. The mechanistic details of atomistic-scale oxide growth under an applied electric field were studied in order to explain the aqueous oxidation of a Ni(111) substrate. We evaluated the structural and morphological features associated with the growth of the oxide film. The charge distributions and selflimiting oxide thicknesses under various electric field intensities have been used to explain the oxidation process as well as the experimentally observed enhancement of the oxidation kinetics and the density of the grown oxide film.

Etotal = E bond + EvdW + ECoulomb + Eval + Etors + Eunder + Eover + E lp + E H‐bond

(1)

These partial contributions correspond to the bond energies (Ebond), undercoordination penalty energies (Eunder), overcoordination penalty energies (Eover), valence angle energies (Eval), torsion angle energies (Etors), lone-pair energies (Elp), hydrogen-bond energies (EH‑bond), and terms to handle nonbonded interactions, namely, van der Waals (EvdW) and Coulomb (ECoulomb) interactions. All but the last two terms are bond-order dependent, that is, their relative contributions depend on the local environment of each atom. The fundamental assumption of the ReaxFF approach is that the bond order between a pair of atoms, BO′ij, can be obtained directly from the interatomic distance rij according to the relation shown in eq 2: BO′ij = BOijσ + BOijπ + BOijππ ⎡ ⎛ rij ⎞ pbo2 ⎤ ⎡ ⎛ rij ⎞ pbo4 ⎤ ⎢ ⎥ ⎢ = exp pbo1 ⎜ σ ⎟ + exp pbo3 ⎜ π ⎟ ⎥ ⎢⎣ ⎝ r0 ⎠ ⎥⎦ ⎢⎣ ⎝ r0 ⎠ ⎥⎦ ⎡ ⎛ rij ⎞ pbo6 ⎤ + exp⎢pbo5 ⎜ ππ ⎟ ⎥ ⎢⎣ ⎝ r0 ⎠ ⎥⎦

(2)

The three pairs of parameters in eq 2, namely, pbo1/pbo2, pbo3/ pbo4, and pbo5/pbo6, correspond to the bond orders of the σ bond, the first π bond, and the second π bond (ππ), respectively. Calculating the bond orders in this way allows us to distinguish between various contributions from the σ, π, and double-π bonds. In this study, the bond orders were given by a general relationship between bond order and interatomic distance and were updated continuously and dynamically. The bond-order expression included a set of parameters that was adjusted to fit relations such as the ones between volume and energy for crystals with coordination numbers 4, 6, and 8. All other nonbonded interactions (short-range Pauli repulsion and long-range dispersion) were included in the van der Waals term (EvdW). The Coulomb energy (ECoulomb) of the system was taken into account for all atom pairs. A shielded Coulomb potential was used to adjust for orbital overlap between atoms that were close together. Atomic charges were calculated using a geometry-dependent charge distribution that was determined using the electronegativity equalization method (EEM) and the QEq method developed by Mortier et al.,29 Rappé and Goddard,30 and Janssens et al.,31 in which individual atomic charges and the corresponding electrostatic energy E(q) change dynamically over time:

2. COMPUTATIONAL METHODS The Ni(111)−H2O interface was modeled using a ReaxFF reactive force field. This potential has been specifically designed to describe chemical reactivity, the dissociation and formation of chemical bonds, surfaces, defects, diffusion, etc. It has been fully fitted to ab initio quantum-mechanical (QM) calculations. ReaxFF is a special type of force field that bridges the gap between the computational quantum-chemistry method and traditional MD. Our simulations were performed using an MD approach based on the ReaxFF force field developed by van Duin et al.21 The total interaction energy is expressed as a function of the bond order, which depends on the interatomic distance. Bonded interatomic interactions are implemented in terms of bond stretching, angle bending, and torsional rotation within a molecule. The ReaxFF energy terms are shown in eq 1:

E(q) =

N

⎡

i=1

⎣

∑ ⎢⎢χi qi + ηiqi 2 + Tap(rij)kc

⎤ ⎥ (rij3 + γij−3)1/3 ⎥⎦ qiqj

(3)

where qi is the charge on ion i (∑Ni=1qi = 0), N is the total number of ions, q denotes the vector of length N containing the charges, kc =14.4 eV·Å·e−2 is the dielectric constant (e represents the electron charge), χi and ηi are the electronegativity and atomic hardness, respectively, of ion i, Tap(r) is a seventh-order taper function, and γij is a parameter for shielding between atoms i and j. The charge values were determined at each time step of the simulation and were dependent on the geometry of the system. This feature made it possible to describe charge transfer in chemical reactions using 11797



Article

ReaxFF. Details of the implementation of this method can be found in ref 32. The terms that do not depend on bond order (ECoulomb and EvdW) were screened by a taper function and shielded to avoid excessive repulsion at short distances. A detailed description of the individual terms can be found in ref 32.

3. RESULTS AND DISCUSSION 3.1. Force Field Development for Ni/O/H. The ReaxFF force field to describe the O/H interaction with nickel and Ni− O oxides is presented in this section. The development of ReaxFF for the Ni−Ni interaction was extensively reported in a previously published article.33 In that work, to ensure that ReaxFF appropriately treats nickel atoms in a range of chemical environments and configurations with different numbers of near neighbors, the authors trained it to reproduce the energies for expansions and compressions of a variety of nickel crystal structures obtained from QM calculations. Also, the bulk modulus, cohesive energy,34 and Ni(111) and Ni(100) surface energies have been fitted with the present potential.33 For the O−H interaction, we used the parameters determined in a previous approach focused on the simulation of the water− silica interface.35 In that work, water properties such as bond angle and distance, charges, pair correlation data, and selfdiffusion constants were compared with experimental and ab initio data. To validate our ReaxFF description of Ni/O/H, direct comparison with QM results was used. All of the periodic QM calculations were performed with the SeqQuest periodic density functional theory (DFT) code and utilized the spinpolarized Perdew−Burke−Ernzerhof (PBE) flavor of DFT and pseudopotentials (see ref 33 for more details). 3.1.1. Nickel−Oxygen Interaction. We first tested the ability of the ReaxFF potential to predict the stabilities of a variety of crystal structures of the oxides of nickel, such as NiO and Ni2O3. For each system, the energy was calculated by ReaxFF and compared with QM data. This energy was obtained for a broad range of compression and expansion of various oxide phases. We see that ReaxFF gave some discrepancies for the nickel oxide phases, most notably the NiO crystal, where ReaxFF predicted a 5% larger equilibrium volume than QM calculations (Figure 1). ReaxFF correctly described the equation of state of the Ni2O3 phase and was in good agreement with the QM calculations (Figure 2). 3.1.2. Ni(OH)2 Bond Dissociation. The bond energy in the reactive potential was optimized using QM-derived values of bond dissociation profiles of Ni(OH)2 clusters. Figure 3 shows

Figure 2. ReaxFF and DFT equation of state plots for the Ni2O3 crystal.

Figure 3. ReaxFF and DFT data for the OH−NiOH bond dissociation curve.

the OH−NiOH bond dissociation curves calculated using DFT and ReaxFF. The dissociation curves were constructed from the equilibrium geometries through single-point calculations by changing the bond length. ReaxFF gave an equilibrium bond length of 1.72 Å, which is in good agreement with the DFT result (1.77 Å). 3.1.3. O2, O, and OH Adsorption on Ni(111). To provide additional support for the validation of ReaxFF, we studied the adsorption of O2, O, and OH on the Ni(111) surface using ReaxFF and compared the results with data from QM calculations. We evaluated the binding energies for various sites of Ni(111) (fcc, hcp, and top). For O2 on Ni, the most favorable binding site was found to be the top site for both DFT and ReaxFF, which gave values of −30.3 and −35.1 kcal/ mol, respectively (Figure 4). For OH on Ni, QM calculations found that OH binding on fcc sites and OH 2-fold sites (Figure 5f,g) have very close energies. ReaxFF showed that the hcp site is the most stable relative to other sites. For O adsorbed on Ni, ReaxFF predicted hcp sites to be more favorable than fcc sites. This was not the case for the QM calculation. 3.1.4. Heats of Formation. To obtain an accurate description of nickel, the ReaxFF parameters were optimized to fit the heats of formation for nickel and the nickel oxides NiO and Ni2O3 (Table 1). For Ni(fcc), ReaxFF gave a value of −103.71 kcal/mol for the heat of formation, whereas DFT gave a value of −102.0 kcal/mol. This result is in good agreement with the experimental result (−102.4 kcal/mol).36 We also calculated the heats of formation of the condensed phase of oxides. DFT gave heats of formation of −24.2 and −82.15 kcal/ mol for NiO and Ni2O3, respectively, while ReaxFF gave values

Figure 1. ReaxFF and DFT equation of state plots for NiO. 11798



Article

Table 1. Heats of Formation (Oxides)/Atomization Energies [kcal/mol] crystal structure

ReaxFF

DFT

NiO Ni2O3 Ni(fcc)

−37.7 −70.4 −103.7

−24.2 −82.15 −102.0

charges of the gas-phase and bulk water. For gas-phase water, the simulation box dimensions varied from 4 Å × 4 Å × 4 Å to 20 Å × 20 Å × 20 Å, and the water was positioned at the center of the box. The bulk water was created by filling a 15.5 Å × 15.5 Å × 15.6 Å box with 125 water molecules. This resulted in a density of 0.99 g·cm−3. The water molecules were assigned random alignments and locations. We optimized the two ReaxFF and VASP systems by minimizing the total energy. These results show that when the length in the box is small, the water molecule reacts with images due to the periodicity of the system. The molecule is polarized when the cell length is shortened. With increasing box size, we observed a convergence of the charges (Figure 6). The water molecule can be considered noninteracting. The charges obtained with ReaxFF and REPEAT for gas-phase water (Figure 6) and bulk water (Figure 7) were in good agreement with experiment. The Bader method overestimated the atomic charges, as mentioned in the literature.44 3.2. Simulation of the Interaction of Water Molecules with the Nickel Surface. 3.2.1. Surface Geometry Simulations. MD simulations were carried out at room temperature (300 K) to investigate the interaction of water

Figure 4. ReaxFF and DFT binding energies of O2, OH, and O on Ni(111) surfaces. The labels refer to the structures displayed in Figure 5.

of −37.8 and −70.15 kcal/mol, respectively. ReaxFF provided a more negative heat of formation for the oxides than DFT. However, the experimental heat of formation (−57.1 ± 0.3 kcal/mol)37 is even lower than the ReaxFF value, indicating that DFT underestimated the stabilities of these oxides. 3.1.5. H2O Charge Calculation. The charge distribution was computed for two H2O systems: bulk water and an isolated gasphase molecule. The results obtained with ReaxFF, Bader,38 and REPEAT39 were compared with experimental results.40,41 With the Bader and REPEAT methods, the charge distribution of the atoms was obtained by performing postprocessing of Vienna Ab Initio Simulation Package (VASP)42,43 output files. These calculations correspond to the oxygen and hydrogen

Figure 5. Plots of O2, OH, and O binding structures: (a) O2 top-top-bound on Ni; (b) O2 on hcp Ni; (c) O2 top-bound on fcc Ni; (d) O2 on fcc Ni; (e) OH on hcp Ni; (f) OH on fcc Ni; (g) OH 2-fold-bound on Ni; (h) O on hcp Ni; and (i) O on fcc Ni. These labels are used in Figure 4. The DFT and ReaxFF simulations were performed on a seven-layer Ni-slab; this figure displays only the top three Ni layers. 11799



Article

surfaces were oriented in the z direction with areas of ∼675 Å2. Periodic boundary conditions were imposed along x and y on the lateral surface, and fixed boundary conditions were applied along z. A reflecting wall was used at the end of the box to confine the water molecules and avoid the reaction of water with the surface below the sample. Fifteen atomic layers were used in the z direction, giving a total of 1800 Ni(111) atoms (Figure 8). The supercell dimension along the z axis was made large enough to mimic a vacuum-exposed surface.45

Figure 8. Schematic diagram showing the initial configuration of the nickel−water system.

All of the simulations were performed using the LAMMPS code.46 The equations of motion were integrated with a time step of 0.1 fs for both short-range and long-range forces. The initial velocities of the systems were chosen from a Maxwell− Boltzmann distribution at 300 K. Canonical (NVT) MD simulations employing the Nosé−Hoover47,48 thermostat technique were used to study room-temperature nickel oxidation in all cases. Additionally, reflecting boundary conditions were imposed on all atoms and molecules that might have reached the upper limit of the simulation box. The atomic charges were updated at every iteration to minimize the electrostatic energy, subject to the constraints imposed by the principle of electroneutrality. Prior to analysis, the system’s energy was minimized using a conjugate gradient method to eliminate thermal fluctuations. The simulation of the interaction of water molecules with Ni(111) at 300 K was carried out using ReaxFF in the absence of an external electric field. No dissociation of water was observed; the average H−O−H bond angle was 107.6°, and the O−H lengths were 0.99 Å. The water molecules were organized in a “bilayer” structure. This structure was observed experimentally by Nakamura and Ito49 via IR reflection absorption spectroscopy on some metal M(111) surfaces (M = Cu, Ni, Pt) and calculated to appear on several flat transitionmetal surfaces.50,51 Figure 9 presents the distribution of water molecules on a clean Ni(111) surface, where 0.0 Å corresponds to the position of the surface nickel atoms. The first two peaks correspond to the water bilayer adsorbed on the nickel surface. The first layer was adsorbed at a distance of ∼1.0 Å and the second at 2.0 Å from the surface. The third peak at 4.8 Å corresponds to a third water layer. There was a gap of 2.8 Å between the second and third layers. This third layer is attributed to a transition state between the water bilayer adsorbed on the surface and the bulk liquid water. These results are in good agreement with the DFT calculations by Levesque et al.52 We also investigated the charge distribution on the interfaces of the Ni(111) system, and the results are shown in Figure 10. From this, we can observe that the charge for bulk nickel was approximately (0 ± 0.02)e. The charges of oxygen and hydrogen in water fluctuated around −0.7e and +0.3e, respectively. Charge transfer between nickel and water molecules was very low, as the charge of the water adsorbed on the surface was in agreement with the charge of bulk water calculated with ReaxFF.35 The nickel atoms on the surface had

Figure 6. Polarization charges of (a) oxygen and (b) hydrogen for gasphase water as functions of simulation box size in ReaxFF, REPEAT, and Bader calculations, and comparison with experimental measurements.

Figure 7. Polarization charges of oxygen and hydrogen in bulk water obtained with ReaxFF, REPEAT, and Bader.

molecules with the Ni(111) surface. The size of the simulation box was 25.3 Å × 26.3 Å × 54.1 Å. The nickel surfaces were obtained by adding two vacuum slabs, one on each side of the fcc nickel surface. The dimension of the vacuum slab along the z axis was chosen in such a way as to give a water density of 0.99 g·cm−3 with 480 water molecules. Water molecules were randomly inserted in the upper vacuum slab. The nickel 11800



Article

electric field in our system, we used the formalism proposed by Chen and Martinez.53 The electrostatic energy was modified to take into account electronegativity differences and charge conservation. The electric field Eel contributes an additional energy term E(q, Eel) associated with the charges on the ions: E(q, Eel) = E(q) − q T·∑ R νEelν ν

(4)

We use the boldface convention for vectors and matrices acting in the space of atomic charge variables (i.e., q and Rν are vectors of length N). E(q) is the electrostatic energy without the electric field (given by eq 3); qT is the transpose of the vector of charges; Rν is the vector of Rν,i values, where Rν,i is the νth spatial component of the position vector Ri of the ith ion; and Eνel is the νth spatial component of the external electric field Eel. The external electric field simply perturbs each atomic electronegativity by an amount Rν,iEνel, which is the potential produced by the external field. Therefore, we can replace χi by χi = χi − ∑νRν,iEνel to obtain the new charge distribution. Details of this method can be found in the paper by Chen and Martinez.53 Finally, we included additional electrostatic forces Fi = qiEel associated with the electric field for all of the systems. We began by testing these modifications on isolated systems such as a single water molecule in the gas phase (Figure 11)

Figure 9. Water distribution on a clean Ni(111) surface at the end of the simulation.

Figure 10. Charge distribution of the system at the end of the simulation (e denotes the electron charge).

an average charge of +0.101e because they were affected by the presence of the water molecules. However, the average charge of the first water layer was −0.110e, and those of the second and third layers were −0.04e and −0.06e, respectively (Table 2). The atomic charges fluctuated around zero on both sides of

Figure 11. Variation of the atomic charges of oxygen and hydrogen in a single water molecule as functions of the electric field intensity.

Table 2. Average Charge (q) of Subsurface and Surface Nickel and First Three Adsorbed Water Layers qNi (e)

and a Ni(111) slab (Figure 12). For the single water molecule, we determined the charge distributions of oxygen and hydrogen as functions of the electric field intensity. Optimization was achieved with ReaxFF and VASP by minimizing the total energy. The oxygen and hydrogen charges changed with changes in the electric field intensity. We observed polarization of the water molecule due to the presence of the electric field. For the free surface of nickel, we calculated the average charge of the surface nickel atoms using ReaxFF MD at 300 K over a period of 20 ps. This revealed that the surface charge increased with the electric field intensity. The surface was charged negatively when the applied electric field was negative and positively when the field was positive. We concluded that the presence of the electric field is the cause of the polarization of the water molecule and the free surface of nickel. We then studied the influence of the electric field on the oxidation of nickel. In all cases, the electric field was applied 3 Å from the

qH2O (e)

subsurface

surface

first layer

second layer

third layer

−0.016

+0.101

−0.110

−0.04

−0.06

the interface between the surface nickel atoms and the first water layer. The metal surface was positively charged whereas the first water layer was negatively charged, giving rise to an electric double layer. This goes to show that electric double layers may be reproduced in ReaxFF. To increase the reactivity of nickel, we applied an external electric field between the metal and the solution. 3.2.2. Effect of an Applied Electric Field. In this part of our study, we investigated the influence of an electric field on oxide growth. The field intensity at the metal−aqueous solution interface varied between 10 and 20 MeV/cm. To introduce the 11801



Article

The OH− was adsorbed at a distance of 1.0 Å from the nickel surface, forming a low hydroxide layer (Figure 14, top left), and

Figure 12. Average charge across the first layer of nickel atoms in a Ni(111) slab at different electric field intensities.

surface. This field affected the surface nickel atoms and the first water layer. Figure 13 shows a side view of the interface between the nickel substrate and the aqueous solution at the end of

Figure 14. Side views of the dissociation of water on the nickel surface and the dissolution of nickel at stages I (0.7 ps) and II (2.4 ps). Red, white, and green spheres represent O, H, and Ni atoms, respectively.

H+ bonded with another H2O to form H3O+ (Figure 14, top right). This subsequently went into solution. The average charge of the surface nickel atoms was +0.29e, while those of OH− and H3O+ were −0.248e and +0.18e, respectively. Stage II. The second stage was characterized by the dissolution of nickel atoms (Figure 14 bottom). Two nickel atoms were observed to leave the surface, marking the formation of an oxide nucleus. These two atoms were located 2.1 Å from the initial surface, and their charges were +0.48e and +0.44e. The average charge of the surface nickel atoms increased relative to the first stage, reaching +0.31e. We also observed a reinforcement of the hydroxide layer on the surface and the H3O+ molecules in the solution (54 OH− and 54 H3O+ at this stage). Stage III. The third stage began with the penetration of oxygen into the nickel substrate at 3.0 ps (30 000 iterations). This oxygen came from the OH− adsorbed on the surface. The hydrogen reacted with H2O to form H3O+. The oxygen did not penetrate the vacuum left by the nickel atom but instead occupied the space between the nickel atoms at the surface. This oxidation process repeated itself throughout the simulation, thereby fostering oxide growth. This mechanism occurred at all electric field intensities. We found that increasing the field intensity sped up the process and the kinetics of the reaction. The same mechanism was previously observed experimentally by Okamato54 and Sato,55 who demonstrated how it could explain the formation of passive films on metallic surfaces. Figure 15 presents the evolution of the maximum position of nickel and the minimum position of oxygen in the oxide throughout the simulation for an electric field intensity of 17 MeV/cm. In this graph, 0.0 Å corresponds to the initial position of the surface nickel atoms. Before 45 ps, the growth of the oxide scale on the nickel surface was primarily related to the movement of anions toward the interior of the metal substrate and the migration of nickel toward the free surface. After 45 ps, a constant spacing between Ni and O was observed. This last clearly indicated the absence of oxide growth. To our point of view, this feature was correlated with a filling of the oxide layer by oxygen. This point was also observed when we looked at the

Figure 13. Side view of the interface between the nickel substrate and water at 350 ps relaxation for various electric field intensities (0 MeV/ cm to 20 MeV/cm). The red, white and green spheres represent O, H and Ni atoms, respectively.

simulations (350 ps relaxation) with increasing electric field intensities. We observe that the presence of an electric field created a chemical reaction between the nickel atoms and water molecules. We analyzed the nickel oxidation process during the simulation with an electric field intensity of 17 MeV/cm. Several stages were observed: Stage I. After hydration of nickel, the first stage that we observed was the dissociation of water molecules:

H 2O → OH− + H+

(5) 11802



Article

Figure 17. Evolution of the simulated oxide film thickness at room temperature for various electric field intensities.

Figure 15. Evolution of the maximum position of Ni and the minimum position of O in the system for an electric field intensity of 17 MeV/cm.

mode of the oxide layer remained the same. Our simulations indicated that the oxide film growth experienced an initial rapid increase followed by a slower growth phase after 50 ps of simulation time. We observed that the oxidation kinetics and the oxide growth were considerably enhanced in the presence of external electric fields in the 10−20 MeV/cm range at room temperature. Increasing the electric field from 10 to 20 MeV/ cm resulted in faster kinetics and an earlier commencement of the slower growth phase. We also noticed a linear dependence of the film thickness on the applied electric field (Figure 18). This result is in agreement with the linear dependence of film thickness on formation potential that was previously observed by Yoneyama et al.56

evolution of the simulated oxide film thickness (see Figure 17). We also analyzed the evolution of the chemical compositions at an electric field intensity of 17 MeV/cm (Figure 16). The

Figure 16. Variation of the chemical compositions in the system over time for an electric field intensity of 17 MeV/cm.

proportion of water molecules decreased because they dissociated into OH− and H3O+, which reacted with the surface. The number of OH− ions decreased quickly at ∼30 ps in the simulation. After the hydroxide dissociated, the oxygen atoms reacted with the nickel, thus forming the oxide film, while the hydrogen reacted with H2O to form H3O+. Some of this hydrogen was released into the solid. We also found that the number of oxygen anions increased over time. The influence of the electric field on the evolution of the simulated oxide film thickness in Ni(111) oxidation is shown in Figure 17. The atoms belonging to the oxide layer were identified as those whose z positions lay between the outermost nickel atom and the innermost oxygen atom in the oxide. One may notice that all four curves indicate similar behavior. Although increasing the electric field intensity led to faster initiation of oxidation, it did not change the other aspects of the kinetics. The only stage that was influenced by the field intensity was stage I, because of the dissociative chemisorption of water and the dissolution of nickel. Nonetheless, the growth

Figure 18. Thickness of the oxide film as a function of electric field intensity.

At the end of the simulation, the charge distributions of nickel, oxygen, and hydrogen were analyzed for the four electric field intensities, as shown in Figure 19. We observed that in the case of nickel and oxygen, the four curves indicate similar behavior. The charge of the nickel atoms within the substrate fluctuated around 0.0e at all field intensities. We observed an increase in the charge of the nickel at the metal−oxide interfacial layer; furthermore, this increase became more significant with increasing field intensity. The maximum charge of the nickel in the oxide layer was approximately +0.6e. As for 11803



Article

4. CONCLUSION The oxidation of Ni(111) in pure water has been studied using a reactive force field, both with and without an externally applied electric field. In the absence of an electric field, no reaction between the nickel and water molecules was observed. The water was adsorbed as two layers. An electric double layer was detected at the interface; this layer was the result of the positive charge of the surface nickel atoms and the negative charge of the first water layer. The presence of the electric field conditioned the oxidation of nickel in aqueous solution. Different stages of oxidation were observed: adsorption, dissociation of water, deprotonation, dissolution, and growth of hydroxide/oxide. The oxide growth was enhanced by nickel dissolution and oxygen ingress. The oxide film thickness was linearly dependent on the electric field.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +33 0380399190. E-mail: [email protected]. Notes

The authors declare no competing financial interest. ⊥ Tel: +33 0380396172. E-mail: olivier.politano@u-bourgogne. fr.

■

ACKNOWLEDGMENTS The authors acknowledge the use of the DSI-CCUB computing facilities at the Université de Bourgogne. They are also very grateful to the Conseil Régional de Bourgogne and the CEASaclay for financial support of O.A.

■

REFERENCES

(1) Mundy, C. J.; Kuo, I. F. W. Chem. Rev. 2006, 106, 1282−1304. (2) Filhol, J. S.; Neurock, M. Angew. Chem., Int. Ed. 2006, 45, 402− 406. (3) Taylor, C.; Kelly, R. G.; Neurock, M. J. Electrochem. Soc. 2006, 153, E207−E214. (4) Gokhale, A. A.; Dumesic, J. A.; Mavrikakis, M. J. Am. Chem. Soc. 2008, 130, 1402−1414. (5) Phatak, A. A.; Delgass, W. N.; Ribeiro, F. H.; Schneider, W. F. J. Phys. Chem. C 2009, 113, 7269−7276. (6) Nakamura, M.; Ito, M. Chem. Phys. Lett. 2004, 384, 256−261. (7) Izvekov, S. I.; Mazzolo, A.; Opdorp, K. V.; Voth, G. A. J. Chem. Phys. 2001, 114, 3248−3256. (8) van Duin, A. C. T.; Merinov, B. V.; Jang, S. S.; Goddard, W. A., III. J. Phys. Chem. A 2008, 112, 3133−3140. (9) Bruesch, P.; Christen, T. J. Appl. Phys. 2004, 95, 2846−56. (10) Zhang, P.; Zheng, W. T.; Jiang, Q. J. Phys. Chem. C 2010, 114, 19331−19337. (11) Davis, J. R. ASM Specialty Handbook; ASM International: Materials Park, OH, 2000. (12) Young, B. A.; Gao, X.; Srivatsan, T. S. J. Nucl. Mater. 2009, 394, 63−66. (13) Liu, X. B.; Fu, G. Y.; Liu, S.; Shi, S. H.; He, X. M.; Wang, M. D. Nucl. Eng. Des. 2011, 241, 4924−4928. (14) Jang, H. J.; Park, C. J.; Kwon, H. S. Electrochim. Acta 2005, 50, 3503−3508. (15) Magnussen, O. M.; Scherer, J.; Ocko, B. M.; Behm, R. J. J. Phys. Chem. B 2000, 104, 1222−1226. (16) Streitz, F. H.; Mintmire, J. W. Phys. Rev. B 1994, 50, 11996− 12003. (17) Elsener, A.; Politano, O.; Derlet, P. M.; Van Swygenhoven, H. Modell. Simul. Mater. Sci. Eng 2008, 16, No. 025006. (18) Garruchet, S.; Politano, O.; Arnoux, P.; Vignal, V. Appl. Surf. Sci. 2010, 256, 5968−5972.

Figure 19. Charge distributions of (a) nickel, (b) oxygen, and (c) hydrogen as functions of the position along z at 350 ps for various electric field intensities.

the oxygen, we observed that its charge fluctuated around −0.8e, which might correspond to oxygen belonging to the water cluster or H3O+ in aqueous solution. The charge of the oxygen increased with decreasing distance from the oxide/ hydroxide−solution interface. The maximum charge of the oxygen in the oxide layer was approximately −0.2e. The atomic charge of the hydrogen within the solution was approximately +0.35e. The hydrogen that had been released into the solid during the hydroxide−oxide transition had a charge of approximately +0.1e. The charge of the nickel and oxygen within the oxide layer increased with increasing electric field intensity. 11804



Article

(19) Sankaranarayanan, S. K. R. S.; Ramanathan, S. Phys. Rev. B 2008, 78, No. 085420. (20) Garruchet, S.; Politano, O.; Arnoux, P.; Vignal, V. Solid State Commun. 2010, 150, 439−442. (21) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III. J. Phys. Chem. A 2001, 105, 9396−9409. (22) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A., III. J. Phys. Chem. A 2003, 107, 3803−3811. (23) Han, S. S.; van Duin, A. C. T.; Goddard, W. A., III; Lee, H. M. J. Phys. Chem. A 2005, 109, 4575−4582. (24) Zhang, Q.; Ç aǧin, T.; van Duin, A.; Goddard, W. A., III; Qi, Y.; Hector, L. G., Jr. Phys. Rev. B 2004, 69, No. 045423. (25) Nielson, K. D.; van Duin, A. C. T.; Oxgaard, J.; Deng, W.; Goddard, W. A., III. J. Phys. Chem. A 2005, 109, 493−499. (26) Pan, T. Chem. Phys. Lett. 2011, 511, 315−321. (27) Jeon, B.; Sankaranarayanan, S. K. R. S.; van Duin, A. C. T.; Ramanathan, S. J. Chem. Phys. 2011, 134, No. 234706. (28) Russo, M. F.; Li, R., Jr.; Mench, M.; van Duin, A. C. T. Int. J. Hydrogen Energy 2011, 36, 5828−5835. (29) Mortier, W. J.; Ghosh, S. K.; Shankar, S. J. Am. Chem. Soc. 1986, 108, 4315−4320. (30) Rappé, A. K.; Goddard, W. A., III. J. Phys. Chem. 1991, 95, 3358−3363. (31) Janssens, G. O. A.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A. J. Phys. Chem. 1995, 99, 3251−3258. (32) Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P. Comput. Phys. Commun. 2008, 178, 73−87. (33) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A., III. J. Phys. Chem. C 2010, 114, 4939−4949. (34) Assowe, O.; Politano, O.; Vignal, V.; Arnoux, P.; Diawara, B. Diffus. Defect Data, Pt. A 2012, 323, 139−145. (35) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C. T.; Pandit, S. A. J. Chem. Phys. 2010, 132, No. 174704. (36) Kittel, C. Introduction to Solid State Physics, 8th ed.; Wiley: New York, 2005. (37) Boyle, B. J.; King, E. G.; Conway, K. C. J. Am. Chem. Soc. 1954, 76, 3835−3837. (38) Bader, R. F. W.; MacDougall, P. J.; Lau, C. D. H. J. Am. Chem. Soc. 1984, 106, 1594−1605. (39) Campana, C.; Mussard, B.; Woo, T. K. J. Chem. Theory Comput. 2009, 5, 2866−2878. (40) Martin, F.; Zipse, H. J. Comput. Chem. 2004, 26, 97−105. (41) Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. J. Chem. Phys. 1973, 59, 2254−2259. (42) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47, 558−561. (43) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15−60. (44) Martin, F.; Zipse, H. J. Comput. Chem. 2005, 26, 97−105. (45) Politano, O.; Garruchet, S.; Salazar, J. M. Mater. Sci. Eng., A 2004, 749, 387−389. (46) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19 (also see http:// lammps.sandia.gov/). (47) Nosè, S. Mol. Phys. 1984, 52, 255−268. (48) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (49) Nakamura, M.; Ito, M. Chem. Phys. Lett. 2005, 404, 346−350. (50) Michaelides, A.; Alavi, A.; King, D. A. Phys. Rev. B 2004, 69, No. 113404. (51) Michaelides, A. Appl. Phys. A: Mater. Sci. Process. 2006, 85, 415− 425. (52) Levesque, M.; Roques, J.; Domain, C.; Perron, H.; Veilly, E.; Simoni, E.; Catalette, H. Surf. Sci. 2008, 602, 3331−3337. (53) Chen, J.; Martínez, T. J. J. Chem. Phys. 2009, 131, No. 044114. (54) Okamato, G. Corros. Sci. 1973, 13, 471−489. (55) Sato, N. Corrosion 1989, 45, 354−368. (56) Yoneyama, H.; Ochi, K.; Tamura, H. J. Appl. Electrochem. 1976, 6, 283−291.

11805


Reactive Molecular Dynamics of the Initial Oxidation Stages of Ni

Recommend Documents