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.

Author: Hoang (Bolton) Tran

Problem Statement

Previously (Post 2), the energetic and tilt angle of hydroxyl (OH) adsorption on Pt (111) was studied. In this study (Post 3), the OH bond stretching vibration \((\nu_{OH})\) is of primary interest. Firstly, the zero-point energy (ZPE) will be considered for adsorption energy. Secondly, the adsorption energy \((E_{ad})\) and \(\nu_{OH}\) will be computed and analyzed at different OH adsorption environments:

• Adsorption sites: Atop vs. bridge (see Fig. 1 in Post 2)
• Surface coverage: 1.0 ML vs. 0.25 ML (Monolayer – number of adatom per surface atom)
• Water solvation: introduce water molecules above the adsorbed molecule.

Methodology

Software:

DFT was implemented in both CASTEP and DMol3, a plane wave basis set software embedded inside BIOVIA Material Studio (MS). [1]

Calculation:

Geometric optimization (CASTEP): Please see Post 2

Vibrational analysis (DMol3):

• Exchange-Correlation Functional: Generalized Gradient approximations – Perdew-Burke-Ernzerhof (GGA-PBE) [2]
• Self-Consistent Field tolerance: 1E-05 eV
• Monkhost-Pack [3] grid size: 10x10x1 (50 irreducible kpoints) for 1×1 supercell & 5x5x1 (13 irreducible kpoints) for 2×2 supercell.
• Self-consistent dipole correction in z-direction (normal to surface)

The procedure of optimizing, cleaving and adsorbing OH onto Pt(111) surface are the same as outlined in Post 2. This procedure was carried out again for each adsorption environments outlined below. The optimized structure of each adsorption environment was then used as input for the vibrational frequency calculation.

ZPE correction:

The adsorption energy is defined as:

\begin{equation} E_{ad}=E_{OH/slab}-E_{OH}-E_{slab} \end{equation}

To calculate the ZPE correction for \(E_{ad}\), ZPE was calculated for both the free OH radical and the adsorbed OH. Since the free OH only has O-H stretching vibration, the ZPE of adsorbed OH also considered only the O-H stretching \(\nu_{OH}\).

Adsorption environments:

1.  Adsorption sites:

The two adsorption sites of interest are the high symmetry atop and bridge sites. Values of \(E_{ad}\) at 1.0 ML coverage were already calculated from previous study. See Post 2 for more details.

2. Surface coverage:

Adsorption energy at the atop site for two difference OH coverage values are considered: 1.0 ML and 0.25 ML.  In Post 2, \(E_{ad}\) was already calculated for 1.0 ML. To achieve 0.25 ML, 2×2 supercell with one absorbed OH was built (Fig. 1).

Figure 1. Visualization of OH surface coverage of 1.0 ML and 0.25 ML (click on tab to view)

3. Water solvation:

Solvation effect, or the interaction of water molecules near surface with the adsorbed OH was explored. Atop sites at 0.25 ML was used with 1 or 2 water molecules floating just above the adsorbed OH (Fig. 2).

Figure 2. Optimized structure of water floating above adsorbed OH (click on tab to view)

The adsorption energy for the system with extra water molecules were calculated differently:

\begin{equation} E_{ad}=E_{OH/slab/water}-E_{OH}-E_{slab/water} \end{equation}

In this case, \(E_{slab/water}\) was calculated by optimizing a bare slab with floating water(s).

Results and discussion

1. ZPE corrections:

The adsorption energies, calculated across different adsorption environments, with and without ZPE are shown in table 1 below.

Table 1. Comparing \(E_{ad}\) with and without ZPE:

The difference in ZPE are all 0.01 eV or less, which is a lot smaller than the difference between \(E_{ad}\) at different environment. Therefore, it is concluded that ZPE does not affect the adsorption energy of OH on Pt (111) significantly.

This result also has significant implementation later.

2. Adsorption environments:

Calculated values of \(\nu_{OH}\) and \(E_{ad}\) for each environments are shown in Figure 3 and 4 below:

Figure 3. Frequency of OH stretching at different adsorption environments and of free OH.

Figure 4. Adsorption energy of OH at different environments.

The hypothesis of this study is: The adsorption energy is stronger (more negative) when the adsorption environment (site, coverage, solvation) allows the adsorbed OH to form stronger hydrogen bonds.

The strength of hydrogen bond is characterized by \(\nu_{OH}\), weaker OH stretching (lower \(\nu_{OH}\)) means stronger hydrogen bonds. The presented results put this hypothesis to the test.

Adsorption sites:

Looking first at the 1.0 ML (blue section) in Fig. 3, it is observed that the OH vibrations at atop and bridge sites are both smaller than that of free OH radical, indicating some degree of hydrogen bonds formed at the sites. The hydrogen bond strength is greater at the atop site. Looking now at 1.0 ML section in Fig. 4, atop adsorption is stronger than that of the bridge site (results of Post 2).

Surface coverage:

Results at “1.0 ML Atop” and “0.25 ML Atop 0w” represent the direct comparison for surface coverage. \(\nu_{OH}\) at 0.25 ML is about the same as free OH, indicating the absence hydrogen bonds. Fig. 4 then shows that OH adsorbs weaker at 0.25 ML.

Water solvation:

Focusing now at the 0.25 ML sections (orange sections). Surprisingly, adding one water (Atop 1w) does not lower \(\nu_{OH}\), indicating that there are still no hydrogen bonding. This, however, makes sense by observing the position of this water molecule in Fig. 2 above. The placement of this water molecule did not seem to produce hydrogen bonding, which would be water’s O atom pointing toward OH’s H atom.

Nevertheless, OH adsorption is still slightly stronger with one water. This can be attributed to the solvation effect on the slab itself, because this effect is the only difference left between the two equations for \(E_{ad}\) (\(E_{slab}\) and \(E_{slab/water}\)).

Considering now the case with two water, \(\nu_{OH}\) is now lower, indicating formation of hydrogen bonding, which is also confirmed by looking at position of the second added water in Fig. 2. Finally, the adsorption with 2 waters is much stronger than with both 0 water and 1 water by looking at Fig. 4.

Conclusion

Despite the above analysis seems to indicate that stronger hydrogen bonds lead to stronger adsorption, the aforementioned hypothesis is false. The most obvious evidence is from the ZPE correction. It was already shown that ZPE does not contribute significantly to the adsorption energy. The ZPE correction only comes from the difference in frequency of OH stretching \(\nu_{OH}\). Therefore, \(\nu_{OH}\), despite varying between different environments, does not contribute significantly to \(E_{ad}\).

In other words, the most that one can conclude from these results is that there is a correlation between hydrogen bond strength \(\nu_{OH}\) and adsorption energy \(E_{ad}\). The reason behind this correlation can be thought of as hydrogen bonds formation lowering the energy of adsorbed OH \(E_{OH/slab}\) or \(E_{OH/slab/water}\), which subsequently lowers adsorption energy \(E_{ad}\). However, it would be wrong to conclude that \(\nu_{OH}\) is solely responsible for the difference in \(E_{ad}\), because the magnitude does not check out.

Finally, another interesting takeaway is that water solvents increase the adsorption tendency, not merely by interacting directly with the adsorbed species (through hydrogen bonding), but potentially through other electronic interaction with the surface as well. This is why understanding of the solvation effects on solid-liquid interfacial processes (e.g. electrocatalysis) is of great interest.

References

[1] “First principles methods using CASTEP”, Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005) S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson, M. C. Payne

[2] J. P. Perdew, K. Burke, M. Enzerhof (1996). Generalized Gradient Approximation Made Simple, Phys, Rev. Lett., 77, 3865.

[3] H. J. Monkhorst and J. D. Pack (1976), Special points for Brillouin-zone integrations, Phys. Rev. B 13, 5188.