Self-Diffusion of Surface Defects at CopperâWater Interfaces - The

Article pubs.acs.org/JPCC

Self-Diffusion of Surface Defects at Copper−Water Interfaces Suresh Kondati Natarajan and Jörg Behler*,† Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, D-44780 Bochum, Germany S Supporting Information *

ABSTRACT: Solid−liquid interfaces play an important role in many fields like electrochemistry, corrosion, and heterogeneous catalysis. For understanding the related processes, detailed insights into the elementary steps at the atomic level are mandatory. Here we unravel the properties of prototypical surface-defects like adatoms and vacancies at a number of copper−water interfaces including the low-index Cu(111), Cu(100), and Cu(110), as well as the stepped Cu(211) and Cu(311) surfaces. Using a first-principles quality neural network potential constructed from density functional theory reference data in combination with molecular dynamics and metadynamics simulations, we investigate the defect diffusion mechanisms and the associated free energy barriers. Further, the solvent structure and the mobility of the interfacial water molecules close to the defects are analyzed and compared to the defect-free surfaces. We find that, like at the copper−vacuum interface, hopping mechanisms are preferred compared to exchange mechanisms, while the associated barriers for hopping are reduced in the presence of liquid water. The water structure close to adatoms and vacancies exhibits pronounced local features and differs strongly from the structure at the ideal low-index surfaces. Moreover, in particular at Cu(111) the adatoms are very mobile and hopping events along the surface are more frequent than the exchange of coordinating water molecules in their local environment. Consequently, adatom self-diffusion processes at Cu(111) involve entities of adatoms and their associated solvation shells.

■

INTRODUCTION Defects like adatoms, vacancies, and steps are omnipresent at metal surfaces and play an important role in their chemical reactivity and properties. A thorough characterization of complex surface-imperfections is already challenging for surfaces in vacuum, in particular if large-scale defects or combinations of different defects are involved. In case of solid− liquid interfaces, the complexity is further increased by the mutual interplay between surface features and the solvent. The presence of solvent molecules can in principle stabilize or destabilize certain types of defects and even changes in the overall structure of the metal surface are possible, like in the prominent example of step bunching.1,2 Further, the kinetics of structural changes and chemical reactions at surfaces can be modified drastically. In turn, also the structural and dynamical properties of the interfacial solvent molecules can be very different from the bulk and may depend on subtle geometric details of the surface. Experimental studies on metal−water interfaces are very challenging due to the need for interface sensitive techniques and the complex properties of the interfacial water layers.3 Several experimental techniques like scanning tunneling microscopy (STM) and atomic force microscopy (AFM) have enabled studying a wide range of solvent coverages, from monomers to the bulk, at metal surfaces. Combining STM with density functional theory (DFT), a number of interesting observations on the nature of metal−water interactions have been made so far.4−10 However, most of the available © 2017 American Chemical Society

information on metal−water interfaces pertains to ideal lowindex surfaces without defects. The self-diffusion of defects like adatoms in the presence of bulk water is a little understood phenomenon that is of high practical relevance in various research areas like heterogeneous catalysis at solid−liquid interfaces, electrochemical energy conversion, and surface growth. The most probable diffusion mechanism can be identified from a set of possible mechanisms by comparing their free energy barriers. However, for the simpler metal−vacuum interface, numerous results from experimental techniques like field ion microscopy (FIM),11 STM,12 and He scattering13 and from theoretical calculations using a variety of interaction models, from simple pair potentials with effective medium theory,14−16 embedded atom methods,17,18 tight binding schemes19,20 up to first-principles calculations,21−27 have aided researchers to understand defect diffusion mechanisms in vacuum in detail.28 The theoretical calculation of energy barriers for diffusion processes at surfaces in vacuum usually involves suitably constrained structural relaxations of several intermediate states of the system along the reaction pathway using electronic structure methods and finding the energy difference between the starting structure and the transition state.14,21 Popular algorithms like the nudged elastic band method29 are routinely Received: December 16, 2016 Revised: January 25, 2017 Published: February 2, 2017 4368

DOI: 10.1021/acs.jpcc.6b12657 J. Phys. Chem. C 2017, 121, 4368−4383

Article

The Journal of Physical Chemistry C

symmetry functions characterizing the atom, fulfills a number of important conditions for machine learning potentials, like rotational, translational and permutation invariance of the structural description. For each atom, the symmetry function vector is then fed into an individual atomic neural network, which connects this environment to the atomic energy contribution Ei. For each chemical element in the system, the atomic NNs are constrained to have the same architecture and parameters, which is another requirement for achieving permutation invariance of the NNP total energy expression. The fitting parameters of the atomic NNs are determined in an iterative gradient-based optimization process using a known set of reference data from electronic structure calculations. Once a set of parameters has been found, which accurately reproduces the reference energies and, if available, the forces, the parameters are frozen and the NNP can be used to predict the energies and analytic energy gradients of similar structures, which can contain much larger numbers of atoms. Numerous validation schemes for assessing the accuracy of NNPs are available, and the errors with respect to the chosen reference electronic structure method is typically in the order of a few meV/atom for total energies and of about 0.1 eV/bohr for the forces. The evaluation of the NNP scales close to linear with respect to the number of atoms in the system, and about 100− 200 atoms can be calculated per second and compute core. The technical details for the construction and application of highdimensional NNPs, which have been applied successfully to systems as diverse as isolated molecules,52,53 semiconductors,54,55 metals,56 molecular57,58 and metal clusters,59−61 as well as to liquid water62,63 and electrolyte solutions,64,65 can be found in several recent review articles.50,66

employed to identify the saddle points between the reactant and product configurations.20 Unfortunately, such techniques are unfeasible for computing free energy barriers at metal−water interfaces using costly firstprinciples methods due to the presence of a large number of mobile solvent molecules that need to be thoroughly sampled at the temperature of interest. This requires the use of molecular dynamics (MD) simulations or even accelerated schemes. In particular, ab initio MD30 is nowadays the goldstandard for simulations of solid−liquid interfaces,31−33 but for large systems including a variety of surface imperfections and a substantial number of solvent molecules simulations enabling a comprehensive sampling are computationally too expensive. Consequently, a lot of effort has been spent on the development of more efficient approaches to calculate the potential energy surface (PES) of these systems.34−36 In recent years, a novel type of accurate PESs for complex systems based on machine learning (ML) methods37−41 has emerged, which offers a promising alternative for describing solid−liquid interfaces with close to first-principles accuracy. In this work we use high-dimensional neural network potentials (NNPs),42 a particular class of ML potentials, based on DFT data, to investigate the self-diffusion mechanisms of adatoms and vacancies at a variety of copper surfaces in contact with bulk water. Apart from determining the free energy diffusion barriers at the solid−liquid interface, we analyze in detail the changes in the structural and dynamical properties of the interfacial water molecules induced by surface defects, which are expected to play a crucial role in the efficiency of catalytic reactions and electrochemistry. Consequently, the main goal of the present study is to gain fundamental insights into the properties of prototypical solvated surface imperfections and how they, in turn, affect the solvent molecules at the solid−liquid interface. Due to the efficiency of the employed NNP, we are able to perform long simulations employing very large supercells to ensure that converged properties are obtained.

■

COMPUTATIONAL DETAILS The NNP used in the present work has been generated previously for investigating solid−liquid interfaces of the ternary Cu/O/H system, and a detailed description as well as tests for its validation are given elsewhere.67 In essence, the atomic NNs for copper, oxygen, and hydrogen consist of two hidden layers with 25 artificial neurons each, and contain 53, 53, and 51 atom-centered symmetry functions,51 respectively, for describing the energetically relevant local atomic environments up to a cutoff radius of 12 bohr. The parameters defining these symmetry functions can be found in the Supporting Information of ref 67. The DFT reference calculations required for the determination of the NNP weight parameters have been carried out using the VASP code68 employing the RPBE functional69 for electronic exchange and correlation. The obtained DFT energies and forces have been corrected to include approximate dispersion interactions using the D3 method of Grimme et al.70 in order to obtain a reliable description of liquid water and the solid−liquid interface.57,62,71,72 Projector augmented waves73 (PAWs) with a plane wave energy cutoff of 700 eV are used to describe the electron−core interaction in the DFT calculations. PBE-based PAW potentials with core radii of 2.3, 1.1, and 0.8 Å for Cu, O, and H, respectively, are used as taken from the VASP pseudopotential library.74 The reference data set contains the single-point DFT energies and force components of in total 10,293 structures with up to 384 atoms and consists of geometries of liquid water, ice, bulk copper, bulk cuprous oxide, and copper−water interface geometries as described in ref 67. The NNP predicts the DFT energies and forces of structures

■

HIGH-DIMENSIONAL NEURAL NETWORK POTENTIALS Neural networks have first been used for the construction of interatomic potentials employing electronic structure data by Doren and co-workers in 1995,43 and in recent years numerous alternative approaches to construct ML potentials have been proposed by many groups in this rapidly developing field.44−47 To date NNPs have been reported in the literature for a variety of systems.48,49 In 2007, the method has been generalized by Behler and Parrinello to high-dimensional PESs42 now enabling simulations of condensed systems including thousands of atoms.50 The central idea of this approach, which is also used in the present work, is to express the total potential energy E of a given atomic configuration as a sum of atomic energy contributions Ei, which depend on the energetically relevant local chemical environments. These environments include all neighboring atoms up to a cutoff radius, which needs to be converged for the system of interest to avoid that important interactions are truncated. Typically, cutoffs between 6 and 10 Å have been found to be sufficient for this purpose. The positions of all atoms in the environment of a particular atom i are described by a vector of many-body atom-centered symmetry functions51 Gi. This vector, which can be obtained by a transformation of the Cartesian coordinates to a set of 4369


Article


control of the CV sampling, which usually increases with time in the case of traditional metadynamics simulations. Here, we use the well-tempered approach to metadynamics to study the low barrier diffusion mechanisms like adatom hopping at clean surfaces and the conventional metadynamics to obtain the other barriers.

not included in the training set with root mean squared errors of about 0.9 meV/atom and 66.3 meV/bohr, respectively. The MD simulations have been performed using the LAMMPS molecular dynamics package75 with an additional pair style providing the neural network energies, forces, and stress tensors for the MD algorithm in every time-step. MD simulations of the clean interfaces and of interfaces with single adatoms or vacancies placed at one side of the slab only have been carried out in a 7-layered slab of a (10,0) × (−5,10) supercell for Cu(111) surface following the notation of Park and Madden,76 a seven-layered slab of a (6√2 × 6√2)R45° supercell of the Cu(100) surface, and a 12-layered slab of a (6 × 8) supercell of the Cu(110) surface. Additionally, MD simulations of interfaces with single adatoms have been carried out for stepped surfaces using a 8-layered (2 × 3) surface cell of the Cu(211) and a 12-layered (2 × 3) surface cell of the Cu(311) surface. These stepped surfaces were generated using the atomic simulation environment (ASE).77 Typically, the resulting systems contain approximately 2000 atoms. A time step of 0.5 fs has been used in all simulations. For all surfaces, the vacuum region between the repeated metal slabs has been completely filled with water resulting in a water film with a diameter of approximately 40 Å. Prior to the simulations, all supercells have first been equilibrated with respect to pressure (1 bar) by attaching a barostat along the surface normal and regarding temperature (300 K) in the NPT ensemble for about 1 ns. Consequently, the density obtained in the center of the water film corresponds to the RPBE-D3 equilibrated density in bulk liquid water as demonstrated in previous work,67 where also detailed convergence tests concerning the required system size can be found. This initial equilibration is followed by 1 ns long production simulations sampling the NVT ensemble of the systems at 300 K. Properties like radial distribution functions (RDF), residence lifetimes, and density distributions as presented in the results section are computed from the resulting NVT trajectories. The free energy barriers (FEB) of the defect diffusion processes have been identified employing metadynamics simulations,78 which have been performed using the colvars package79 in LAMMPS. In order to activate a particular defect diffusion mechanism, a set of collective variables (CV) specific to the corresponding diffusion pathway has been chosen and sampled. There is in principle a wide range of possible CVs to choose from, e.g., distances or distance components between atoms, coordination numbers, and interatomic angles. Additionally, constraints can be applied to the CVs in order to enhance the sampling of the interesting regions in the chosen CV space. Details about the CVs chosen in the present work can be found in the Supporting Information. In metadynamics simulations, a history-dependent potential Vmeta acting on the CVs is constructed step-by-step during the course of the simulation by adding repulsive Gaussian potentials at certain time intervals to facilitate a thorough sampling of the CVs. The height and width of the Gaussians depends on the particular mechanism studied, whereas the Gaussian centers are placed at the visited CV coordinate values encountered so far in the simulation. At the end of the simulations, for each investigated mechanism, the accumulated Vmeta is analyzed to identify the transition states along the diffusion pathway and to quantify the FEBs of these processes. A slightly modified scheme called welltempered metadynamics80 has been proposed that enables attaining convergence in metadynamics simulations in the long time limit by allowing the user to modify the temperature

■

RESULTS AND DISCUSSION Adatom Mobility. As a first step of the investigation of adatom diffusion at the low-index surfaces of copper in the presence of bulk water, we have performed a qualitative analysis of the mobility of single adatoms by plotting the density distributions at the Cu(111), Cu(100), and Cu(110) surfaces as obtained in the 1 ns NVT MD simulations at 300 K. The results plotted in Figure 1 show a significant mobility at the Cu(111) surface (Figure 1a) between different hollow sites, which are connected via pathways along the bridge sites. The individual site preferences of the adatom at Cu(111) are given in Figure 2, which shows the relative probability of finding the adatom at the fcc, hcp, bridge, and atop sites as a function of the averaged distance from the nearest copper neighbors in the respective site. As expected, the fcc sites have the highest probability to be populated by the adatom followed closely by the hcp sites, and clearly reduced probabilities for bridge sites and atop sites are observed. The residence lifetimes of the adatoms at these sites are consistent with these probabilities as the lifetimes are substantially longer at the fcc and hcp sites (1.9 and 1.6 ps, respectively) than at the bridge (0.4 ps) and atop (0.1 ps) sites. Contrarily, the adatom at the Cu(100) surface (Figure 1b) is essentially immobile at 300 K and remains in the initial 4-fold hollow site, while the adatom at the Cu(110) surface made only a single hop along the trench between the closely packed atomic rows resulting in a double maximum in the density distribution (Figure 1c). This is a consequence of the high free energy barriers for adatom hopping at these surfaces as discussed below. Properties of the Interfacial Water. The properties of the interfacial solvent molecules and possible changes induced by the adatoms are of high interest. In Figure 3a−c the bulk normalized distributions of the water oxygen atoms along the surface normal at the ideal low-index interfaces are compared with the corresponding distributions in the presence of an adatom for Cu(111), Cu(100), and Cu(110), respectively. Due to the large supercells used in the present work, the oxygen distribution is dominated by the surface regions far from the adatoms, and consequently, the global differences are only small. In essence, the double peak structure of the oxygen distribution in the first hydration layer of the surface, which is discussed in detail in ref 67 for the ideal surface, is maintained with a small reduction in the height of the first subpeak region I due to the exclusion volume of the adatom. Further, the presence of the adatom results in a slight increase in the oxygen density in the second subpeak region II. In addition to these changes in the average water oxygen distribution, we have determined the residence lifetimes of the water molecules in the first hydration layer to shed light on the influence of the adatoms on the mobility of the water molecules in the first hydration layer of the surface. The residence lifetime τ has been computed by fitting the correlation function cresidence(t ) = 4370

⟨a(0)ab(t )⟩ ⟨a(0)⟩

(1) DOI: 10.1021/acs.jpcc.6b12657 J. Phys. Chem. C 2017, 121, 4368−4383

Article


Figure 2. Relative site preferences of the copper adatom at the Cu(111) surface as computed from a 1 ns NVT MD simulation at 300 K. The horizontal axis shows the average distance ⟨D⟩ of the adatom from its neighboring copper atoms. The residence lifetimes of the adatom at the different sites are given in parentheses in the legend. We find that the adatom resides longest at fcc sites and that the lifetime decreases in the following order: fcc > hcp > bridge ≫ atop.

Cu(111) < Cu(100) < Cu(110), which is a consequence of the stronger interaction of the water molecules with more open surface structures. It is worth mentioning that the mobility within the first hydration layer is much higher, i.e., the lifetime in subregion I is shorter, than the exchange of solvent molecules in the first hydration layer with the bulk liquid. In general, the presence of an adatom increases the lifetimes both in subpeak region I as well as in the total first hydration layer, which is consistent with the general phenomenon that an increased surface roughness results in stronger interaction. Again, however, the plotted lifetimes represent averages over the full simulation cells, and consequently, the effect close to the adatoms is expected to be much stronger than suggested by Figure 3d,e. This is further analyzed in Figure 4, which shows the average residence time correlation functions of the oxygen atoms of the water molecules in the first solvation shell around the adatom (a) and surface atoms of the ideal low-index surfaces (b). For this, we have again used the SSP approach where the water molecules within 2.25 Å from the adatom are considered to be within the “well” region, and the “barrier” is considered to be 3.2 Å from the adatom. As expected, the water molecules around the adatoms are clearly less mobile than at the ideal surface. With and without adatom, the water molecules have the longest residence time at the Cu(110) surface underlining again the relatively strong interaction of the molecules with this open surface. The Cu(111) and Cu(100) interfaces are very similar but exhibit a reversed order for the adatom and the first layer surface atoms. We have already concluded from Figure 1 that the adatom is very mobile at the Cu(111) surface. However, we have extracted the lifetimes of neighboring water molecules from the time correlation functions in Figure 4 by fitting the curves to an exponential function of type e−t/τ, where τ is the lifetime. The resulting lifetimes in Table 1 show that for Cu(111) the water molecules in the first solvation shell of the adatoms reside an order of magnitude longer close to these adatoms than the adatoms stay at a given surface site (see Figure 2). Thus, it can be argued that these water molecules remain attached to the adatom while the adatom is diffusing. To analyze the local structural changes of the water molecules close to the adatoms in more detail, we have calculated the copper−oxygen and copper−hydrogen RDFs of the first layer copper atoms at the ideal low-index surfaces and of the copper adatoms, which are plotted in Figure 5. Due to

Figure 1. Adatom densities (probability distributions) of the copper adatoms at different copper−water interfaces as obtained in 1 ns MD trajectories at 300 K. The black circles show the average positions of the surface copper atoms, and the white boxes represent the simulation cells. The adatom at the Cu(111) surface is very mobile due to a low diffusion barrier (cf. Table 3). The adatom at Cu(100) is basically immobile at the surface, whereas for the adatom at the Cu(110) surface just one hop to a neighboring hollow site has been observed.

to an exponential function of type e−t/τ using stable state pictures81 (SSP) approach. We have used the same definitions of “well” and “barrier” as reported in ref 67. Here, a(0) is 1 if the water molecule is located within the “well” at the correlation time origin and 0 otherwise. Similarly, ab(t) is 1 if the water molecule does not cross the “barrier” until time t and 0 otherwise. The obtained residence lifetimes are compared for the interfaces with and without an adatom in Figure 3d,e. In the absence of the adatom, the residence lifetime is within subpeak I, very close to the surface, as well as the overall lifetimes in the full first hydration layer (regions I + II) increase in the order 4371


Article


Figure 3. Panels a−c show the normalized oxygen distributions of the oxygen atoms along the surface normal in the case of the ideal surfaces (black lines) and surfaces with an adatom (red lines) for Cu(111), Cu(100), and Cu(110), respectively. The data has been normalized to the density in bulk liquid water (blue lines) at the equilibrium density. Due to the low coverage by adatoms in the large simulation cells the oxygen density close to the surface is only slightly decreased. Further the residence lifetimes of the water molecules in the first peak region of the first hydration layer (I) (d) and in the complete first hydration layer (I + II) (e) are plotted for the different interfaces with and without the adatom. Water molecules in these regions reside longer at the Cu(110) surface compared to Cu(100) and Cu(111) surfaces. In general, the presence of adatom slightly increases the residence lifetimes.

Figure 4. Resident time autocorrelation functions of the oxygen atoms in the first solvation shell around the adatom in Figure 5 for the different interfaces (a) and at the ideal low-index surfaces without adatoms (b).

Figure 5. Panels a−c show the radial distribution functions g(r) of the oxygen and hydrogen atoms around the first layer copper atoms of the ideal Cu(111), Cu(100), and Cu(110) surfaces in contact with water. Due to the exclusion volume of the neighboring copper atoms the number of close oxygen and hydrogen atoms is very small for the compact Cu(111) surface and increases in case of more open surfaces like Cu(110). In panels d−f the corresponding radial distribution functions around single copper adatoms at these surfaces are plotted showing that the adatoms are, as expected, significantly more exposed to the solvent molecules. The number of water molecules in the first solvation shell of the adatoms is analyzed in Table 2.

Table 1. Residence Lifetime of Water Molecules within the First Solvation Shell of Adatoms and Surface Atoms of Low Index Copper Surfacesa τ [ps] surface

surface atom

adatom

Cu(111) Cu(100) Cu(110)

2.9 3.5 13.1

27.5 20.0 143.1

positions. Consequently, the adatoms, which are well accessible to the solvent molecules, are coordinated by many water molecules. The trend in the peak heights is opposite to the firstlayer copper atoms at the ideal low-index surfaces because the adatom at the Cu(111) surface is more exposed to the solvent than the adatoms at Cu(100) and even more pronounced at Cu(110). This can also be rationalized by the vertical distance of the adatoms to the first surface layer, which is 1.99 Å for Cu(111), 1.76 Å for Cu(100), and only 1.19 Å for Cu(110). For comparison, we also computed the optimized vertical distances of the adatoms from the first layer surface atoms in vacuum using DFT and found that the values are 1.88, 1.65, and 1.22 Å for the Cu(111), Cu(100), and Cu(110) surfaces. Still, this comparison has to be taken with a pinch of salt as the distances for the solid−liquid interfaces are averaged separations including thermal fluctuations of the MD simulation. Like for the ideal surfaces, in case of the adatoms,

These lifetimes are computed by fitting the lifetime correlation functions in Figure 4 to an exponential function of type exp(−t/τ), where τ is the residence lifetime in ps. a

the exclusion volume of the neighboring copper atoms, the heights of the first peaks for the Cu(111), Cu(100), and Cu(110) surfaces shown in panels a, b, and c, respectively, are rather low. For the same reason, the first peak heights increase in the order Cu(111) < Cu(100) < Cu(110) as a consequence of the more open surface structures. For all surfaces, the oxygen peak is closer to the surface copper atoms than the hydrogen peak. The oxygen and hydrogen RDFs of the copper adatoms at these surfaces shown in Figure 5d−f, respectively, exhibit much higher peaks, in which there are no significant shifts in the peak 4372


Article


Figure 6. Contour plots showing the correlation between the dipole angle β and the OH angle γOH of the water molecules within the first solvation shell around the adatom at the surface as shown in Figure 5d−f for Cu(111), Cu(100), and Cu(110), respectively. The dipole angle β is defined as the angle between the dipole vector and the surface normal, while γOH is the angle enclosed by the OH bonds and the surface normal.

first solvation shell followed by three and two water molecules in 6.3% and 5.4% of the time, respectively, resulting in an average adatom coordination by 2.1 water molecules. At the Cu(100) interface, the probability of finding two water molecules in the first solvation layer is reduced to 67.7%, while in 30.4% of the simulation only one water molecule is present. For the Cu(110) interface the coordination by one water molecule only is clearly dominant (93.7%). The spatial distributions of the water molecules, specifically the densities of the oxygen atoms, in the first hydration layer of the Cu(111)−, Cu(100)−, and Cu(110)−water interfaces are shown in Figures 7−9, respectively, for the defect-free ideal interface, an adatom, and a surface vacancy. The oxygen densities are shown for the full first solvation layer in the panels in the first row, as well as for the regions I and II according to Figure 3 in rows two and three. All contour plots are normalized to the oxygen density in equilibrated bulk liquid water corresponding to a relative density of 1.0. At the defect-free Cu(111) surface shown in Figure 7a, the oxygen atoms in the first hydration layer (I + II) prefer the atop sites, and from the detailed analysis of the contributions of regions I and II, it is evident that this preference is dominated by the first subpeak region I. With an adatom present, the oxygen atom density obtained from the MD trajectory is basically indistinguishable from that of the clean surface since the adatom diffuses rapidly at the Cu(111) surface thus smearing out any additional features. For this reason we have carried out an additional simulation with the adatom frozen in an fcc site resulting in the oxygen atom density shown in Figure 7b. It can be observed that the oxygen atoms close to the adatom form a shell-like distribution. As there is no specific site preference with respect to the underlying copper surface, this distribution seems to be governed mainly by the adatom. Unraveling the contributions of subpeak regions I and II shows that this structure is present mainly in region II, which is rather distant from the surface thus explaining the missing influence of the surface copper atoms. In the plot of region I it can be seen that the adatom creates a substantial exclusion volume without water molecules at the surface, and beyond this water-free region the oxygen distribution is very similar to the adatom-free surface. The vacancy at the Cu(111) surface in Figure 7c did not show any diffusion away from its initial location for the duration of the simulation, and a high oxygen atom density is localized at the vacancy within the first hydration layer. Opposite to the adatom case, the increase in the oxygen density is mainly concentrated in subpeak region I close to the defect,

the oxygen peaks are closer than the hydrogen peaks, but the relative heights are reversed compared to the first layer copper atoms. For the adatom case, the hydrogen peak is lower than the oxygen peak, which is likely to be a consequence of the competition between the angular alignment of the water molecules enforced by the adatom and by the neighboring first layer copper atoms. Figure 6 shows the correlation between the dipole angle β, defined as the angle between the molecular dipole moment and the surface normal, and the angle γOH between the OH bonds and the surface normal of the molecules in the first solvation shell around the adatoms of the studied interfaces. For the Cu(111) surface a strong maximum is found with both angles having values between approximately 80 and 120° indicating that the water molecules around the adatom are essentially aligned parallel to the surface. Further, there is a very small peak with dipole angles around 45°. This peak becomes slightly larger at the Cu(100) interface, and the dipole angle distribution is more spread out. In case of the Cu(110) interface, the dipole angle distribution is very different that most molecular dipoles point away from the surface with angles up to at most 40° with respect to the surface normal. The OH angle distribution also differs from the other two interfaces as the OH angles at Cu(110) show a maximum around 60°, which is also consistent with a more upright orientation of the water molecules around the adatoms. In addition to the RDFs and the angular orientations, we have quantified the distribution of possible coordinations of the adatoms by solvent molecules within the first solvation layer in Table 2. In case of the Cu(111) surface, in around 88.3% of the simulation time only two water molecules are found within the Table 2. Analysis of the Number of Water Molecules N Present in the First Solvation Shell of the Copper Adatoms at the Different Copper−Water Interfaces Shown in Figure 5a Fraction of time steps [%] surface

N=1

N=2

N=3

Navg

Cu(111) Cu(100) Cu(110)

5.4 30.4 93.7

88.3 67.7 6.23

6.3 1.9 0.07

2.1 1.7 1.1

a

Given are the fractions of molecular dynamics time steps for all coordinations N = 1−3 and the average number of water molecules Navg in the first solvation shell for each surface. 4373


Article


Figure 7. Contour plots showing the distribution of the oxygen atoms of the water molecules in regions I, II, and (I + II) of the first hydration layer of the Cu(111)−water interface. The first column shows the ideal surface (a), the second column refers to the surface with an adatom (b), and the third column contains the data for a surface with a vacancy (c). The black circles represent the average position of the surface copper atoms. As the adatom exhibits significant diffusion at this surface, the plots in column (b) have been obtained for an adatom constrained at a fcc site for the entire simulation. The values are normalized to the density of oxygen atoms in liquid bulk water simulation (O atom density = 1.0).

while region II shows a depletion of water molecules. Consequently, the vacancy is dragging the water molecules closer to the vacancy site. The increased water density in region I is mainly associated with top sites in the immediate environment of the vacancy and on top of the vacancy itself, while the second nearest top sites exhibit a slightly decreased oxygen density.

The oxygen density for the Cu(100) surface is shown in Figure 8. For the ideal defect-free surface in panel (a) like for the Cu(111) surface the oxygen atoms are primarily located at the atop sites, and again, this structure is governed by subpeak region I. Still, region II shows a more pronounced preference for hollow and bridge sites. The general trends for the adatom and vacancy cases are qualitatively very similar to the Cu(111) 4374


Article


Figure 8. Contour plots showing the distribution of the oxygen atoms of the water molecules in regions I, II, and (I + II) of the first hydration layer of the Cu(100)−water interface. The first column shows the ideal surface (a), the second column the refers to the surface with an adatom (b), and the third column contains the data for a surface with a vacancy (c). The black circles represent the average position of the surface copper atoms. The vacancy exhibits significant diffusion at this surface, and thus, plots in column (c) have been obtained for a vacancy constrained to a single surface atom site for the entire simulation. The values are normalized to the density of oxygen atoms in liquid bulk water simulation (O atom density = 1.0).

surface. However, for the adatom the shell-like increased oxygen density in region II is less spherical and instead shows density maxima mainly at bridge and hollow sites (Figure 8b). Just beyond these rings of bridge site and hollow site O atoms, the density resembles to that of the clean interface. In contrast to the Cu(111) surface, the vacancy at Cu(100) is rather mobile. For the analysis in Figure 8c we have therefore carried out an additional MD simulation with a frozen vacancy

position. In this way we have obtained similar patterns as observed for the Cu(111) surface. For the defect-free Cu(110) surface shown in Figure 9a, all the atop sites and the long bridge sites exhibit high oxygen densities in the first hydration layer with the region I contributing the high atop site density and region II the high density in the long bridge sites. In the MD simulation with an adatom, the adatom performed one hop, and consequently, the averaged position of the adatom shown as black circle is close 4375


Article


Figure 9. Contour plots showing the distribution of the oxygen atoms of the water molecules in regions I, II, and (I + II) of the first hydration layer of the Cu(110)−water interface. The first column shows the ideal surface (a), the second column the refers to the surface with an adatom (b), and the third column contains the data for a surface with a vacancy (c). The black circles represent the average position of the surface copper atoms. The adatom on this surface made only one hop along the trench to the adjacent 5-fold site in the entire simulation, while the vacancy did not move. The values are normalized to the density of oxygen atoms in liquid bulk water simulation (O atom density = 1.0).

findings regarding regions I and II are also similar to the other two surfaces. Free Energy Barriers. As we have seen in the previous section, there are strong differences in the mobilities of the defects, either adatoms or vacancies, depending on the specific surface. Therefore, to investigate the free energy barriers related to the diffusion mechanisms, we have performed metadynamics simulations to overcome even high energy barriers efficiently. The investigated adatom diffusion mechanisms at the (lowindex) copper−water interfaces are summarized in Table 3,

to a short bridge site and the oxygen plotted in Figure 9b is a superposition of both adatom positions. The high oxygen atom density observed just above the adatom position has its origin mainly in region II and a considerable smearing of continuous high density is also observed between the close packed metal rows adjacent to the adatom. The vacancy at the Cu(110) surface has been found to be immobile in the entire simulation, and high oxygen atom densities in Figure 9c are localized at the position of the vacancy similar to the other surfaces, and the 4376


Article


Table 3. Mechanisms and Free Energy Barriers of Adatom Self-Diffusion Processes at Different Copper−Water Interfacesa

a

The NN values have been calculated for the solid−liquid interface using metadynamics, while the DFT values correspond to the barriers in vacuum as obtained from nudged elastic band calculations. The reported “count” refers to the number of events observed spontaneously in a 1 ns MD trajectory at 300 K. The energy barriers for the Cu(111) hopping process refer to the diffusion from fcc to hcp (and from hcp to fcc in parentheses).

Figure 10. Contour plot showing the free energy barrier associated with the adatom exchange processes at the Cu(100) surface. The vertical displacements Z of the adatom (dark blue, Z1) and the replaced surface atom (light blue, Z2) are used as collective variables in the metadynamics simulations. Additionally, a constraint is introduced for the distance between these two atoms to ensure the correct exchange of the selected atoms. The constraining potential has been tested not to influence the accuracy of the free energy calculation.

4377


Article


Table 4. Compilation of the Energy Barriers Associated with Adatom Self-Diffusion Processes at the Cu(100), Cu(111), and Cu(110) Surfaces in Vacuum as Determined in the Present Work and as Reported in the Literaturea Cu(100) Method DFT EM14 EM16 CEM15 EAM25 EAM82 EAM18 EAM17 EAM83 TB19,84 TB85 TB20 DFT

CPMD22 Experiment

NN

hop present work, NEB (1 e−corr.) (no corr.) (MD/MC-CEM) (CEM)

exchange

482 420 350 425 520 470 500

978 220 1120 1100 430 730

480 530 530 430

800 790 760 700

Cu(111) hop fcc

hop hcp

Adatom in Vacuum 41 31 130

present work

550 670 690 280 360

190

hop short bridge

exchange

1264 480

275 180

997 630

244 250

292

826

419

240 280

1150

300 310

36 28 44

1120

43 520

hop long bridge

53

41

Me-D/Hy-D MD MD NEB GGA (PW91)25 GGA (PW91)86 GGA (PW91)26 LDA26 LDA He scattering13 low E Ion scattering87 STM88

Cu(110) exchange

250 320

1330

300 380

322

1049 1049

351 351

140

950

160

1455

960 890 940 970

1130

37 ± 5 Adatom in Water 23 11

1000

a

EM = effective medium method, CEM = corrected effective medium method, TB = tight binding, Me-D = metadynamics, Hy-D = hyperdynamics, NEB = nudged elastic band, EAM = embedded atom method. All energy values are in meV. For comparison, also the free energy barriers obtained with the NN at the copper−water interface are given. The mechanisms of adatom self-diffusion are shown in Table 3. The “fcc” and “hcp” hops of the adatom at the Cu(111) surface refer to diffusion from fcc to hcp and hcp to fcc, respectively. In some theoretical studies and experiments, the reported barriers for the adatom hop at Cu(111) surface are not explicitly mentioned whether it is from fcc to hcp or from hcp to fcc. In those cases, we have placed the values in the middle of both columns.

which also contains the counts of the spontaneous processes observed in conventional 1 ns MD trajectories and the free energy barriers obtained in the NN-based metadynamics simulations as well as, for comparison, the calculated DFT barriers for the clean surfaces in vacuum. In the schematic plots in the table, three panels are provided for each diffusion mechanism describing, from left to right, the starting state, the transition state, and the final state. For each of the mechanisms, suitable CVs have been chosen along the corresponding adatom diffusing pathway based on the initial and final state (See SI). Then, history-dependent bias functions have been added to the PES at different values of the CVs to drive the adatom from the starting structure to the final structure through the transition state, which does not need to be known beforehand. As an example, we show the procedure involved in extracting the FEB of the adatom exchange mechanism at the Cu(100) interface in Figure 10, where the adatom is colored in dark blue and the atom it exchanges with is colored in light blue. The distances traveled by these atoms along the surface normal direction (Z in this case) have been taken as the two collective variables that have been sampled during the metadynamics simulation. An additional constraint restricting the two atoms from moving too far away from each other is also included to ensure a successful exchange process between these two atoms. The right panel in Figure 10 shows the free energy profile of

the adatom exchange process at the Cu(100)−water interface possessing two minima (blue region) located at (0 Å, 2 Å) and (2 Å, 0 Å) representing the starting state and the final state, respectively. The transition point is located approximately at (1.5 Å, 1.5 Å), and the free energy at this point is found to be 1.13 eV. A similar procedure has been followed for identifying the FEB of all other adatom diffusion mechanisms in Table 3. For fast processes, i.e., low barrier mechanisms, like the adatom hopping mechanism at the Cu(111), Cu(100), and Cu(110) surfaces and the exchange mechanism at the Cu(110) surface, we have employed the well-tempered version of metadynamics in order to obtain smooth free energy profiles. Finally, we have also computed the energy barriers of these diffusion mechanisms in vacuum using DFT nudged elastic band (NEB) calculations29 in VASP with a 400 eV energy cutoff and compared the obtained data to the results available in the literature in Tables 4 and 5. In order to validate that the NN potential provides the correct free energy barriers in the metadynamics simulations of the solid−liquid interfaces, a comparison of the DFT and NN energies of representative structures along the diffusion pathway would be desirable. Unfortunately, due to the large system size we employed in the NN metadynamics, this is impossible even for a few single-point calculations. Therefore, we have split the validation process into two steps. First, we have validated the NN potential for smaller systems, which are 4378


Article


For all adatom hopping mechanisms, the presence of water results in a lowering of the FEBs compared to the barriers in vacuum, although the effect is often quite small. For the exchange mechanism, the picture is not so clear as for Cu(111) and Cu(110) the FEBs are lowered, while for Cu(100) a slightly higher barrier is found. The reason is currently unclear, and it can only be speculated that the spatial solvent structure (cf. Figures 7b and 8b) in region II gives rise to different energetics of displacing the solvating water molecules along with the exchange process. In particular the exchange process in (110) is facilitated since the exchanging atom is not buried in the surface like at the Cu(111) and Cu(100) surfaces. The mechanisms and FEBs of the vacancy hopping processes are compiled in Table 6 and compared to the DFT values in vacuum. Several energy barriers available in the literature can be found in Table 5. Vacancy hopping at the Cu(111) surface has twice the barrier as Cu(100) surface mainly due to the closely packed arrangement of surface atoms at the Cu(111) surface. At Cu(110) surface there are two distinct vacancy hopping mechanisms, one along the close packed row ([1̅10] direction) and the other across the close packed row ([001] direction). The hop along the close packed row is 200 meV more favorable than the hop across it since the transition states in the hop along [1̅10] and hop along [001] directions are at the long bridge site and short bridge site of the layer beneath the surface, respectively. Overall, the barriers in the presence of water are only marginally different from those in vacuum in case of vacancy diffusion. Finally, we have investigated diffusion processes at the stepped Cu(211) and Cu(311) surfaces with and without water. The obtained barriers are compiled in Table 7. It is observed that, as expected, the FEBs of adatom hop across the step edge are approximately twice as large as the hop along the step edge for both surfaces. Interestingly, the presence of water results in lowered barriers for diffusion across step edges, while the barriers for diffusion along step edges are increased.

Table 5. Compilation of the Energy Barriers Associated with Vacancy Self-Diffusion Processes at the Cu(100), Cu(111) and Cu(110) Surfaces in Vacuum as Determined in the Present Work and as Reported in the Literaturea

DFT EM16 EAM18 EAM83 NN

Cu(100)

Cu(111)

hop

hop

Vacancy in Vacuum present work, NEB 276 607 437 618 506 355 582 510 Me-D/Hy-D 440 Vacancy in Water present work 300 600

Cu(110) hop [11̅ 0]

hop [001]

550 921

865

600

800

a

EM = effective medium method, Me-D = metadynamics, Hy-D = hyperdynamics, NEB = nudged elastic band, EAM = embedded atom method. All energy values are in meV. For comparison, also the free energy barriers obtained with the NN at the copper−water interface are given. The mechanisms of vacancy diffusion are illustrated in Table 6.

still accessible by DFT, by recomputing the energies of structures visited in the simulations. For instance, in Figure 11

Figure 11. Comparison of representative DFT and NN relative energies of frames extracted from a metadynamics simulation featuring the adatom exchange mechanism in Cu(100)−water interface shown in Figure 10. For this purpose a 255 atom system comprising a 7layered (2√2 × 2√2)R45° supercell of Cu(100) with eight atoms per layer (additional adatom in hollow site) and 66 water molecules is chosen. These structures have not been used in the training process of the NN, but still there is a good agreement between the DFT and NN predicted energy values as the errors are within 2 meV/atom.

■

CONCLUSION A neural network potential capable of describing copper−water interfaces with DFT accuracy has been used to study the properties of surface defects like adatoms and vacancies at various copper surfaces and the role of these defects on the solvation properties of the surfaces. In 1 ns MD simulations at 300 K rapid diffusion of adatoms at the Cu(111) surface has been found, while adatoms at Cu(100) and Cu(110) are basically immobile. Water molecules close to adatoms are significantly less mobile than at the defect-free ideal low-index surfaces. The distribution of water molecules around the adatoms as observed in radial distribution functions can be rationalized by geometric arguments as the water coordination of the adatom decreases from the closely packed surface Cu(111) to more open Cu(110) via the Cu(100) surface, and an opposite trend is observed for the first layer atoms of the defect-free surfaces. Moreover, the water molecules in the first solvation layer around the adatom have been found to reside substantially longer in that region when compared to those around the first solvation layer of surface atoms. The density distribution of oxygen atoms in the first hydration layer of the interfaces with and without adatoms and vacancies has been analyzed in detail. Clean Cu(111) and Cu(100) surfaces have similar density profiles exhibiting high densities in the atop metal sites, whereas in clean Cu(110) high densities are observed in the atop metal sites as well as the long bridge sites.

the DFT and NN energies for a 7-layered (2√2 × 2√2)R45° supercell of Cu(100) with an adatom and 66 water molecules are compared, which show a close agreement with errors of only about 2 meV/atom, which is the typical error of the NN potential. In a second step, the simulations are then converged with respect to system size as described in ref 67, which is a second possible source of error for the energy barriers. The barriers compiled in Table 3 show that with the exception of the Cu(110) surface, where two hopping mechanisms are possible, one via the long bridge (LB) site and the other via the short bridge (SB) site, the hopping mechanisms generally have substantially lower barriers than the exchange mechanisms, both at the solid−liquid interface and in vacuum. At Cu(110) the very unfavorable hop across the short bridge site has a very high energy barrier, while the hop across the long bridge and the exchange mechanism are energetically very similar. In general, our DFT results for the vacuum interface are in good agreement with the data compiled in Table 4, but it should be noted that most of the compiled values refer to methods typically less accurate than DFT and that also a comparison with available DFT data is complicated by the use of different exchange correlation functionals. 4379


Article

The Journal of Physical Chemistry C Table 6. Free Energy Barriers for Vacancy Diffusion Processes at Different Copper−Water Interfacea

a

The NN values have been calculated for the solid−liquid interface using metadynamics, while the DFT values correspond to the barriers in vacuum as obtained from nudged elastic band calculations. The reported “count” refers to the number of events observed spontaneously in a 1 ns MD trajectory at 300 K.

Table 7. Free Energy Barriers for Self-Diffusion Processes of Adatoms at Stepped Copper Surfacesa

a

The NN values have been calculated for the solid−liquid interface using metadynamics, while the DFT values correspond to the barriers in vacuum as obtained from nudged elastic band calculations. In case of the processes across the step edges the given values correspond to the shown mechanism, while the barriers in parentheses refer to processes in the reverse direction. The reported “count” refers to the number of events observed spontaneously in a 1 ns MD trajectory at 300 K.

The water molecules around the adatoms form pronounced shell structures, while in case of the vacancies in all the studied interfaces we found strongly increased water densities very close to the surface. Finally, free energy barriers of typical defect diffusion processes at low-index and stepped copper surfaces were computed by metadynamics simulations. At the low-index surfaces, the adatom hopping process has a lower free energy diffusion barrier compared to the exchange mechanism. At the Cu(110) surface, however, the exchange mechanism is favored

over hopping through short bridge sites, while the hopping through long bridge site was found to be the most probable hopping mechanism. Vacancy hopping mechanisms on flat copper surfaces have been studied, and hopping at the Cu(100) surface is twice as probable as hopping at Cu(111) and Cu(110). Finally, in stepped surfaces like Cu(311) and Cu(211), hopping along the step edges was found to be more probable than hopping across the step edge. 4380


Article


(10) Mehlhorn, M.; Schnur, S.; Groß, A.; Morgenstern, K. MolecularScale Imaging of Water Near Charged Surfaces. ChemElectroChem 2014, 1, 431. (11) Ayrault, G.; Ehrlich, G. Surface Self-Diffusion on an Fcc Crystal: An Atomic View. J. Chem. Phys. 1974, 60, 281. (12) Van Gastel, R.; Somfai, E.; Van Albada, S. B.; Van Saarloos, W.; Frenken, J. W. M. Nothing Moves a Surface: Vacancy Mediated Surface Diffusion. Phys. Rev. Lett. 2001, 86, 1562. (13) Ernst, H.-J.; Fabre, F.; Lapujoulade, J. Nucleation and Diffusion of Cu Adatoms on Cu(100): A Helium-Atom-Beam Scattering Study. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 1929. (14) Hansen, L.; Stoltze, P.; Jacobsen, K. W.; Nørskov, J. K. SelfDiffusion on Copper Surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 44, 6523. (15) Perkins, L. S.; DePristo, A. E. Self-Diffusion Mechanisms for Adatoms on Fcc (100) Surfaces. Surf. Sci. 1993, 294, 67. (16) Stoltze, P. Simulation of Surface Defects. J. Phys.: Condens. Matter 1994, 6, 9495. (17) Liu, C. L.; Cohen, J. M.; Adams, J. B.; Voter, A. F. EAM Study of Surface Self-Diffusion of Single Adatoms of Fcc Metals Ni, Cu, Al, Ag, Au, Pd, and Pt. Surf. Sci. 1991, 253, 334. (18) Karimi, M.; Tomkowski, T.; Vidali, G.; Biham, O. Diffusion of Cu on Cu Surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 52, 5364. (19) Evangelakis, G. A.; Papanicolaou, N. I. Adatom Self-Diffusion Processes on (001) Copper Surface by Molecular Dynamics. Surf. Sci. 1996, 347, 376. (20) Bulou, H.; Massobrio, C. Mechanisms of Exchange Diffusion on Fcc (111) Transition Metal Surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 205427. (21) Feibelman, P. J. Diffusion Path for an Al Adatom on Al(001). Phys. Rev. Lett. 1990, 65, 729. (22) Lee, C.; Barkema, G. T.; Breeman, M.; Pasquarello, A.; Car, R. Diffusion Mechanism of Cu Adatoms on a Cu(001) Surface. Surf. Sci. 1994, 306, L575. (23) Stumpf, R.; Scheffler, M. Ab Initio Calculations of Energies and Self-Diffusion on Flat and Stepped Surfaces of Al and Their Implications on Crystal Growth. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 4958. (24) Yu, B. D.; Scheffler, M. Anisotropy of Growth of the ClosePacked Surfaces of Silver. Phys. Rev. Lett. 1996, 77, 1095. (25) Boisvert, G.; Lewis, L. J. Self-Diffusion of Adatoms, Dimers, and Vacancies on Cu(100). Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 56, 7643. (26) Fordell, T.; Salo, P.; Alatalo, M. Self-Diffusion on Fcc (100) Metal Surfaces: Comparison of Different Approximations. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 65, 233408. (27) Yildirim, H.; Rahman, T. S. Diffusion Barriers for Ag and Cu Adatoms on the Terraces and Step Edges on Cu(100) and Ag(100): An Ab Initio Study. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 235413. (28) Antczak, G.; Ehrlich, G. Jump Processes in Surface Diffusion. Surf. Sci. Rep. 2007, 62, 39. (29) Henkelman, G.; Uberuaga, B. P.; Jonsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901. (30) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: Cambridge, 2009. (31) Izvekov, S.; Mazzolo, A.; VanOpdorp, K.; Voth, G. A. Ab Initio Molecular Dynamics Simulation of the Cu(110) - Water Interface. J. Chem. Phys. 2001, 114, 3248. (32) Liu, L.-M.; Zhang, C.; Thornton, G.; Michaelides, A. Structure and Dynamics of Liquid Water on Rutile TiO2(110). Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 161415. (33) Groß, A.; Gossenberger, F.; Lin, X.; Naderian, M.; Sakong, S.; Roman, T. Water Structures at Metal Electrodes Studied by Ab Initio Molecular Dynamics Simulations. J. Electrochem. Soc. 2014, 161, E3015.

The present work demonstrates that, employing NNPs, large-scale simulations of solid−liquid interfaces can be carried out with DFT accuracy to investigate the properties of surface imperfections like adatoms and vacancies in the presence of water. While here we have restricted our studies to simple local defects, this work can be considered as a first step toward the general structural and dynamical characterization of large-scale solid−liquid interfaces with possible applications in electrochemistry, corrosion, and heterogeneous catalysis.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b12657. Settings for metadyanamics calculations and the resulting free energy profiles for defect diffusion mechanisms in various copper−water interfaces (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Suresh Kondati Natarajan: 0000-0002-7018-5253 Present Address †

Institut für Physikalische Chemie, Georg-August-Universität Göttingen, Tammannstr. 6, D-37077 Göttingen, Germany. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Cluster of Excellence RESOLV (EXC 1069) funded by the Deutsche Forschungsgemeinschaft and by the DFG projects Be3264/5-1 and Be3264/6-1 (Heisenberg fellowship to J.B.). The authors thank Matti Hellström for fruitful discussions.

■

REFERENCES

(1) Baier, S.; Ibach, H.; Giesen, M. The Instability of Vicinal Electrode Surfaces Against Step Bunching I: Experiment. Surf. Sci. 2004, 573, 17. (2) Ibach, H.; Schmickler, W. The Instability of Vicinal Electrode Surfaces Against Step Bunching II: Theory. Surf. Sci. 2004, 573, 24. (3) Zaera, F. Probing Liquid/Solid Interfaces at the Molecular Level. Chem. Rev. 2012, 112, 2920. (4) Cerda, J.; Michaelides, A.; Bocquet, M.-L.; Feibelman, P. J.; Mitsui, T.; Rose, M.; Fomin, E.; Salmeron, M. Novel Water Overlayer Growth on Pd(111) Characterized with Scanning Tunneling Microscopy and Density Functional Theory. Phys. Rev. Lett. 2004, 93, 116101. (5) Michaelides, A.; Morgenstern, K. Ice Nanoclusters at Hydrophobic Metal Surfaces. Nat. Mater. 2007, 6, 597. (6) Lee, J.; Sorescu, D. C.; Jordan, K. D.; Yates, J. T., Jr Hydroxyl Chain Formation on the Cu(110) Surface: Watching Water Dissociation. J. Phys. Chem. C 2008, 112, 17672. (7) Mehlhorn, M.; Carrasco, J.; Michaelides, A.; Morgenstern, K. Local Investigation of Femtosecond Laser Induced Dynamics of Water Nanoclusters on Cu(111). Phys. Rev. Lett. 2009, 103, 026101. (8) Kumagai, T.; Kaizu, M.; Okuyama, H.; Hatta, S.; Aruga, T.; Hamada, I.; Morikawa, Y. Symmetric Hydrogen Bond in a WaterHydroxyl Complex on Cu(110). Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 045402. (9) Forster, M.; Raval, R.; Hodgson, A.; Carrasco, J.; Michaelides, A. c(2 × 2) Water-Hydroxyl Layer on Cu(110): A Wetting Layer Stabilized by Bjerrum Defects. Phys. Rev. Lett. 2011, 106, 046103. 4381


Article

The Journal of Physical Chemistry C (34) Kohlmeyer, A.; Witschel, W.; Spohr, E. Molecular Dynamics Simulations of Water/Metal and Water/Vacuum Interfaces with a Polarizable Water Model. Chem. Phys. 1996, 213, 211. (35) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C. T.; Pandit, S. A. A Reactive Molecular Dynamics Simulation of the SilicaWater Interface. J. Chem. Phys. 2010, 132, 174704. (36) große Holthaus, S.; Köppen, S.; Frauenheim, T.; Ciacchi, L. C. Atomistic Simulations of the ZnO(1−210)/Water Interface: A Comparison between First-Principles, Tight-Binding, and Empirical Methods. J. Chem. Theory Comput. 2012, 8, 4517. (37) Handley, C. M.; Behler, J. Next Generation Interatomic Potentials for Condensed Systems. Eur. Phys. J. B 2014, 87, 1. (38) Li, Z.; Kermode, J. R.; De Vita, A. Molecular Dynamics With On-The-Fly Machine Learning of Quantum-Mechanical Forces. Phys. Rev. Lett. 2015, 114, 096405. (39) Botu, V.; Ramprasad, R. Learning Scheme to Predict Atomic Forces and Accelerate Materials Simulations. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 094306. (40) Behler, J. Perspective: Machine Learning Potentials for Atomistic Simulations. J. Chem. Phys. 2016, 145, 170901. (41) Botu, V.; Batra, R.; Chapman, J.; Ramprasad, R. Machine Learning Force Fields: Construction, Validation, and Outlook. J. Phys. Chem. C 2017, 121, 511. (42) Behler, J.; Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. Phys. Rev. Lett. 2007, 98, 146401. (43) Blank, T. B.; Brown, S. D.; Calhoun, A. W.; Doren, D. J. NeuralNetwork Models of Potential-Energy Surfaces. J. Chem. Phys. 1995, 103, 4129. (44) Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, Without the Electrons. Phys. Rev. Lett. 2010, 104, 136403. (45) Rupp, M.; Tkatchenko, A.; Müller, K.; von Lilienfeld, O. A. Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning. Phys. Rev. Lett. 2012, 108, 058301. (46) Mills, M. J. L.; Popelier, P. L. A. Polarisable Multipolar Electrostatics from the Machine Learning Method Kriging: An Application to Alanine. Theor. Chem. Acc. 2012, 131, 1137. (47) Caccin, M.; Li, Z.; Kermode, J. R.; Vita, A. D. A Framework for Machine-Learning-Augmented Multiscale Atomistic Simulations on Parallel Supercomputers. Int. J. Quantum Chem. 2015, 115, 1129. (48) Handley, C. M.; Popelier, P. L. A. Potential Energy Surfaces Fitted by Artificial Neural Networks. J. Phys. Chem. A 2010, 114, 3371. (49) Behler, J. Neural Network Potential-Energy Surfaces in Chemistry: A Tool for Large-Scale Simulations. Phys. Chem. Chem. Phys. 2011, 13, 17930. (50) Behler, J. Representing Potential Energy Surfaces by HighDimensional Neural Network Potentials. J. Phys.: Condens. Matter 2014, 26, 183001. (51) Behler, J. Atom-Centered Symmetry Functions for Constructing High-Dimensional Neural Network Potentials. J. Chem. Phys. 2011, 134, 074106. (52) Gastegger, M.; Kauffmann, C.; Behler, J.; Marquetand, P. Comparing the Accuracy of High-Dimensional Neural Network Potentials and the Systematic Molecular Fragmentation Method: A Benchmark Study for All-Trans Alkanes. J. Chem. Phys. 2016, 144, 194110. (53) Kolb, B.; Zhao, B.; Li, J.; Jiang, B.; Guo, H. Permutation Invariant Potential Energy Surfaces for Polyatomic Reactions Using Atomistic Neural Networks. J. Chem. Phys. 2016, 144, 224103. (54) Behler, J.; Martoň aḱ , R.; Donadio, D.; Parrinello, M. Metadynamics Simulations of the High-Pressure Phases of Silicon Employing a High-Dimensional Neural Network Potential. Phys. Rev. Lett. 2008, 100, 185501. (55) Artrith, N.; Morawietz, T.; Behler, J. High-Dimensional NeuralNetwork Potentials for Multicomponent Systems: Applications to Zinc Oxide. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 153101.

(56) Artrith, N.; Behler, J. High-Dimensional Neural Network Potentials for Metal Surfaces: A Prototype Study for Copper. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 045439. (57) Morawietz, T.; Behler, J. A Density-Functional Theory-Based Neural Network Potential for Water Clusters Including van der Waals Corrections. J. Phys. Chem. A 2013, 117, 7356. (58) Kondati Natarajan, S.; Morawietz, T.; Behler, J. Representing the Potential-Energy Surface of Protonated Water Clusters by HighDimensional Neural Network Potentials. Phys. Chem. Chem. Phys. 2015, 17, 8356. (59) Artrith, N.; Hiller, B.; Behler, J. Neural Network Potentials for Metals and Oxides- First Applications to Copper Clusters at Zinc Oxide. Phys. Status Solidi B 2013, 250, 1191. (60) Artrith, N.; Kolpak, A. M. Understanding the Composition and Activity of Electrocatalytic Nanoalloys in Aqueous Solvents: A Combination of DFT and Accurate Neural Network Potentials. Nano Lett. 2014, 14, 2670. (61) Elias, J. S.; Artrith, N.; Bugnet, M.; Giordano, L.; Botton, G. A.; Kolpak, A. M.; Shao-Horn, Y. Elucidating the Nature of the Active Phase in Copper/Ceria Catalysts for CO Oxidation. ACS Catal. 2016, 6, 1675. (62) Morawietz, T.; Singraber, A.; Dellago, C.; Behler, J. How van der Waals Interactions Determine the Unique Properties ofWater. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 8368. (63) Cheng, B.; Behler, J.; Ceriotti, M. Nuclear Quantum Effects in Water at the Triple Point: Using Theory as a Link Between Experiments. J. Phys. Chem. Lett. 2016, 7, 2210. (64) Hellström, M.; Behler, J. Concentration-Dependent Proton Transfer Mechanisms in Aqueous NaOH Solutions: From AcceptorDriven to Donor-Driven and Back. J. Phys. Chem. Lett. 2016, 7, 3302. (65) Hellström, M.; Behler, J. Structure of Aqueous NaOH solutions: Insights from Neural-Network-Based Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2017, 19, 82. (66) Behler, J. Constructing High-Dimensional Neural Network Potentials: A Tutorial Review. Int. J. Quantum Chem. 2015, 115, 1032. (67) Natarajan, S. K.; Behler, J. Neural Network Molecular Dynamics Simulations of Solid-Liquid Interfaces: Water at Low-Index Copper Surfaces. Phys. Chem. Chem. Phys. 2016, 18, 28704. (68) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169. (69) Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved Adsorption Energetics Within Density-Functional Theory Using Revised Perdew-Burke-Ernzerhof Functionals. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 7413. (70) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (71) Tonigold, K.; Groß, A. Dispersive Interactions in Water Bilayers at Metallic Surfaces: A Comparison of the PBE and RPBE Functional Including Semiempirical Dispersion Corrections. J. Comput. Chem. 2012, 33, 695. (72) Forster-Tonigold, K.; Groß, A. Dispersion Corrected RPBE Studies of Liquid Water. J. Chem. Phys. 2014, 141, 064501. (73) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953. (74) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758. (75) Plimpton, S. Fast Parallel Algorithms For Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117, 1. (76) Park, R. L.; Madden, H. H., Jr. Annealing Changes on the (100) Surface of Palladium and Their Effect on CO Adsorption. Surf. Sci. 1968, 11, 188. (77) Bahn, S. R.; Jacobsen, K. W. An Object-Oriented Scripting Interface to a Legacy Electronic Structure Code. Comput. Sci. Eng. 2002, 4, 56. 4382


Article

The Journal of Physical Chemistry C (78) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562. (79) Fiorin, G.; Klein, M. L.; Hénin, J. Using Collective Variables to Drive Molecular Dynamics Simulations. Mol. Phys. 2013, 111, 3345. (80) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. (81) Laage, D.; Hynes, J. T. On the Residence Time for Water in a Solute Hydration Shell: Application to Aqueous Halide Solutions. J. Phys. Chem. B 2008, 112, 7697. (82) Marinica, M.-C.; Barreteau, C.; Desjonquères, M.-C.; Spanjaard, D. Influence of Short-Range Adatom-Adatom Interactions on the Surface Diffusion of Cu on Cu(111). Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 70, 075415. (83) Bal, K. M.; Neyts, E. C. Merging Metadynamics into Hyperdynamics: Accelerated Molecular Simulations Reaching Time Scales from Microseconds to Seconds. J. Chem. Theory Comput. 2015, 11, 4545. (84) Evangelakis, G. A.; Papageorgiou, D. G.; Kallinteris, G. C.; Lekka, C. E.; Papanicolaou, N. I. Self-Diffusion Processes of Copper Adatom on Cu(110) Surface by Molecular Dynamics Simulations. Vacuum 1998, 50, 165. (85) Ndongmouo, U. T.; Hontinfinde, F. Diffusion and Growth on Fcc (110) Metal Surfaces: A Computational Study. Surf. Sci. 2004, 571, 89. (86) Sathiyanarayanan, R.; Einstein, T. L. Ab Initio Calculations of Interactions Between Cu Adatoms on Cu(110): Sensitivity of Strong Multi-Site Interactions to Adatom Relaxations. Surf. Sci. 2009, 603, 2387. (87) Breeman, M.; Boerma, D. O. Migration of Cu Adatoms on a Cu(100) Surface, Studied with Low-Energy Ion Scattering (LEIS). Surf. Sci. 1992, 269, 224. (88) Repp, J.; Meyer, G.; Rieder, K.-H.; Hyldgaard, P. Site Determination and Thermally Assisted Tunneling in Homogenous Nucleation. Phys. Rev. Lett. 2003, 91, 206102.

4383


Self-Diffusion of Surface Defects at CopperâWater Interfaces - The

Recommend Documents

Self-Diffusion of Surface Defects at CopperâWater Interfaces - The