Author Archives: Stephen Holoviak

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

Leave a reply

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)

Evaluating Surface Energy Calculations of Pt(111) for Different Slab Model Parameters.

5 Replies

By Stephen Holoviak

This post will attempt to calculate the surface energy of the Pt(111) surface using a slab model and plane wave DFT calculations. Slab models were constructed using different cell parameters, such as number of layers, vacuum spacing and size of supercell. The surface energy for these models was calculated and compared to other models and experimental results.

Introduction:

The existence of a surface on a solid will increase the energy of the atoms at the surface when compared with those atoms in the bulk of the crystal due to the change in coordination atoms at the surface experience. Surfaces can be approximated in DFT calculations through the creation of slab models. The energy of these surfaces, \sigma, can be calculated with the following equation, where A is the total surface area of the slab, and E is the system energy of the slab and bulk respectively, and n is the number of atoms in the slab model.

\sigma = \frac{1}{A}[E_{slab} - nE_{bulk}] (1)

Experimentally, the surface energy of Pt(111) is found to be 0.159 [\frac{eV}{\r{A}}](2)

Since the experimental value of the surface energy is fairly small, all calculations were converged to a limit of 5E-3[eV], this limit is depicted on the graphs testing convergence by a red line.

Calculations:

Functional: GGA-PBE(3)

Pseudopotential: “On the fly generation” of ultrasoft psuedopotentals in CASTEP(4) for Pt with the following parameters:

  • Core Radius: 2.4[Bhor] ~=1.27[Å]
  • 32 electrons in valance with (4f14 5s2 5p6 5d9 6s1) configuration.

All k-points assigned with a Monkhorst-Pack method (5).

1-Layer Slab:

In an attempt to model a platinum surface, a slab model with 1 layer of Pt(111) atoms and 10[Å] of vacuum spacing was constructed.

Convergence with respect to the cutoff energy was calculated with a fixed (8x8x1) k-point mesh. The energy of the calculations were then compared to the previous calculation to find the difference. The level of convergence for the cutoff energy was within the previously described limit, which is off of the scale of the graph, indicated by the arrows under the line.

The calculations were found to be well converged for all cutoff energies above 450[eV].

Convergence with respect to the k-point mesh was calculated with a fixed energy cutoff of 450[eV]. The energy of the calculations were then compared to the previous calculation to find an energy difference.

A 1 atom slab of molecules is a very unstable system and not an ideal way to model a surface. Because of this instability changing the number of k-points included in the calculation always changed the energy of the system significantly. There was no consistent value of energy found with sufficient # of k-points.

The number of k-points in the vacuum direction was also increased, generating 8x8x2 and 8x8x4 meshes, to test the effect on the calculations, however there was no effect, changing the slab energy by less than 1E-4[eV].

The amount of vacuum space was also increased from 10[Å] to 20[Å] to test the effect on the calculation, however there was no effect , changing values of slab energy by less than 1E-4[eV].

5-Layer Slab:

A slab model with five layers of Pt(111) and 10[Å] of vacuum space was constructed. This model should represent a solid much better than the 1 layer model.

Convergence for number of k-points and cutoff energy was tested in the same way as the 1 layer model (see figures below). At cutoff of 550[eV] and a k-point mesh of 16x16x1 the calculations were found to be converged to within 5E-3[eV]. The energy of the slab was calculated with these parameters.

\sigma = (\frac{1}{53.3323[A^2]})(-65253.5044[eV] - 5*-13050.9532[eV] ) = 0.024[eV]

The energy of the bulk crystal was then calculated using the same cutoff and a 16x16x16 k-point mesh. Using the calculated energies, the area of the slab, and n = 5 for our slab, the surface energy was found to be 2.4E-2[eV]. This value is much smaller than the experimental value of 0.159[eV].

 

In order to find the source of this discrepancy, the calculations were repeated for a larger slab with 5 layers, spanning 2×2 and 3×3 unit cells. The surface energies calculated for these cells were similar to the original model, changing only on the order of 1E-3[eV] in final energy.

Changes to the amount of vacuum space were also made to the model and found to have no significant effect on the calculated surface energy, changing the value on the order of 1E-5.

Another attempt to improve the surface energy calculation for this model was to change the functional from the GGA PBE to the LDA CA-PZ(6, 7). It was reported in Vitos (8) that when compared to LDA, GGA predicts surface energies 7-16% lower, and when compared to experimental results predicts values up to 29% lower. The reasoning for this is described as the GGA functional underestimating the exchange energy. When the LDA functional was used for the calculations the surface energy was calculated as 3.4E-2[eV].

\sigma = (\frac{1}{53.3323[A^2]})(-65306.89886[eV] - 5*-13061.74095[eV] ) = 0.034[eV]

This value is a significant improvement over the GGA value, however it is still only a quarter of the experimental value.

8-Layer Slab:

The final attempt to improve the calculated surface energy was done by constructing an 8-layer slab with 10[Å] of vacuum space. The calculations were tested for convergence by the same method as the 1 and 5 layer models. The surface energy calculated with 550[eV] cutoff, 16x16x1 k-point mesh, and LDA CA-PZ functional was 3.4E-2 [eV]. This result is comparable to that of the 5-layer slab.

\sigma = (\frac{1}{53.3323[A^2]})(-104492.1244[eV] - 8*-13061.74095[eV] ) = 0.034[eV]

The vacuum space was changed to 20[Å] in order to test the effects on the calculations. An insignificant shift in energy on the order of 1E-4 [eV] was found.

The supercell size was also changed to 2×2 to test the effects on the calculations. The surface energies calculated for these cells were similar to the original model, changing only on the order of 1E-3[eV] in final energy.

Conclusion:

In order to calculate the surface energy of a metal a common approach is to use slab models. This post found 5-layer and 8-layer slab models resulted in essentially the surface energy. A 1-layer model was found to be a very unstable system and poor physical model. A 10[Å] vacuum spacing was found to be sufficient for all of the slabs tested, and increasing the vacuum spacing had no effect on the calculated energies. It was also demonstrated that the GGA PBE functional underestimates the surface energy more so than the LDA CA-PZ functional.

All models tested significantly underestimated the value of the surface energy when compared to the experimental results. The source of this error is currently unknown.

References:

(1) Sholl, D.S. Steckel, J.A. (2009) “Density Functional Theory: A Practical Introduction.” Wiley, 96-98.

(2) Miedema, A.R. (1978). “Surface energies of solid metals”. Zeitschrift fuer Metallkunde, 69(5), 287-292.

(3) Perdew, J. P. Burke, K. & Ernzerhof, M.  (1996) “Generalized Gradient Approximation Made Simple.” 3865–3868.

(4) Clark, S.J., et. al.(2005)”First principles methods using CASTEP”, Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570

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

(6) Ceperley, D. M. Alder, B. J. “Ground State of the Electron Gas by a Stochastic Method”, Phys. Rev. Lett., 45, 566-569 (1980)

(7) Perdew, J. P. Zunger, A. “Self-interaction correction to density-functional approximations for many-electron systems”, Phys. Rev. B, 23, 5048-5079 (1981).

(8) Vitos L. et. al. (1998) “The surface energy of metals.” Surface Science, 411(1-2), 186-202

By Stephen Holoviak.

This post will overview an attempt to predict the structure of platinum crystals. Calculations were made of the energy of Pt crystals in three different crystal structures: simple cubic (sc), face-centered cubic (fcc), and hexagonal close-packed (hcp). The lattice parameters of the crystals were varied and the number estimate 0[K] energy was calculated using the CASTEP[1] implementation in Materials Studio.

Calculations:

Functional: GGA-PBE[2]

Pseudopotential: OTFG ultrasoft

  • Core Radius: 2.4[Bhor] ~=1.27[Å]
  • 32 electrons in valance with (4f14 5s2 5p6 5d9 6s1) configuration.

Cutoff Energy: 350[eV]

# k-points: 8x8x8

Testing for Convergence:

Cutoff Energy:

To ensure the calculations were well converged, tests were run on the experimentally observed fcc structure. The default Pt structure was loaded into materials studio with an fcc structure and a lattice parameter of 3.9239[Å]. First, the cutoff energy was varied using the default 6x6x6 set of k-points.

The calculations were found to be well converged, with crystal energies within 0.05[eV] of each other, at cutoff energy of 350[eV] or higher. The number of k points being used in the calculation was then evaluated for convergence, using the 350[eV] cutoff energy and varying the number of k-points being used.

# of k-points:

The calculations were found to be well converged, within 0.05[eV] of calculated crystal energy of each other, for a minimum of 8x8x8 k-points.

Structure Calculations:

Simple Cubic:

In order to determine the lattice constant for the simple cubic structure, the energies of several initial estimates were found. The data from these estimates were then fitted to a 2nd-degree polynomial and the minimum value was found. Several more lattice parameters were tested on and around this minimum. These energies were then fitted to another 2nd-degree polynomial and another minimum was calculated. The energy of the structure was then calculated at this minimum. Since the energy differences between calculations after two iterations are well below the energy difference used to determine convergence the iterations were stopped here.

The minimum energy lattice parameter was found to be at a = 2.62[Å] with an energy of -13050.541[eV]. 

Face Centered Cubic:

Determining the lattice constant for the face-centered cubic crystal was very similar to finding the sc lattice constant. The energies of several initial estimates were found. The data from these estimates were then fitted to a 2nd order polynomial and the minimum value was found. Several more lattice parameters were tested on and around this minimum. These energies were then fitted to another 2nd-degree polynomial and another minimum was calculated. The energy of the structure was then calculated at this minimum. Since the energy differences between calculations after two iterations are well below the energy difference used to determine convergence the iterations were stopped here.

The minimum energy lattice parameter was found to be at a = 3.97[Å] with an energy of -13050.940[eV].

Hexagonal Close-Packed:

Determining the lattice constant for the hexagonal close-packed crystal had several key differences. The hcp crystal has two lattice constants that must be specified, a and c. The best packing to fill space is a ratio of c/a = 1.633[3], this ratio was fixed in the calculations and all changes were made in terms of the lattice parameter a. Also, the different lengths of the lattice constants mean that the number of k-points must be adjusted to keep a similar sampling density in k-space, for the hcp calculations the number of k-points was adjusted to 8x8x5. Another important adjustment is the fact that the primitive hcp unit cell has two atoms in it, where the sc and fcc cells were primitive, so the calculated 0[k] energies must be normalized per atom in order to compare them to the previous calculations. After these changes were made, the energies of several initial estimates were found. The data from these estimates were then fitted to a 2nd order polynomial and the minimum value was found. Several more lattice parameters were tested on and around this minimum. These energies were then fitted to another 3rd degree polynomial and another minimum was calculated. The energy of the structure was then calculated at this minimum. Since the energy differences between calculations after two iterations are well below the energy difference used to determine convergence the iterations were stopped here.

The minimum energy lattice parameters were found to be a = 2.81[Å], c = 4.59[Å], with an energy of -13050.87[eV]

Conclusion:

The lowest energy phase of platinum crystals was calculated to be the fcc crystal structure with a lattice constant of 3.97[Å] a summary of all of the calculated results can be seen in the table below.

Structure:Lattice Constant[Å]:Calculated 0[K] Energy[eV]:
Experimental (fcc)3.92N/A
Simple Cubic2.62-13050.54
Face Centered Cubic3.97-13050.94
Hexagonal Close-Packeda = 2.18
c = 4.59		-13050.87

The experimentally determined structure for platinum is the fcc crystal structure with a lattice constant of 3.92[Å][4], which is slightly smaller than the lattice constant calculated.

 

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.) “Generalized Gradient Approximation Made Simple.”, Perdew, J. P., Burke, K. & Ernzerhof, M. 3865–3868 (1996).

3.)”Density Functional Theory: A Practical Introduction” D. Sholl and J. Steckel, (Wiley 2009)

4.) “Materials Science and Engineering”,  Callister, (Wiley 1994)