Category Archives: Transition State

Possible transition processes of Pt adatom on Pt(100) surface

by Hepeng Ye

Introduction:

Surface diffusion is a common process for solid materials, and in our case, the Pt adatom could diffuse along the surface of Pt(100). Diffusion is governed by the thermal energy from the environment and the property of such behavior is important for understanding surface phase formation, heterogeneous catalysis and other aspects.

In our case, Pt adatom is hypothesized to undergo a hopping between two adjacent four-fold sites or undergoes concerted substitution. To study the process, a model needs to be built and in our case the model is fairly simple since neither diffusion process has linear interpolation issue, and as consequence, one just need an initial and final two-state model for the transition state study. By connecting the initial and final state of the system, one could locate the TS (transition state) of that process and use the energy barrier from the transition search to further study the kinetic of process.

Experimental:

First cleave the Pt(100) lattice, and create the slab and a five-layer slab of 2X2 super-cell with 5Å vacuum space was set as the starting point. As mentioned the goal of the experiment is to identify the more likely transition process from the two, and in order to get reasonable comparison between the two processes, one need to first optimized the system(slab) to make sure that the super-cell is large enough for boundary conditions and enough vacuum space to prevent interaction from the next periodic slab along the Z-axis. Layers of the slab will not be optimized due to time limit, and based upon experience, five layers should be enough.

For vacuum space optimization, build three slabs with 5Å, 10Å and 15Å vacuum space, and calculate energy with CASTEP, functional used is GGA PBE, pseudo-potential is OTFG-ultrasoft, and all other parameters are set to default (272.1eV energy cutoff, 6x6x1 K-points, core radii is 1.27A with electron configuration to be 5s2 5p6 5d9 and 6s1).

Figure 1. Example of the slab used for vacuum space optimization. (5layer, with 10A vacuum)

Use the lowest energy vacuum space to continue slab size optimization.

To do this, one will first identify the most stable height of the Pt adatom on the surface. Sitting on top of the fourfold site, it is obvious that the adatom cannot be either too close or too far from the surface where too close to each other there will be repulsion and too far to each other will make Pt adatom isolated from the surface and could interact with the next periodic slab. Consider the short range dispersion, there should be a local energy minimum for the Pt to sit on top of the surface and bring down the system’s energy. So the height is tuned manually by putting the adatom between 1Å to 3Å above the center of the four-fold site.

When Pt adatom transit from one four-fold site to another, the energy for it being reactant should be the same as it being the product. Using this idea, I will separately calculate the energy for Pt adatom sitting at different fourfold sites, for different size of slabs.

The lowest energy height optimized previously will be used, and start with 3×3 super-cell.

Figure 2. Position of initial and final state of fourfold site.

 

 

As shown in the Figure 2, in 3×3 super-cell, the starting point of Pt, say, at the center, and could diffuse into the nearby fourfold sites. Due to symmetry, the energy barriers for all four direction should be the same, so just consider one of the pathways. In this case, energy of the system is calculated for Pt on the center site and on the edged site.

To get energies of each initial and final state, all calculations are performed with CASTEP, functional used is GGA PBE, pseudo-potential is OTFG-ultrasoft, and all other parameters are set to be the same as above when doing slab-size optimization.

After calculations are done, and comparing the two energies and check the relative difference, and based upon that difference, one will decide whether it necessary to keep using 3×3 super-cell or increase the slab size to 4×4 super-cell or even larger.

Then, in order to find out the transition state, one need to build the initial and final state or the reactant and product state respectively. There are two possible diffusion processes, one is the hopping between two fourfold sites, another is the concerted substitution of Pt adatom with Pt on the top layer. And in the calculations, two assumptions will be made to simplify the problem and reduce computational cost:

  1. In either process, Pt atom will only diffuse with the closest Pt atom or fourfold site.
  2. Due to the fact that system used has already optimized its size, and due to the symmetry, only consider Pt adatom diffuse in one direction, instead of four.

Based upon these assumptions, two sets of reactant and product slabs will be built for corresponding diffusion pathways.

 

Figure 3. Transition process scheme for hopping between fourfold states.

Figure 4. Transition process of concerted substitution of Pt adatom with Pt atom on the surface.

After determining the initial and final state of the process, DMol3 TS-search will be used to find the transition state/energy for each process while comparing transitional energy calculated with different search protocols.

Data analysis and discussion:

For the slab size optimization, calculations were performed for 3×3 super-cell where the Pt adatom sits in the fourfold site in the middle, and the fourfold site on the edge, and the energies are close enough to say they are in the same environment.

For the manually adjusted Pt adatom positions, three-point-adjust method is used to iteratively adjust and narrow down the height range and get to the lowest energy and its corresponding height/distance. One could see from the plot that at 1.8Å above the surface, the energy is the lowest.

Figure 5. Energy diagram of Pt adatom with relative heights to the slab surface.

After that, a 5 layer, 3×3 super slab is applied for the Transition-State study.

In Castep transition state search, there are some different synchronous transit methods[3], and here, some of these methods were used and were compared their differences. For the fourfold to fourfold hopping, the transition energies for different TS-search protocols are shown in the Table 1 below:

Table 1. Transition energy of hopping between fourfold sites calculated from different TS-search protocol.

For concerted substitution, the transition energy for each TS-search protocol is also listed in the table2.

Table 2. Transition energy of concerted substitution calculated from different TS-search protocol.

Different transition searching protocols have different calculation process, but for each process, they are all in agreement with the transition energies are on the same magnitude. While for Pt substitution process, the transition energy barrier is ten times larger than that of hopping process, which makes sense that such process requires breaking Pt metallic bonds with the surrounding Pt atoms and then another Pt atom comes in and forms new bonds, there should be a lot more energy required for such reaction.

Conclusion:

Based upon the energy barrier calculated from the Transition-Search, it is obvious that the concerted substitution has way larger energy barrier than the fourfold sites hopping, which means at same conditions(say same temperature) it will be much less likely for Pt adatom diffuse into the Pt(100) surface and undergoes the concerted substitution process.

There are some experimental setups that could be improved, say the slab I used has 5 layers, which I did not optimize. Also, I stopped at 3×3 supercell size and did not look at 4×4 supercell simply because the energy between two fourfold sites are very close to each other(~0.005 hartree for the total system) for 3×3 supercell, but maybe 4×4 will give even more identical energy, I do not know and maybe someone could try 4×4 disregarding the computational cost.

 

Reference:

  1. J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 77 (1996)

2. J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais,      Atoms, Molecules, Solids, And Surfaces – Applications of the Generalized Gradient  Approximation for Exchange and Correlation, Phys. Rev. B, 46 (1992)

3. CASTEP GUIDE, Transition state search task. Date of access: April. 27. 2019

https://www.tcm.phy.cam.ac.uk/castep/documentation/WebHelp/content/modules/castep/tskcasteptss.htm

Calculating transition state and activation energy of Ag diffusion via adjacent threefold sites on Cu(111) surface

Author: Stephen Holoviak

Introduction:

Diffusion of adatoms on a surface is an important phenomena in many processes; from thin film deposition (1) to catalysis. This study is an attempt to understand the activation energy of an Ag atom diffusing between the adjacent threefold sites on the Cu(111) surface.

Figure 1: Threefold sites on Cu(111) surface.

Method:

All calculations were carried out in the Dmol3(2) package in materials studio.

Exchange-correlation functional used was LDA-PWC (3).

Kpoint density: 1x1x1, all meshes were distributed via the Monkhorst-pack method (4). If more computational resources were available, more k-points would be added.

SCF Tolerance: 10E-4 [eV]

Geometry Optimization Tolerance: 10E-4[Ha] = 0.0272[eV]

Orbital Cutoff: 3.4[Å] with an all electron core treatment and a double numeric (DN) basis set (5).

Convergence Testing:

The k-point mesh and SCF tolerance used in these calculations lead to unstable results that do not provide very accurate energies. Unfortunately, the computational resources available for this project were extremely limited and in order for the calculations to reach any result at all these extremely coarse parameters had to be used. The instability was confirmed by comparing the calculations against partial results with a SCF tolerance of 10E-6[eV] and a 4x4x1 k-point mesh. The system energies between the coarse calculations and the finer calculations varied by nearly 0.19[eV].

Calculations:

The reactant and product states of the reaction were defined as the hcp and fcc sites respectively.

A 3-layered 1×1 slab of Cu(111) with a vacuum spacing of 10[Å] was used for the calculations. An Ag atom was placed at the hcp site and the geometry of the Ag and all Cu layers were optimized in order to define the reactant state. An Ag atom was placed at the fcc site and the geometry of the system was optimized in order to define the product state. With more computational resources, this cell should probably be increased in size, with more unit cells and more layers, fixing the geometry of the bottom layers to represent the bulk material.

Figure 2: Reactant and product vacuum slabs used for transition state searching.

The transition state of the reaction was found using the complete LST/QST method as implemented in Dmol3 (6).

The energy of each image of the reaction was defined as follows:

E = E_{Ag|Cu(111)} - E_{Cu(111)} - E_{Ag}

Figure 3: Complete LST/QST transition state search.

The transition state was found to be the bridge site between the two threefold sites. While the exact numbers of the frequencies and energy at the transition state will not be very accurate, there was one imaginary frequency along the reaction coordinate that supports its identity as a transition state.

Figure 4: Reactant state, transition state, and product state structures.

The energy barrier found for the hcp to fcc transition was 0.051[eV], but the level of convergence for these calculations is on the same order of magnitude as this barrier, so the exact barrier cannot be determined from this study. The height of this barrier reported in literature is 0.023[eV] (7), giving the same order of magnitude as this study.

The energy of the two threefold sites was found to be slightly different, with the fcc site being 0.003[eV] more stable, however the calculations from this study cannot tell which one is actually more stable.

Conclusion:

Due to computational resource limitations, an exact value for the activation energy of Ag diffusing between adjacent threefold sites on Cu(111) could not be found, however the order of magnitude was determined to be 10E-3[eV]. The location of the transition state was found to be the bridge site in between the threefold sites, this location is supported by the presence of an imaginary frequency along the reaction coordinate. It was also found that the energy of the threefold sites is different, due to the presence(hcp) or absence(fcc) of a Cu atom directly below the site in the next layer, but it could not be determined which of these sites is more energetically favorable.

References:

(1)     M. Ohring (2001) “Materials Science of thin Films”, Section 7.4: “Kinetic Processes in Nucleation and Growth”, 386-400

(2)     B. Delley (2000) “From molecules to solids with the DMol3 approach” J. Chem. Phys. 113

(3)     J. Perdew, Y. Wang (1992). “Accurate and simple analytic representation of the electron-gas correlation energy”. Phys. Rev. B. 45 (23): 13244–13249.

(4)     H. Monkhorst, J. Pack (1976) “Special points for Brillouin-zone integrations.” Physical Review B, 13(12): 5188-5192

(5)     B. Delley (1990). “An All-Electron Numerical Method for Solving the Local Density Functional for Polyatomic Molecules”. J. Chem. Phys. 92: 508–517.

(6)     T. Halgren, W. Lipscomb (1977) “The synchronous-transit method for determining reaction pathways and locating molecular transition states.” Chem. Phys. Letters 49: 225-232

(7)     A. Kotri, E. El Koraychy, M. Mazroui, Y. Boughaleb (2017) “Static investigation of adsorption and hetero‐diffusion of copper, silver, and gold adatoms on the (111) surface” Surface and Interface Analysis, 49 (8)

Calculation of Activation Energy and Frequencies for Ag Diffusion between Adjacent Threefold Sites on Cu (111) Surface.

Author : Shyam Deo

Introduction

In order to understand formation of nanostructures, it is necessary to understand the underlying microscopic processes, such as adsorption, desorption, and diffusion. In the current study, we focus on the dynamics of a single adatom of silver (Ag) for diffusion between adjacent three-fold sites on Cu (111) surface. The activation energy for diffusion has also been calculated along with the frequencies of diffusion.

Methods

Plane-wave Density Functional Theory (DFT) calculations were carried out using the Vienna Ab Initio Simulation Package (VASP), version 5.4.4. The electron-electron exchange and correlation energies were computed using the Perdew, Burke, and Ernzerhof functional (PBE) [1]. The projector augmented-wave (PAW) [2] method was used to represent the ion-core electron interactions. For Cu, the cutoff radius for the above psuedopotential is 2.3 Bohr (1.22 Å) with 11 valence electrons in the following configuration 3d10 4s1. For Ag, the cutoff radius for the above pseudopotential is 2.4 Bohr (1.27 Å) with 11 valence electrons in the following configuration 4d10 5s1. The structural convergence criteria were 0.05 eV Å-1 for all unconstrained atoms, while the convergence criteria defining self-consistency of the electron density was 10-5 eV.

Figure 1: (a) Cu (111) 3×3 unit cell with 5 layers and vacuum. (b) Convergence test for k-points. (c) Convergence test for ENCUT

A Cu FCC experimental bulk lattice constant of 3.615 Å was used for building Cu (111) surface [3]. A five layered 3 x 3 unit slab of Cu (111) was used for surface calculations where the bottom three layers were frozen during the optimization and the top two layers were unconstrained (Figure 1a). A plane wave energy cutoff of 450 eV and Monkhorst-Pack [4] k-point mesh of 5 x 5 x 1 was used for all surface slab calculations after proper convergence tests described in the next section. To minimize spurious interslab dipole interactions between the periodic slabs, a vacuum space of 15 Å was used and dipole corrections were added in the direction perpendicular to the surface. Transition state was located using the climbing image nudged elastic band (CI-NEB) method [5] (5 images) by relaxing the reaction tangent force below 0.05 eV Å-1.Transition state was verified to contain a single imaginary frequency along the reaction coordinate. These vibration frequencies were calculated by the tag IBRION = 7 implemented in VASP. It determines the Hessian matrix (matrix of the second derivatives of the energy with respect to the atomic positions) using density functional perturbation theory [6-7].

Convergence Tests 

Convergence test for K-points mesh sampling the brilliouin zone was done at a plane wave energy cut off of 550 eV by varying the k-points against the single point energy of Cu (111) surface. The relative energies are plotted against the number of irreducible k-points as shown in Figure 1b. The energy varied by an amount 0.006 eV between 13 and 24 irreducible K-points. However, the energy fluctuated when the K-points were raised further. To reduce the computation cost, lower K-points were chosen due to limited availability of computational resources which should be a reasonable starting point because we are mostly concerned with a difference in energies between different systems and not the absolute numbers themselves. Then for ENCUT convergence, the cut-off energy for plane wave was varied from 350 eV to 550 eV as plotted in Figure 1c. The energy changed only by 0.017 eV between 450 eV and 550 eV. Finally, the k-points mesh of 5x5x1 with 13 irreducible K-points and a plane wave energy cut off of 450 eV was considered suitable for all subsequent surface calculations of adsorption energies.

Results and Discussions

Adsorption of Ag Atoms over Cu (111)

A Cu (111) surface has two types of three-fold binding sites available for adsorption by a silver atom – hollow-fcc and hollow-hcp sites shown in Figure 2a. The figure shows adjacent hollow-fcc and hollow-hcp sites. A hollow fcc site is a three fold site with no atom directly below  but with an atom in the layer next to it. A hollow hcp site is a three fold site with an atom directly below  and  with no atom in the layer next to it. We probed the adsorption of silver atom over 3×3 Cu (111) unit cell over these two adjacent sites.

The adsorption energy for Ag was calculated as:

Eads,Ag = EAg* – EAg – E*

where Eads,Ag is the adsorption energy of Ag, EAg* is the energy of Ag bound to the surface, EAg is the energy of single Ag atom and E* is the energy of the bare Cu(111) surface. The adsorption energies for the two adjacent three-fold sites are labelled in Figure 2b. Both sites show the same adsorption energy of 2.38 eV (Figure 2b and 2c). These two states obtained were then studied for diffusion between the two adjacent sites and the activation energy was probed for the same.

Figure 2: a) Cu (111) surface showing two adjacent hollow-fcc and hollow-hcp sites of adsorption, b) Ag atom adsorbed on a hollow-fcc site on Cu (111) and c) Ag atom adsorbed on a hollow-hcp site on Cu (111). Adsorption energies for the two states of adsorption are also shown for comparison.

Activation Energy for Diffusion

The two states with adsorption of Ag over hollow-fcc site and hollow-hcp site was used as two end points for NEB calculation. The former was used as the initial state while the later was taken as the final state after diffusion from the hollow-fcc site to the hollow-hcp site. The intermediate images for NEB run were obtained using linear interpolation to give five images along the straight line trajectory for diffusion between the two adjacent sites. The energy profile for images along the reaction coordinate obtained from BEB run is shown in figure 3. The transition state is also shown. The activation energy obtained is 0.0328 eV for diffusion from hollow fcc site to hollow hcp site while the activation energy for reverse diffusion obtained is 0.0305 eV. A bridge site adsorption of Ag was obtained as the transition state (one with the highest energy). This was confirmed for single imaginary frequency (with value -1.89 THz) along the reaction coordinate (from hollow-fcc to hollow-hcp and back vibration). The other two vibration frequencies obtained were 4.12 THz (vibration towards the plane of surface and away) and 1.01 THz (vibration parallel to the plane of adsorption). A single imaginary frequency along the reaction coordinate confirms the transition state for diffusion from hollow-fcc to the adjacent hollow-hcp site. All vibrational frequencies are reported in Figure 4 (a).

Figure 3: Energy profiles computed using NEB (5 images) for diffusion of Ag atom from hollow-fcc site to hollow-hcp site of adsorption. The TS state is labelled and the transition state structure is shown.

The frequencies of hopping was then calculated by the following expression as a function of temperature (Figure 4) where ΔE is the activation energy for hopping, ν1 to νN are  the vibrational frequencies associated with the minimum initial state (hollow-fcc site adsorption) and ν†1 to ν†N-1 are  the real vibrational frequencies associated with the transition state. The rate of hopping has been plotted in Figure 4b as a function of temperature.

Figure 4: (a) Vibrational frequencies in THz for Ag atom adsorption over hollow fcc site and transition state (bridge site)  on Cu (111) and (b) hopping rate for Ag atom on Cu (111) as a function of temperature (K)

Conclusion

NEB calculations for diffusion from the hollow-fcc site adsorption of Ag to the adjacent hollow-hcp site was studied and an activation barrier of 0.0328 eV was obtained. These two sites showed same adsorption energies of 2.38 eV. As a further study, diffusion coefficient can be calculated by ab initio methods from the calculated activation energies and adsorption energies as done earlier by Minkowski et al. for Cu adatoms over Cu (111) surfaces [8]. The calculated adsorption energies and activation barrier for Ag adatom matches reasonably with the reported values of 2.43 eV and 0.023 eV, respectively [9]. The coefficients for surface diffusion and rate for hopping can help understand aggregation of atoms on surfaces to form clustered particles of metals on the catalyst support.

References

  1. Perdew, J.P. and Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy. Physical Review B, 1992. 45(23): p. 13244-13249.
  2. Blöchl, P.E., Projector augmented-wave method. Physical Review B, 1994. 50(24): p. 17953-17979.
  3. Straumanis, M.E. and L.S. Yu, Lattice parameters, densities, expansion coefficients and perfection of structure of Cu and of Cu-In [alpha] phase. Acta Crystallographica Section A, 1969. 25(6): p. 676-682.
  4. Monkhorst, H.J. and J.D. Pack, Special points for Brillouin-zone integrations. Physical Review B, 1976. 13(12): p. 5188-5192.
  5. Nudged-elastic band used to find reaction coordinates based on the free energy. The Journal of Chemical Physics, 2014. 140(7): p. 074109.
  6. S. Baroni, P. Giannozzi and A. Testa, Phys. Rev. Lett. 58, 1861 (1987).
  7. X. Gonze, Phys. Rev. A 52, 1086 (1995); X. Gonze, Phys. Rev. A 52, 1096 (1995).
  8. Mińkowski, M. and M.A. Załuska-Kotur, Diffusion of Cu adatoms and dimers on Cu(111) and Ag(111) surfaces. Surface Science, 2015. 642: p. 22-32.
  9. Kotri, A.El koraychy, E.Mazroui, M., and Boughaleb, Y. ( 2017Static investigation of adsorption and hetero‐diffusion of copper, silver, and gold adatoms on the (111) surfaceSurf. Interface Anal.49705– 711. doi: 10.1002/sia.6211.

 

Comparing Two Transition Processes for Platinum Adatom Diffusion on Pt (100)

By Nathan McKee

Introduction

This post examines the diffusion of a platinum adatom on the (100) surface of platinum. Platinum has an fcc structure, so the (100) surface has a fourfold hollow site to which a platinum adatom can bind. An adatom in this position has multiple ways to diffuse to an adjacent fourfold hollow site, two of which are considered in this post. One is the simple hopping of the adatom to an adjacent site. The other is the atomic exchange of the adatom with an atom on the top layer of the surface, such that the surface atom travels to the adjacent site while it is replaced by the adatom. These processes are illustrated in Figure 1. DFT calculations were performed to find the transition states for both of these processes, and to estimate the frequency of each. The DFT calculations were carried out with the plane-wave based code CASTEP. The GGA PBE functional was used1, as well as on-the-fly generated (OTFG) ultrasoft pseudopotentials2. These pseudopotentials include 32 valence electrons per platinum atom, in the 4f14 5s2 5p6 5d9 6s1 configuration, with a cutoff radius of 1.27 Å. The convergence tolerance was set at 2.0*10-4 eV per atom.

Figure 1: Illustrations of surface diffusion in the cases of a) hopping, and b) exchange. All atoms shown are Pt, with some being highlighted to indicate that they play the biggest part in the transition process.

Cell Construction

To represent the (100) surface of platinum, a vacuum slab was used, with a vacuum distance of 12 Å. The thickness of the slab was chosen to be 2 layers, so as to represent the surface sufficiently without surpassing the computational limits of the system. In choosing the x and y dimensions of the cell, the interaction between adatoms must be considered. The use of plane-wave DFT necessitates the construction of a periodic crystal, but the goal of this post is to examine the diffusion of only a single adatom. This requires the cell to be large enough that the adatoms in adjacent cells do not affect each other significantly. A 3 x 3 cell was chosen to this end. The resulting cell is shown in Figure 1.

Cutoff Energy and k Points

To ensure that the calculations converge properly, an analysis of the selection of k points and the cutoff energy was performed on the slab, including the adatom. The k points were chosen to be an NxNx1 Monkhorst-Pack grid3 of evenly spaced points in reciprocal space, as is conventional for slab models in which the “a” and “b” lattice constants are equal. This analysis revealed that using a 6x6x1 k point grid in conjunction with a cutoff energy of 350 eV yields the energy within 0.1 eV. While this uncertainty is substantial, it is less than the difference between the transition state energies of the diffusion processes being considered, as will be shown later.

Transition State Search

Before considering transition states, the product and the reactant had to be identified through geometry optimization. The bottom layer of the slab was held in place to represent the bulk platinum crystal, while the top layer and the adatom were allowed to relax and find the optimal configuration. The force tolerance was set to 0.05 eV/Å, and the energy tolerance was set to 2.0 * 10-4 eV per atom. In general, one would expect to have to perform this geometry optimization for both the reactant and the product, but in the case the reactant and product are identical except for a constant translation. Therefore, only the configuration of the reactant was optimized, and the product was created by manually shifting the reactant atoms along the x axis.

With the reactant and product complete, the search for transition states begins by matching corresponding atoms in the reactant and product configurations. For the hopping transition, the adatoms in the reactant and product were matched with each other, while the surface atoms were matched according to position so that they barely moved. For the atomic exchange transition, the adatom in the reactant was matched with one of the surface atoms nearest to the hopping path in the product. Likewise, that surface atom in the reactant was matched with the adatom in the product.

The algorithm used in the transition state searches was the complete LST/QST approach. This algorithm begins with a linear synchronous transit (LST)4, followed by a conjugate gradient mnimization. The state found in this way is then used as the intermediate state in a quadratic synchronous transit (QST)4. If necessary, the last two steps are repeated until a transition state is found.

Results

These calculations yield a transition state for each diffusion process. When the reactant energies are subtracted from the transition state energies, the result is called the activation energy, ΔE. The transition rate can then be calculated using the Arrhenius equation5:

where k is the transition rate, and the prefactor nu is the atomic vibration frequency.

The activation energies obtained from the calculations are 1.04 ± 0.10 eV and 0.82 ± 0.10 eV for the hopping transition and the exchange transition, respectively. Estimating the vibrational frequency to be 5E-12 Hz, predictions can be made for the transition rates at any temperature. At 300 K, for instance, the rates are calculated to be:

khop = 1.7E-29 Hz

kexch = 9.3E-26 Hz

From these results, it is apparent that a small difference in activation energy can have a large impact on transition rates. To predict transition rates more precisely, calculations with smaller uncertainties would need to be performed. Nonetheless, these calculations indicate with some confidence that the atomic exchange transition is more frequent than the hopping transition. This preference is supported experimentally6, albeit for a slightly different exchange process in which both atoms move along the same diagonal direction.

References

  1. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868 (1996).
  2. Gonze, X. & Finocchi, F. Pseudopotentials Plane Waves–Projector Augmented Waves: A Primer. Phys. Scr. 2004, 40 (2004).
  3. Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 13, 5188–5192 (1976).
  4. Halgren, T. A.; Lipscomb, W. N. Chem. Phys. Lett.49, 225 (1977)
  5. Arrhenius, S. A. (1889). “Über die Dissociationswärme und den Einfluß der Temperatur auf den Dissociationsgrad der Elektrolyte”. Z. Phys. Chem. 4: 96–116.
  6. G. L. Kellogg and P. J. Feibelman, Phys. Rev. Lett. 64 (1990) 3147.

Determination of the transition state and the activation energy for Ag atom diffusion on Cu111 surface via the adjacent threefold sites

Author : Junseok Kim

 

Introduction

Understanding a diffusion process of an adatom on a metal surface is very important in many catalytic fields, especially for nanocrystal growth. One way for studying the diffusion process is to find a transition state with corresponding activation energy. Regarding to this, we studied Ag atom diffusion on Cu111 surface in this project. Especially, we focused on the hopping of Ag between adjacent threefold sites (fcc and hcp hollow) on Cu111 surface. And, we used Density functional theory (DFT) for all energy calculations since it is a powerful tool for calculating the energy of an atomic surface or slab [1]. Lastly, we determined the transition state for the hopping of Ag atom and the activation energy.

 

Method

All DFT energy calculations were performed using the Vienna Ab initio Simulation Package (VASP) [2, 3]. For an exchange-correlation functional, the generalized gradient approximation (GGA) – Perdew Burke Ernzerhof (PBE) was used [4]. Projector augmented wave (PAW) was also employed to describe the interaction of ionic core and valance electron [5]. To generate the pseudopotential, 11 valance electrons (3d10 4s1 of the electron configuration) with 2.3 Bohr cutoff radius was used for Cu, while 11 valence electrons (4d10 5s1 of the electron configuration) with 2.4 Bohr cutoff radius was used for Ag.

We employed 3.615 Å from experiment [6] as the lattice constant of bulk Cu for building Cu111 surface. We used a four layered 2 x 2 unit slab for Cu111 surface with fixing two bottom layers and 15 Å of vacuum spacing during the surface optimization and energy calculation, and a single Ag atom was placed on the surface. Due to asymmetric slab, self-consistent dipole correction was included for all energy calculations.

The energy convergence test was carried for the bare Cu111 surface without Ag atom with respect to k-point grid and cut off energy. K-point grid was varied from 5x5x1 to 12x12x1 and cut off energy performed varying the cut off energy from 300 to 550 eV. The convergence test was continued until the difference with respect to the result at the highest k-point grid and cut off energy is less than 0.02 eV.

Then, the geometry optimization for Ag atom on adjacent threefold sites (fcc and hcp sites) of Cu111 was performed with tolerances of 0.01 eV/Å for the force and 10^-6 eV for the energy before determining the transition state. To find the transition state, we used the climbing image nudged elastic band (CI-NEB) method with the reaction tangent force less than 0.05 eV [7]. Also, we used linear interpolation to get 5 intermediate images. To confirm that the obtained transition state is valid, we calculated the vibration frequency for the transition state.

 

Result

Convergence test

      

Table 1. The relative total energy of the bare Cu111 surface varying k-point grid (Left) and cut off energy (Right)

Table 1 shows the convergence test for the bare Cu111 surface with respect to k-point grid (Left) and cut off energy (Right). From the convergence test, we obtained the optimal values of k-point grid (8 x 8 x 1) and cut off energy (450 eV), which are used for all calculations.

 

Geometry optimization of Ag atom for adjacent threefold sites of Cu111

 

Figure 1. The optimized configuration of Ag on fcc (Left) or hcp (Right) sites on Cu111 surface ; Cu (Orange) and Ag (Light blue).

After the geometry optimization, it is found that there is no significant difference for the Ag-Cu bond length of fcc (2.593 Å) and hcp (2.595 Å) sites as well as their total slab energy, -58.436 eV for fcc site and -58.434 eV for hcp site. Even though the energy difference for Ag on two adjacent sites exists, it is hard to conclude which one is more favorable than the other, since the difference is too small. We may need to calculate the energy for other sites that Ag atom could be also placed on, such as atop or bridge, to determine the most favorable site, which is not important for this project. We used these two optimized states as the end point of transition state for Ag atom diffusion on Cu111 surface.

 

Transition state and vibration frequencies

Figure 2. Reaction energy profiles using Cl-NEB for the hopping of Ag atom between the adjacent threefold site of surface, fcc and hcp sites.

Figure 3. The configuration (Bridge) of Ag atom on Cu111 surface at the transition state ; Cu (Orange) and Ag (Light blue).

Figure 2 shows the reaction energy profiles for the hopping of Ag atom between the adjacent fcc and hcp sites of Cu111 using Cl-NEB. From Cl-NEB calculation with 5 images along the reaction coordinate, we could find the prospective transition state that has the highest energy. To confirm it is the transition state, we calculated the vibration frequencies for the state, only allowing Ag atom to move during the calculation. A single imaginary frequency (-1.049 THz) was obtained with other two real vibration frequencies, 3.676 and 1.776 THz. As a result, we finally determined the transition state where Ag atom is placed on bridged site of Cu111 surface as shown in Figure 3.

It is found that the activation energy is 0.0568 eV for the hopping from fcc to hcp site, while the reverse activation energy for the hopping from hcp to fcc is 0.0545 eV. The activation energies are quite larger than the reported values (0.023 eV) [8], Shyam’s work (0.0328 eV) and Sharad’s work (0.035 eV). The main reason for these deviations seems because of using a different size of unit cell. We used a relatively smaller size of unit cell (2×2), which leads to a higher surface coverage of Ag than that of others. Even though the coverage effect on Ag atom diffusion on Cu111 surface has not been reported, a high coverage may negatively affect the hopping of Ag on Cu111 and thus lead to the increase in the activation energy [9]. With respect to this, it is worth to note that our value would be a good case for the comparison about the coverage effect on Ag atom diffusion on Cu111. Furthermore, it has been expected that Cu atom can diffuse fast on Cu111 even at room temperature since the activation energy for the self-diffusion (0.03 ~ 0.04 eV) has not much difference with the Boltzmann’s factor (kBT ~ 0.025 eV) at room temperature [10]. Regarding to this, we can also expect that Ag atom would easily diffuse as well as the self-diffusion on Cu111 at room temperature.

 

Conclusion

In this project, we optimized the adsorption of Ag atom on fcc and hcp sites of Cu111 surface with a very small total energy difference between two sites. And, we successfully determined the transition state for the diffusion of Ag atom on Cu111 surface with Cl-NEB method using DFT. Finally, the activation energy for Ag atom hopping from fcc to hcp site on Cu111 surface is found to be 0.0568 eV, which is not well matched to the reported value but is meaningful. Thus, we believe that this DFT calculation would contribute to understanding the diffusion process of Ag atom on Cu111 surface.

 

Reference

[1] “Density-Functional Theory of Atoms and Molecules”, R.G. Parr, W. Yang, Y. Weitao (1994), Oxford University Press.

[2] G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50.

[3] G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.

[4] Perdew, J. P; Burke, K; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865-3868

[5] P. E. Blo¨chl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.

[6] “Crystallography and Surface Structure: An Introduction for Surface Scientists and Nanoscientists. 2nd edition”, Hermann, K. Wiley-VCH: Berlin, 2016.

[7] G. Henkelman, B.P. Uberuaga, H. Jónsson, the Journal of Chemical Physics, 113 (2000) 9901-9904.

[8] Kotri, E. El koraychy, M. Mazroui, Y. Boughaleb, Surface and Interface Analysis. 2017, 49, 705–711

[9] Inderwildi et al, “Coverage dependence of oxygen decomposition and surface diffusion on
rhodium (111): A DFT study”, J. Chem. Phys. 122, 034710 (2005)

[10] Jian Wang et al, “Diffusion barriers on Cu surfaces and near steps”, 2004 Modelling Simul. Mater. Sci. Eng. 12 1209

HCN to CNH isomerization

Author: Naveen Agrawal

Introduction:  Nudged Elastic Band calculation is one of the methods to search for a transition state which considers a chain of images attached through elastic bands between reactant and product.[1] This chain should ideally represent a set of states along the reaction path after the NEB( Nudged Elastic Band) calculation has converged to the minimum energy path. However,we have to guess these set of images to initialize the search. One way to generate these first set of images is a linear interpolation of atomic coordinates between the reactant and product. This scheme provides a faster way to create intermediate states when the reaction itself does not involve a complex rearrangement of atoms, I am presenting a case of isomerization where using linear interpolation creates chemically unreasonably intermediate states and initializing the NEB calculation with these images would not lead to desired reaction trajectory. Specifically for the isomerization of HCN, the linear interpolation leads to unreasonable geometric configurations as described in figure 2 which would be unreasonably difficult or rather impossible to optimize to desired reaction path as described in figure 4.

 Methods: DFT (Density Functional Theory) calculations were performed within the Vienna Ab initio simulation package (VASP, version 5.4.4), using the periodic supercell approach [2]. The projector augmented wave (PAW) method was used for electron-ion interactions [3,4]. The Perdew-Burke-Ernzerhof functional described the electron-electron exchange and correlation energies [5]. A plane wave basis set is used with an energy cutoff of 450 eV. For geometry optimization, the convergence criteria of the forces acting on atoms were 0.02 eV Å-1, while the energy threshold-defining self-consistency of the electron density was 10-5 eV. A cubic unit cell with 20A dimension was used for periodic supercell calculations. Considering that it is a unimolecular rearrangement and the electronic states should be discreet and localized therefore only gamma point was included for the Brillouin zone sampling. Zero Point Vibrational energy (ZPVE) corrections were included in the electronic energy calculated through VASP. Transition state search was performed using the Climbing Image Nudge Elastic band Method (CI-NEB). CI-NEB is an improved version of NEB method which improves the definition of reaction coordinate near the highest energy image.[6]

Geometry optimization of Isomers and linear interpolation: HCN and CNH molecules are optimized using periodic supercell approach as mentioned in the previous section. The optimized structures are linear and the distinction between the two structures is the position of hydrogen atom as attached to Carbon or Nitrogen. The optimized structures as obtained are shown in Figure 1.

Figure 1: Reactant and Product state for HCN to CNH isomerization

A linear interpolation of geometric co-ordinates between the reactant and product state resulted in a series of chemically unreasonable intermediate states as shown in Figure 2. Linear interpolation prompts the transfer of hydrogen  from carbon to the nitrogen through the carbon-nitrogen bond itself which would basically optimize to dissociated atoms due to huge repulsion during the NEB calculation.

Figure 2: Reaction trajectory through linear interpolation from 00 to 05

Transition state search and confirmation: A more reasonable set of intermediate images are considered as presented in Figure 3 as an initial guess to the Climbing Image Nudge Elastic Band (CI-NEB) method. These images were created based on nature of the geometry for the transition state presented in the literature [7] allowing the reasonable transition of hydrogen through C-H bending.

Figure 3: Initial guess of intermediate images for NEB calculation

This sequence of intermediate states resulted into a converged NEB calculation with image 03 as the transition state as indicated by the change of sign of tangent forces along the reaction path and the converged normal forces perpendicular to the path for all the images to be less than 0.02 eV/A as can be seen in the attached snapshot of summary of VASP output attached in table below. Ideally the tangent forces should change sign across the transition state only once, however in the calculation summary the initial images show very small positive sign. This is possibly due to the close proximity of the image 01 and 02 to the reactant state, and the gradient of energy at such small energy differences (0.02 eV) is probably not accurate. However, it is still highly likely that indeed 03 is the transition state as it is the maxima along a reasonable reaction trajectory as indicated by figure 4 and table below. The relative energy of the highest energy image  is 1.99 eV (activation barrier) without ZPVE corrections. I also increased the number of intermediate images to seven while gradually bending the C-H bond to allow better sampling of reaction path and it was possible to see a slightly better transition to the steep barrier of 1.99 eV as can be seen in the attached animation in appendix.

Image numberEnergy relative to initial stateTangent forcesNormal forces
00
10.0228620.00040.008946
20.0225670.005150.007189
31.992472-0.0060650.011774
40.7722631.0239810.019798
50.616651

 

Figure 4: Converged NEB reaction trajectory

A transition state should only have one imaginary normal mode of vibration and should lead to the product or reactant state considered.To confirm that 03 indeed is a valid transition state, a vibrational analysis was performed using IBRION =5 in VASP (finite difference energy derivative with harmonic approximation).Vibrational analysis revealed a high frequency unstable normal mode (negative curvature of energy surface) for the 03 image as shown in the first animation below. There were two additional low frequency mode obtained from vibrational analysis which represented rotation of C-N bond and are very weakly associated with reaction coordinate of H shuttle through rotation of C-H bonds.  Therefore, the activation barrier and the transition state geometry are precise within the convergence bounds of calculation mentioned.

Conclusion: HCN to CNH isomerization involves shuttling of hydrogen atom through bending of C-H bond in the transition state as can be seen in Figure 5 and Figure 4. Linear interpolation between the reactant and the product state provides very unreasonable guess for the reaction trajectory because of the linear geometry of both product and reactant state. Allowing the C-H bending in initial guess images provides a well converged NEB calculation with 1.84 eV of kinetic barrier with ZPVE corrections included. Previous theoretical studies are in agreement with the barrier calculated. [8]The presence of first order saddle point was confirmed with vibrational analysis. 

References:

1. Henkelman, G., Jóhannesson, G., & Jónsson, H. (2002) (pp. 269-302). Springer, Dordrecht.

2. Kresse, G. and Furthmüller, J.Physical review B54(16), p.11169.

3. P. E. Blochl, Phys Rev B, 1994, 50, 17953-1797

4. G. Kresse.et.al Phys Rev B, 1999, 59, 1758-1775

5. J. P. Perdew. et.al , Phys Rev B, 1992, 46, 6671-6687.

6. Henkelman.et.al  , The Journal of chemical physics, 2000,113(22), 9901-9904.

7. Koch.et.al., The Journal of Physical Chemistry C, 2007,111(41), 15026-15033.

8.Gong.et.al, The Journal of chemical physics,2005, 122(14), 144311.

Appendix:

Calculating the Barrier of Vinyl Alcohol Tautomerization to Acetaldehyde

Tautomerization of Vinyl Alcohol

In this post we wish to determine the reaction energy, barrier, and transition state structure for the tautomerization reaction of vinyl alcohol converting to acetaldehyde using DFT. This reaction is shown below in Fig.1. Tautomers are constitutional isomers; having the same number of each atom, but with a different connectivity. Most often, as in this case, tautomerization reactions are intramolecular proton exchange reactions.

Fig. 1 Tautomerization reaction of vinyl alcohol (left) to acetaldehyde (right).

It is known that acetaldehyde is more stable than vinyl alcohol at ambient conditions [1]. Thus, we should expect to find that the reaction energy is negative (forward direction of Fig. 1), which we can verify after we have performed our calculations.

Setting Up the Calculation

We will use DMol3 to geometry optimize the reactant and product state as well as locate the transition state. In addition, we will calculate the vibrational modes and frequencies of each the reactant, product, and transition state for reasons we will return to later.

In order to begin a transition state search in DMol3, the user must supply initial geometries for the reactant and product. With DMol3, the transition state procedure includes a routine for optimizing the reactant and product states before finding the TS, so these initial structures don’t need to be very accurate. Below are figures showing the initial reactant and product structures used for this example.

Fig. 2 Initial geometries of reactant (left) and product (right) to DMol3.

As TS searches follow atoms through spatial translations, it is important to tell DMol3 which atoms from the reactant and product are “equivalent”. This can be done using the reaction preview tool. Once a reaction preview has been made between the initial and final states, we can enter the calculation parameters we wish to use.

Calculation Details

Calculations were performed using the DMol3 Module as available in Materials Studio. The TS search method used was the Complete LST/QST Method with conjugate gradient optimization of reactants and products using a RMS convergence of 0.00185 Ha/Å and a maxinum of 5 QST steps. Exchange and correlation were treated within the generalized gradient approximation using the PBE functional with the DNP basis set and explicit core electrons (all electron calculation) [2-5].  The SCF tolerance was set to 1.0E-5 Ha with linear smearing using a value of 0.005 Ha.

Results and Discussion

First we present the optimized structures and compare the geometries to experimental values available from NIST [6,7].

Fig. 3 Optimized reactant and product from DMol3.

Table 2 Calculated and Experimental Geometries of Vinyl Alcohol
CalculatedExperimental
C1-C21.3361.326
C2-O31.3711.372
C2-H71.0911.097
C1-H51.0881.079
C1-H61.0931.086
O3-H40.9740.960
C1-C2-O3126.7126.2
C1-C2-H7122.9129.1
C2-C1-H6122.0121.7
C2-C1-H5119.9119.5
C2-O3-H4108.6108.3
CalculatedExperimental
C1-C21.5021.501
C2-O31.2181.216
C1-H41.0971.086
C2-H71.1211.114
C1-C2-O3124.7123.9
H4-C1-H6110.3108.3
C1-C2-H7114.9117.5

Below is a plot showing the evolution of the TS search. This plot includes the LST, QST, and GC (conjugate gradient) paths along the reaction coordinate. Note the transition state is also labeled on this plot–indicating a transition state has been found with the calculation parameters.

Fig. 4 DMol3 TS Search Energy Trajectories

DMol3 also outputs the geometry of the TS. One should confirm that the TS found (by computational means) matches intuition about how the reaction proceeds. Below the TS is shown with atom labels.

Fig. 5 TS structure found using DMol3.

The transition state structure is what we might expect: somewhere in between vinyl alcohol and acetaldehyde. Here we can see that at the TS, the C-C double bond hasn’t quite broken while the O-H bond hasn’t quite formed. This leads us to believe we may have found the correct transition state, however, we must verify by examining the vibrational modes/frequencies of the TS. Transition states occur at saddle points in the potential energy surface (PES) and, therefore, should have one negative (imaginary) frequency corresponding the the vibrational mode along the reaction path.

The vibrational analysis tool can be used on the TS to verify these requirements are met. Indeed, we find the TS has one negative vibrational frequency: -2000. To verify this mode is along the reaction path we may animate it and confirm visually. This animation is reproduced below as a series of frames.

Fig. 6 Imaginary mode for the TS found using DMol3.

Thus we can be confident that we have found a real TS corresponding to the tautomerization of vinyl alcohol. Of special importance are the reaction energy and reaction barrier which DMol3 calculates during the search procedure. The energy for this reaction was -11.2 kcal/mol while the barrier was 51.9 kcal/mol. Note that our prediction of a negative reaction energy was correct. We also notice this reaction has a relatively high barrier (inaccessible at room temperature).

So, we have found a TS structure for the reaction of vinyl alcohol to acetaldehyde and have calculated the reaction energy and barrier. Next we need to compare the structure and energies with literature values.

Juan Andres et. al. did a comprehensive ab inito study on this very reaction and their results are readily available online [8]. In their results they report that DFT methods found the reaction energy to be in the range 8.5-14.3 kcal/mol while the barrier was found to fall within the range 87.9-93.5 kcal/mol. Thus there is great agreement with our reaction energy but the barriers are in severe disagreement. This, however, is most likely a symptom of basis set selection and the differences in DFT methods for finding TS between our attempt and those of Andres.

To verify that our structures are the same we reproduce their table of geometry data and include our own. The labels in the tables below are consisted with Fig. 4.

Table 1 TS Geometry
This WorkLiterature
C1-C21.407
1.428
C1-O32.2282.299
C1-H41.5001.541
O3-H41.3111.305
C1-C2-O3111.093110.500
H4-C1-C266.73767.200
C2-O3-H475.82778.600
C1-C2-H7130.421131.100
H5-C1-C2-O3155.067147.4
H6-C1-C2-O363.97677.400

From the table above we see there is great agreement in the TS geometry, with the greatest errors showing in the dihedral angles. This leads us to conclude that we have indeed found a TS and that the departures from literature values (especially the reaction barrier) are purely from differences in the level of theory/basis set used in the calculation.

References

[1]      R.D. Johnson III. “CCCBDB NIST Standard Reference Database”. Retrieved 2014-08-30.

[2]     B. Delley, J. Chem. Phys. 92, 508 (1990).

[3]     
B. Delley, J. Chem. Phys. 113, 7756 (2000).

[4]     
J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation                                 made simple. Phys. Rev. Lett., 77:3865, 1996.

[5]     J. P. Perdew, K. Burke, and M. Ernzerhof. Erratum: Generalized gradient                                             approximation made simple. Phys. Rev. Lett., 78:1396, 1997.

[6]    Kuchitsu (ed.), Landolt-Bornstein: Group II: Atomic and Molecular Physics Volume 15:                     Structure Data of Free Polyatomic Molecules. Springer-Verlag, Berlin, 1987.

[7]    H. Hollenstien, Hs. H. Gunthard, “Solid State and gas infrared spectra and normal                           coordinate analysis of 5 isotopic species of acetaldehyde” Spec. Acta 27A, 2027 (1971)

[8]      Andrés, J. , Domingo, L. R., Picher, M. T. and Safont, V. S. (1998), Comparative                                   theoretical study of transition structures, barrier heights, and reaction energies for the                 intramolecular tautomerization in acetaldehyde/vinyl alcohol and                                                       acetaldimine/vinylamine systems. Int. J. Quantum Chem., 66: 9-24.