by Brandon Bocklund

Introduction

Many of the observable properties that makes materials interesting are linked to the electronic structure, so it is desireable with to perform DFT calculations with the most accurate computational methods that are affordable. Calculations within DFT scale as \( \mathcal{O}(n^3) \) for n electrons [1], so many technologically important materials, including transition metals such as Fe, are computationally demanding to calculate when considering all the electrons. This is exacerbated for core electrons, since the wave functions near the nuclei tend to oscillate, requring a large number of plane waves required to describe the wave functions near the nuclear cores in extended systems. Since many key properties are influenced by the outermost valance electrons that control bonding, it is sensible to fix the core electrons and treat them approximately to make calculations of the same level of theory more affordable or to make more complex calulations affordable. Treating core electrons as fixed is known as the frozen core approximation.

In this post, different cutoff radii will be used to generate PAW datasets for Ti. Differences in the generated PAW datasets will be demonstrated by comparing the cutoff energy required to converge the energies of hcp Ti and by comparing the formation energy of TiO2.

Treating core electrons

Core and valence electrons

Typically there are two main parameters to consider when treating electrons as core or valence that effect how the wave functions will be constructed: the core radius and the number of core vs. valence electrons. The core radius is the spherical distance from the center of the atom that forms the outer bound for the approximation to the wave functions. The number of electrons in the core vs. valance is simply the partition of how many electrons of charge density an atom can contribute to bonding. For example, Ti has 22 electrons that must be partitioned between core and valance. Both of these parameters are typically chosen by the person developing the new description.

A basic description of some methods for treating core electrons in extended systems follows.

Pseudopotentials

Pseudopotentials replace the core wave functions with a pseudo wave function with the goal of defining a smooth function to represent the effective potential of the core. The pseudo wave function are constructed so that the value and the derivative of the pseudo wave function matches the value and derivative of the true wave functions at the core radius. The advantage of pseudopotentials is their formal simplicity, however the pseudo wave function prevents the recovery of the full core electron wave functions. James Goff published a post on pseudopotentials [2] that describes their construction in a detailed, introductory way.

LAPW

Another method for treating core electrons is the Linear Augmented Plane Wave (LAPW) method. The main idea of this method is that two partial waves, which can be used to describe the all-electron wave functions, are joined at the core radius and matched to the all-electron solution. Thus, the main advantage is the core wave functions can be recovered, but this method comes at a higher computational cost than the pseudopotential method. A more complete description of the pseuodpotential and LAPW methods can be found in [3].

PAW

The Projector Augmented Wave (PAW) method presented by Blöchl [4] as a formalism that generalizes and combines favorable aspects of the pseudopotential and LAPW methods through projector functions. The basis of the PAW method is the combination of partial wave functions from the isolated atom with pseudo partial waves. The PAW method defines a transformation operator that maps the pseudo wave functions onto the all-electron wave functions,
\begin{equation} \left| \Psi \right> = \mathcal{T} \left| \tilde\Psi \right> \end{equation}
for the all-electron wave function, \( \left| \Psi \right> \), the transformation operator, \( \mathcal{T} \), and the pseudo wave functions \( \left| \tilde\Psi \right> \). The transformation operator is defined as
\begin{equation} \mathcal{T} = 1+ \sum_i \left( \left| \phi \right> – \left| \tilde\phi \right> \right) \left< \tilde p_i \right|\end{equation}
where \( \left| \phi \right> \) and \( \left| \tilde\phi \right> \) are the all-electron partial waves and pseudo partial waves, respectively, and \( \left< \tilde p_i \right| \) are the projectors that are constructed to match the all-electron solution. In the PAW method, PAW datasets are generated that can recover the full core electron wave functions, like the LAPW method, but at a lower cost than the LAPW method, closer to the computational cost of pseudopotentials. More information can be found for the PAW method in [4-6] and more detailed description of the relationship between pseudopotentials, APW, and the PAW methods can be found in [4].

Calculation details

Calculations are performed using plane waves and the PBE functional implemented by Quantum Espresso [8] in the pwscf code. PAW datasets [4] are used to describe the core states. For Ti, PAW datastets are generated as part of this work with the ld1 code distributed with Quantum Espresso, based on the PS Library 1.0 parameterizations [9]. The PS Library 1.0 parameterization for O is used directly. The reference states for Ti and O are chosen as hcp Ti and the O2 molecule in a 10x10x10 fcc box. The size of the O2 box was converged to within 0.1 mRy. Detailed discussion of the convergence of the box size is out scope for this post.

The cutoff energy convergence of the different potentials are compared in the following sections below for Ti and TiO2. Oxygen cutoff energy was converged to 0.1 mRy at a wavefunction cutoff energy of 90 Ry and a kinetic energy cutoff of 720 Ry. For the calculation of formation energy, the maximum converged cutoff energy for each set of potentials is used across all structures. In all cases, Ti required a larger cutoff energy for convergence than O, so it was used for the calculations of TiO2.

The kpoint grids were converged to within 0.1 mRy. For hcp Ti, an 11x11x7 Monkhorst-Pack grid was used. For the O2 box, a 3x3x3 Monkhorst-Pack grid, and for rutile TiO2 a 4x4x7 grid. Showing the detailed kpoints convergence tests are out of scope for this study. Even though the different PAW datasets can lead to different electronic structures and kpoints convergence, a single set of kpoints for each set of structures were used across all of the different PAW datasets for simplicity.

For comparing the cutoff energy of hcp Ti, calculations with each PAW potential used the same structure, a hexagonal primitive cell with a lattice parameter of 5.8 a.u. and a c/a ratio of 1.58 (corresponding to the experimental c/a ratio). All static calculations were performed using tetrahedron integration and an energy convergence of 1e-8 Ry. For calculations of the formation energy, the cells and ionic positions for all cells (O, Ti for each potential, TiO2 for each potential) were relaxed with a force convergence criteria of 1e-9 Ry/bohr. The Brillouin zone integration for relaxations was treated with Gaussian smearing with a smearing width of 0.001 Ry.

Generating PAW datasets for Ti

The ld1 program was used to generate the PAW datasets. The documentation can be found online [10]. The starting point for the tests here is the input file for generating a PAW dataset for the spin polarized, scalar relativistic [12] Ti from the PS Library version 1.0.0 [9]. The input file for the PS Library starting point follows below. The settings here create a PAW dataset for Ti with Troullier-Martins pseudization [11] via Bessel functions. The the key inputs that will be changed is `rcore`, which is the matching radius in atomic units for the core charge smoothing. The bottom section describes the number of wave functions (6) that will be pseudized and which are described by the lines below. Taking the first line (3S) as the example, the first number is the principle quantum number (n=1, starting with 1 for the lowest s, 2 for the lowest p), then the angular momentum quantum number (l=0), followed by the number of occupying electrons, the energy used to pseudize the state (=0.00), the matching radius for the norm-conserving pseudopotential (=0.85), the matching radius for the ultrasoft pseudopotential (=1.30) and finally the total angular momentum (=0.0). Changing the norm-conserving or ultrasoft pseudopotential matching radii do not change the PAW dataset, since the radii are controlled by the input parameters above. Note that the energies for the unbound states (4p=4 Ry and unoccupied 3d=0.05 Ry) need to be non-zero.

&input
title='Ti',
zed=22.,
rel=1,
config='[Ar] 4s2 4p0 3d2',
iswitch=3,
dft='PBE'
/
&inputp
lpaw=.true.,
pseudotype=3,
file_pseudopw='Ti-PSL.UPF',
author='ADC',
lloc=-1,
rcloc=1.6,
which_augfun='PSQ',
rmatch_augfun_nc=.true.,
nlcc=.true.,
new_core_ps=.true.,
rcore=0.8,
tm=.true.
/
6
3S 1 0 2.00 0.00 0.85 1.30 0.0
4S 2 0 2.00 0.00 0.85 1.30 0.0
3P 2 1 6.00 0.00 0.80 1.60 0.0
4P 3 1 0.00 4.00 0.80 1.60 0.0
3D 3 2 2.00 0.00 0.80 1.50 0.0
3D 3 2 0.00 0.05 0.80 1.50 0.0

To demonstrate how the cutoff radius affects the calculation accuracy and expense (through the cutoff energy), two new PAW datasets will be generated with core radii of 0.4 a.u. (much smaller, further referred to as “small core”) and 2.0 a.u. (much larger, further referred to as “large core”). Note that these are note expected to be accurate, but instead to demonstrate how the parameters affect might affect calculations and impress the care and detail that should be given to generating proper potentials.

Energy cutoff convergence of Ti and formation energy TiO2

Figure 1a and 1b compre the relative and absolute energy convergence with respect to the wave function energy cutoff for the differently generated Ti PAW datasets. In Figure 1a, the PS Library converges at a wave function cutoff energy of 120 Ry, the large core at 110 Ry, and the small core does not converge to less than 0.1 mRy/atom even at 140 Ry. Computational demands and the instability of the PAW dataset limited the full convergence.

The cutoff energies follow the trend that would be expected. Compared to the default PS Libary PAW dataset, increasing the core radius allows for more of the wave function to be described by the core, and therefore the dataset is softer (requires smaller energy cutoffs). In addition, the small core is not converged even at 140 Ry because the smaller radius requires a harder potential with more plane waves being required to describe the wave function. Figure 1a also shows unusual behavior for the energy to not decrease variationally with the cutoff energy and instead shows that there was some osciallation in the energy with increasing energy cutoff. Note that the large core and the the PS Library default are quite close in energy, so the differences may be hard to see in both figures. Figure 1b aims to demonstrate that the energies for different pseudopotentials are different even for the same calculation settings, so they cannot be compared.

Figure 1a: Energy cutoff convergence relative to the energy calculated at a cutoff energy of 140 Ry. Energy differences are taken from 10 Ry to 140 Ry. Calculations not shown are outside the bounds of the plot.

Figure 1b: Energy convergence for cutoff energies ranging from 10 Ry to 140 Ry. Missing results at lower energy cutoff are off the scale of the plot.

As discussed previously, energies between different pseudopotentials cannot be directly compared, therefore a metric that can cancel the differences is desired. The formation energy of TiO2 can cancel the differences caused by the different treatment of Ti and give some insight into the practical implications of different treatments of the electronic structure. Using the cutoff energies found for Ti (120 Ry for PS Library, 140 Ry for small core, 110 Ry for the large core) the energies per atom of hcp Ti, O2 in a vaccuum box, and TiO2 are calculated. The formation energy of TiO2 is calculated in the usual way:
\begin{equation} E^\mathrm{form} = E_{\mathrm{TiO_2}} – (E_{\mathrm{Ti}} + E_{\mathrm{O_2}}) \end{equation}
Since in each case, the cutoff energies are convered for all calculations within each formation energy calculation, the formation energies can be safely compared even with different cutoff energies used between them.

The total energies and formation energies are shown in Table 1. The differences in the formation energy are larger than the convergence criteria (0.1 mRy/atom, 0.0001 Ry/atom) for the formation energies of the PS Library to the small core PAW dataset and the PS Library to the large core PAW dataset. The largest difference is 0.0026 Ry/atom for the small core and PS Library formation energies, which corresponds to above 3 kJ/mole-atom. Differences that large could impact the relative stability of TiO2 in a more comprehensive study that used inaccurate PAW datasets.

Conclusion

PAW datasets were reviewed in context of other types of frozen core approximations, then PAW datasets were generated for Ti with larger and smaller core radii than the ones parameterized in the PS Library v1.0.0. It was shown that the different core radii can affect the energy cutoff convergence, the total energy, and the formation energy of a compund, TiO2.

References

[1] Section 10.5, D.S. Sholl, J.A. Steckel, Density Functional Theory, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2009. doi:10.1002/9780470447710.
[2] J. Goff, Transferability of Cu Pseudopotentials in Cu/CuO Systems, Density Funct. Theory Pract. Course. (2019). https://sites.psu.edu/dftap/2019/03/31/transferability-of-cu-pseudopotentials-in-cucuo-systems/.
[3] D.J. Singh, L. Nordstrom, Planewaves, Pseudopotentials and the LAPW Method, Second Edi, Springer US, 2006. doi:10.1007/978-0-387-29684-5.
[4] P.E. Blöchl, Projector augmented-wave method, Phys. Rev. B. 50 (1994) 17953–17979. doi:10.1103/PhysRevB.50.17953.
[5] X. Gonze, F. Finocchi, Pseudopotentials Plane Waves-Projector Augmented Waves: A Primer, Phys. Scr. T109 (2004) 40. doi:10.1238/physica.topical.109a00040.
[6] R.M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, 2004. doi:10.1017/CBO9780511805769.
[7] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (1996) 3865–3868. doi:10.1103/PhysRevLett.77.3865.
[8] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials., J. Phys. Condens. Matter. 21 (2009) 395502. doi:10.1088/0953-8984/21/39/395502.
[9] A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Comput. Mater. Sci. 95 (2014) 337–350. doi:10.1016/j.commatsci.2014.07.043.
[10] http://web.mit.edu/espresso_v6.1/amd64_ubuntu1404/qe-6.1/Doc/INPUT_LD1.html
[11] N. Troullier, J.L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B. 43 (1991) 1993–2006. doi:10.1103/PhysRevB.43.1993.
[12] Koelling, D. D. & Harmon, B. N. A technique for relativistic spin-polarised calculations. J. Phys. C Solid State
Phys. 10, 3107–3114 (1977).

by Angela Nguyen

1. Introduction

The purpose of this project is to calculate the activation barrier for the diffusion of a Ag atom on a Cu(111) surface via adjacent three-fold hollow sites through the use of DFT with a plane wave basis set as implemented by the CAmbridge Serial Total Energy Package (CASTEP) module [1] in Materials Studio. To accomplish this task, the transition state for the diffusion must be determined and then verified via a vibrational frequency calculation and visual observation of the transition state. The vibrational frequency calculation was implemented using the Dmol module [2] in Materials Studio. The Cu(111) surface was chosen as it is the most stable surface of Cu. For the purpose of this project, the Ag atom will be diffusing from the fcc hollow site to the hcp hollow site on the Cu(111) surface.  From the DFT calculations, the transition state was determined to be where the Ag atom is crossing the bridge site on the Cu(111) surface. The activation barrier was then calculated to be 0.31 eV with an imaginary vibrational frequency at -31.37 \(cm^{-1}\).

2. DFT Parameters

Listed below are the parameters used for the optimization of the Cu(111) surface as implemented in CASTEP.

For Cu, the cutoff radius for the above psuedopotential is 2.0 Bohr (1.06 Å) with 11 valence electrons in the following configuration \(3d^{10} 4s^{1}\). For Au, the cutoff radius for the above pseudopotential is 1.59 Bohr (0.84 Å) with 19 valence electrons in the following configuration \(4s^{2} 4p^{6} 4d^{10} 5s^{1}\).

Below is a list of parameters for the vibrational calculations of the transition state as implemented in Dmol.

3. Methodology

3.1 Construction of Cu(111) Vacuum Slab and Convergence Check

A vacuum slab of Cu(111) consisting of 5 layers was constructed from bulk Cu with a vacuum spacing of 12 Å between each slab in which the bottom 3 layers of the slab was constrained to model the bulk. The construction of 5 layers for the vacuum slab was chosen to ensure that the bulk could be modeled by the bottom 3 layers. A more rigorous study in the future could be performed to determine whether varying the number of layers in the slab affects the activation energy. A vacuum spacing of 12 Å was used to ensure that the slabs were isolated from each other due to the periodic nature of DFT. Single point energy calculations were performed on the vacuum slab to ensure convergence with respect to k points and the cutoff energy. A kpoint mesh of 2x2x1 (2 irreducible kpoints) and an energy cutoff of 400 eV was chosen to be sufficient for convergence of future calculations. The Cu(111) surface was then geometrically optimized using the chosen kpoint grid and cutoff energy.

Figures 1 and 2 Figure 2 display plots of the cohesive energy per atom as a function of the number of irreducible k points and cutoff energy respectively as determined by Materials Studio. The black dotted lines represent the +/- 0.01 tolerance used to determine which kpoint grid and energy cutoff to use.

Figure 1. Cohesive energy (eV/atom) vs. Number of irreducible k points. The blue circles represent the data points obtained from DFT calculations in Materials Studio. The blue line is to help guide the eye across the plot. The black dashed lines represent the +/- 0.01 eV tolerance used to determine the optimal number of irreducible k point mesh. The red square represents optimal k point mesh that will be used in future calculations.

Figure 2. Cohesive energy (eV/atom) vs. Energy cutoff (eV). The blue circles represent the data points obtained from DFT calculations in Materials Studio. The blue line is to help guide the eye across the plot. The black dashed lines represent the +/- 0.01 eV tolerance used to determine the optimal energy cutoff. The red square represents optimal energy cutoff that will be used in future calculations.

Two different versions of the optimized Cu(111) surface was constructed with the Ag atom in either the fcc hollow site or the hcp hollow site. These will represent the initial and final states for the following transition state search respectively.

3.2 Transition State Search and Vibrational Frequency Calculation 

With the initial and final states for the diffusion of Ag on Cu(111) built, a linear interpolation of 10 images between the two states was performed. A transition search was then conducted using the above outlined CASTEP parameters using the “Complete LST/QST” algorithm [8] in Materials Studio. An option was chosen in Materials Studio that the initial and final states should be optimized before any transition state search calculation was performed. This was done because linear interpolation between optimized initial and final states failed. Once the transition state was found, a vibrational frequency calculation was performed on the transition state to determine if the transition state was reasonable or not. If it is reasonable, the vibrational calculation should result in one imaginary (negative) frequency in which Ag atom is moving along the reaction coordinate towards the hcp hollow site. In addition, visual inspection of the transition state was done to make sure it was reasonable.

4. Results

Figures 3 and 4 below depict the optimized initial and final states of Ag adsorbed onto Cu(111) on the fcc and hcp hollow site respectively.

Figure 3. Ag adsorbed onto the fcc hollow site of Cu(111). The light blue sphere represents Ag and the orange spheres represents Cu.

Figure 4. Ag adsorbed onto the hcp hollow site of Cu(111). The light blue sphere represents Ag and the orange spheres represent Cu.

Figure 5 shows the transition state for the diffusion of Ag on Cu(111) via adjacent threefold hollow sites. From the figure, it can be seen that the transition state for the diffusion of Ag across Cu(111) is located along the reaction coordinate where the Ag atom is crossing the bridge site. It is reasonable for the transition state to be located here, since the Ag atoms will have to maneuver over two Cu atoms before settling in the adjacent hcp hollow site. The activation energy was determined to be 0.31 eV.

Figure 5. Transition state for the diffusion of Ag on Cu(111). The Ag atom is located along the bridge site of the surface. The light blue sphere is Ag and the orange spheres are Cu.

Lastly, in order to confirm the above transition state, a vibrational calculation was performed, which resulted in an imaginary frequency at -31.37 \(cm^{-1}\). To verify that this was a plausible vibrational frequency for the transition state, a visualizer in Materials Studio was used to view this vibration, which can be seen in Figure 6. From the figure, it can be seen that the Ag atoms is moving along the reaction coordinate towards the hcp hollow site while crossing the bridge site. This confirms that a proper transition state was found for the diffusion of Ag on Cu(111).

The table below shows the activation energy and imaginary frequency as calculated in this post.

5. Conclusions

For the diffusion of Ag on Cu(111) via adjacent hollow sites, it was determined that the transition state is located at the bridge site with an activation energy of 0.31 eV. This was further confirmed through visual observation and by checking for one vibrational frequency for the transition state. One imaginary vibrational frequency was located at -31.31 \(cm^{-1}\) corresponding to the Ag atoms crossing the bridge site towards the adjacent hollow site.

From literature, Kotri et al. stuided the hetero-diffusion of metal atoms on other metal surfaces via adjacent hollow sites [9]. One such surface that they tested was the hopping of Ag on Cu(111) using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [10]. To determine the activation energy, Kotri et al. used the drag method in which the adatom is dragged from its inital site to the nearest empty site [9]. At each step, the adatom is then allowed to relax with the respect to the plane perpendicular to its motion. In their work, it was determined that the activation energy for Ag on Cu(111) was 0.023 eV, an order of magnitude less than the activation energy determined by this project. Interesting to note is that the activation barrier for the Ag atom hopping from the fcc hollow site to the hcp hollow site and vice versa was the same. For this project, the activation energy varied by 0.08 eV depending on which path the Ag atom took.

Reasons for the deviation between activation energies can be due to the difference in methodologies on how the transition state was obtained, how much “optimized” the initial and final states were, and the number of kpoints sampled. The “Complete LST/QST” algorithm could not be the most accurate method to use to find the transition state as it only optimizes the image with the highest energy, while leaving the other images untouched. Even though the transition state was “confirmed” to be correct through vibrational frequency analysis, it would be useful to run another transition state search method to verify if the current activation energy is correct. One such algorithm to use the nudged elastic band (NEB) method in which all images are optimized rather than the image with the maximum amount of energy in the “Complete LST/QST” method. When viewing the optimized inital and final images (Figures 3 and 4), it can be seen that the Ag atom is not located in the middle of each hollow site, where the highest symmetry of the site is. This could lead to higher energy values calculated as this would be the Ag would have to diffuse a longer distance between the fcc and hcp hollow sites. Lastly, a number of kpoints sampled in the Brillouin zone could be increased to sample any missed phenomenon in the system that was not detected before.

6. References

[1] Clark Stewart J et al., “First principles methods using CASTEP ,” Zeitschrift für Kristallographie – Crystalline Materials , vol. 220. p. 567, 2005.
[2] B. Delley, J. Chem. Phys. 92, 508 (1990).
[3] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett., vol. 77, no. 18, pp. 3865–3868, Oct. 1996.
[4] J. P. Perdew et al., “Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation,” Phys. Rev. B, vol. 46, no. 11, pp. 6671–6687, Sep. 1992.
[5] D. D. Koelling and B. N. Harmon, “A technique for relativistic spin-polarised calculations,” J. Phys. C Solid State Phys., vol. 10, no. 16, pp. 3107–3114, Aug. 1977.
[6] D. Vanderbilt, “Soft self-consistent pseudopotentials in a generalized eigenvalue formalism,” Phys. Rev. B, vol. 41, no. 11, pp. 7892–7895, Apr. 1990.
[7] H. J. Monkhorst and J. D. Pack, “Special points for Brillouin-zone integrations,” Phys. Rev. B, vol. 13, no. 12, pp. 5188–5192, Jun. 1976.
[8] Halgren, T. A.; Lipscomb, W. N. Chem. Phys. Lett., 49, 225 (1977).
[9] A. Kotri, E. El koraychy, M. Mazroui, and Y. Boughaleb, “Static investigation of adsorption and hetero-diffusion of copper, silver, and gold adatoms on the (111) surface,” Surf. Interface Anal., vol. 49, no. 8, pp. 705–711, Aug. 2017.
[10] LAMMPS User Manual 30, Sandia National Laboratories (Oct 2014) 2.

Author: Sharad Maheshwari

Introduction

Surface diffusion can play an important role in creating a topological modification on the surface as in the case of nanoparticle growth. In the following work we aim to identify the kinetic barrier for the diffusion of a Silver (Ag) atom as it hops from a three-fold fcc site to the adjacent three-fold hcp site on (111) surface of Copper (Cu). We will use density functional theory-based methods to evaluate the activation barrier for the hopping process. We will also calculate the frequency of this hopping process.

Calculation Details

Electronic Calculations

Electronic structure calculations were performed using the Vienna Ab initio Simulation Package (VASP)[1, 2] a plane wave basis set pseudo-potential code. We used the projector augmented wave (PAW) [3] method for core-valence treatment. The exchange and correlation energies were calculated using the Perdew, Burke, and Ernzerhof (PBE) [4] functional described within the generalized gradient approximation (GGA) [5]. A plane-wave basis set cutoff energy of 450 eV was used. The ionic convergence limit was set to 0.01 eV Å^-1 while the electronic convergence limit was set to 10^-5 eV. The Fermi level was smeared with the Methfessel-Paxton [6] scheme using a smearing width (σ) of 0.2 eV.

Surface Model

A 3 x 3 surface slab model was used to construct periodic surfaces of Cu (111) using an experimental bulk lattice constant of 2.556 Å. The slab models were comprised of 4 layers of metal atoms.  The top two layers of the slab along with the adsorbate were allowed to relax until the force on atoms in these layers was within 0.01 eV Å^-1 while the bottom two layers were kept fixed to imitate their bulk arrangement.  A vacuum region of 10 Å was inserted in the models to exclude periodic interaction between the slabs. Dipole corrections were also added in the direction normal to the surface.

The sampling of the Brillouin zone for all 3 x 3 surface cells was conducted with a k-point mesh of 5 x5 x 1 generated automatically using the Monkhorst Pack method[7]

Transition States

To locate the transition state, a series of 5 linearly interpolated images were formed between the reactant and the product state. The nudged elastic band method (NEB) [8] was used to locate the transition state. Once an approximate transition state was obtained, it was further refined using the climbing image nudged elastic band method (CI-NEB) [9]. Vibrational frequency calculations were then performed using IBRION = 7 in VASP to confirm the transition state found the first order saddle point. It was further ensured that the single imaginary frequency obtained was along the reaction pathway.

Results

Figure 1 illustrates the reactant, transition and product state for the hopping of Ag atom from the fcc to hcp site of Cu(111) surface. The activation barrier for the hopping process is 0.035 eV.

Figure 1. Top and Side view of the a) reactant b) transition c) product states for the Ag hopping from fcc site to the hcp site on Cu (111)

For the barrier evaluated using the CI-NEB method, the reactant and the product states were used to create a series of images. The energy profile obtained from the calculation for the images is shown below in Figure 2. (All energies are referenced to the reactant state and does not include zero- point vibrational energy corrections).

Figure 2. Reaction energy diagram for the hopping of Ag atom from fcc to hcp site on Cu(111)

One important thing that was observed in a trial calculation is that the NEB method missed the transition state when 4 linearly interpolated intermediate images were used as shown in Figure 3. Using small even number of images can result in missing the saddle point and thus if using the NEB method, a higher density of images are needed near the saddle point. This highlights the need for caution when evaluating the activation barrier using the NEB method and using linear interpolation to obtain the images.

Figure 3. Reaction energy diagram obtained from the calculations using NEB method with 4 images obtained through linear interpolation

On an absolute scale, the zero-point vibrational energy (ZPVE) corrections were found to be very small. However, as the activation barrier itself is very small, we corrected the activation barrier by including the ZPVE corrections to the initial and the transition state. The corrected activation barrier after including the ZPVE corrections is 0.027 eV.

To evaluate the frequency of this hopping, we use the following equation:

where υ_i and υ_i^†  are the frequencies evaluated for reactant and transition state respectively, k_B is the Boltzmann’s constant, ΔE_act is the evaluated activation barrier and T is the temperature in Kelvin.  As the transition state is first-order saddle point, it is energy minimum in all direction except for the reaction coordinate, along which it is maximum. Figure 4 shows the hopping frequency of Ag atom from fcc to hcp site on Cu(111) surface at different temperatures.

Figure 4. Frequency of hopping of Ag atom from fcc to hcp site on Cu(111) surface as a function of temperature

Conclusion

We used DFT methods to estimate the kinetic barrier for the hopping of Ag atom from the fcc to hcp site on Cu(111) surface. The activation barrier for this process after including the ZPVE corrections is found to be 0.027 eV. The value found is in reasonable agreement with the reported values [10]. This work illustrates the DFT methods can be used to evaluate kinetic parameters of the process. It also highlights that while using methods to evaluate the activation barrier, one should be careful and must ensure that the transition state found is the correct transition state. Vibrational calculations must be performed and a single imaginary frequency along the reaction coordinate must be confirmed to ensure that the transition state found is indeed the transition state we were looking for.

References

[1] G. Kresse, J. Furthmuller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 6 (1996) 15-50.

[2] G. Kresse, J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B, 54 (1996) 11169-11186.

[3] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B, 59 (1999) 1758-1775.

[4] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 77 (1996) 3865-3868.

[5] 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) 6671-6687.

[6] M. Methfessel, A.T. Paxton, High-precision sampling for Brillouin-zone integration in metals, Physical Review B, 40 (1989) 3616-3621.

[7] H.J. Monkhorst, J.D. Pack, Special Points for Brillouin-Zone Integrations, Phys. Rev. B, 13 (1976) 5188-5192.

[8] G. Henkelman, H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, The Journal of Chemical Physics, 113 (2000) 9978-9985.

[9]  G. Henkelman, B.P. Uberuaga, H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, The Journal of Chemical Physics, 113 (2000) 9901-9904.

[10] Kotri, E. El koraychy, M. Mazroui, Y. Boughaleb, Static investigation of adsorption and hetero-diffusion of copper, silver, and gold adatoms on the (111) surface, Surface and Interface Analysis, 49 (2017) 705-711.