# Solution to exercise 2

Contents

### 1 Aluminum surface and adsorbates

First, we have to find the lattice vectors and stacking sequence
of an fcc (111) surface.
Two vectors that span the *xy*-plane of an fcc surface can be found using
this drawing:

The two lattice vectors are:

e_{1}= (1/2^{1/2}, 0, 0)

e_{2}= (-1/8^{1/2}, (3/8)^{1/2}, 0)

The distance between two fcc planes is 1/3 of the body diagonal, therefore
it is *a*/3^{1/2}.

The coordination *Z* of a bulk atom is *Z* = 12 with six nearest
neighbors in the plane, three in the plane above and three in the
plane below. As the plane above is removed, a surface atom has a
coordination of *Z* = 9.

When adding the function `FCC111` to the script Build.py, we
have to think about the stacking sequence. We need a smart way to
make sure that the first atom sits at coordinates (0, 0) (in terms of **e**_{1} and **e**_{2}), the second atom at (2/3, 1/3),
the third atom at (1/3, 2/3), the fourth atom at (0, 0) and so on.
This can be done with the modulo operation, e.g. if *n* is the number
of the layer (starting from zero), the position of the atom is

[(2n% 3)/3, (n% 3)/3, -nd/h]

where *d* = *a*/3^{1/2} is the spacing of the layers and *h* is
the total height of the unit cell in the *z*-direction.

A comment about whether positions should be given in scaled (in terms of the unit cell vectors) or cartesian coordinates: If the coordinates are given before the unit cell is set in the ListOfAtoms, they have to be given in scaled coordinates. If the coordinates of atoms are given at a time, when the unit cell has been set, they have to be given in cartesian coordinates.

#### 1.1 Aluminium fcc(111) Surface Slabs

One can estimate the surface tension from a simple effective medium theory description

Aσ ≃ [1 - (Z/Z_{0})^{1/2}]E_{coh}

The coordination numbers are *Z* = 9 for the surface and *Z*_{0}
= 12 for the bulk. The cohesive energy is *E*_{coh} = *E*_{atom} - *E*_{B} = 3.34 eV. The surface area per surface
atom is calculated as the shaded area in the figure above. It is

A=a^{2}3^{1/2}/4 = 7.1 Å^{2}

With these values, one obtains σ = 63 meV/Å^{2}.

The definition of the surface tension is

E_{N}= 2Aσ +NE_{B}

where *E*_{N} is the total energy of a slab with *N* layers, *A*
the area of the surface unit cell (the factor 2 because the slab has
two surfaces). The limit *N* -> ∞ corresponds to the macroscopic
crystal termination.

Instead of calculating the bulk energy, one can calculate the surface
tension from energy differences of different slabs. This has the advantage
that there will be a certain error cancellation, as we use the same
parameters. One can calculate *E*_{N-1} by substituting *N* by
*N* - 1 in the equation above. By calculating *NE*_{N-1} and
(*N* - 1)*E*_{N} and subtracting them one can show

σ = (NE_{N-1}- (N- 1)E_{N}) / (2A)

Using a **k**-point setup of 6 x 6 x 1 one obtains the surface energies

N |
σ [meV/Å^{2}] |
---|---|

3 | 54.9 |

4 | 43.2 |

5 | 55.8 |

6 | 57.5 |

One can see that there is an anomaly for *N* = 4 but that the surface energy
converges for an increasing number of layers. The EMT estimate
turns out to be satisfied quite well.

#### 1.2 Adsorption of H on the Al(111) surface

On an fcc (111) surface, there are four inequivalent adsorption sites, as shown:

The two hollow sites differ in this way that in the layer below the hcp hollow site there is an atom directly below the hollow site, whereas for the fcc hollow site there is first an atom directly below the hollow site two layers further down.

In scaled (in terms of the unit vectors) coordinates, the four different adsorption sites and their coordination are:

site | coordinates | coordination |
---|---|---|

ontop | (0, 0) | 1 |

bridge | (1/2, 1/2) | 2 |

hcp hollow | (2/3, 1/3) | 3 |

fcc hollow | (1/3, 2/3) | 3 |

In the program, the position of the adsorbate has to be given in
cartesian coordinates. The *x*- and *y*-coordinates can be calculated
using the scaled positions and the unit vectors. For the
*z*-coordinate one has to find a good guess, e.g. from the covalent
radii *r*_{Al} = 1.18 Å and *r*_{H} = 0.37 Å. One
possibility is to assume the Al-H bonds to be the sum of the covalent
radii *d*_{AlH} = *r*_{Al} + *r*_{H}. Then the
*z*-coordinate of the adsorbate can be calculated in the following way

site | z |
---|---|

ontop | 1.55 Å |

brigde | (d_{AlH}^{2} - a^{2}/8)^{1/2} = 0.59 Å |

hollow | (d_{AlH}^{2} - a^{2}/6)^{1/2} |

Plugging the values in for the hollow site, one notices that one
obtains an imaginary number, i.e. even in the plane the H atom is
already further away from the Al atoms than for the covalent radii.
However, doing the simulations one notices that the H atom actually
ends up sitting ~ 1 Å above the surface. This shows that the rule
of thumb with the covalent radii is not always useful. A better
approximation for *r*_{Al}
is to use half the distance between the aluminium atoms in the bulk, which
is *a*/2/2^{1/2} = 1.43 Å. Adding this to *r*_{H} gives
*d*_{AlH} = 1.80 Å. With this it is

site | z |
---|---|

ontop | 1.80 Å |

brigde | (d_{AlH}^{2} - a^{2}/8)^{1/2} = 1.09 Å |

hollow | (d_{AlH}^{2} - a^{2}/6)^{1/2} = 0.71 Å |

which turns out to be better guesses.

Running these calculations with a **k**-point sampling of 8 x 8 x 1
(probably one can do with less) one gets the relative energies listed
here:

site | total energy [eV] | relative energy [eV] |
---|---|---|

ontop | -127.93 | 0 |

bridge | -128.05 | -0.12 |

hcp hollow | -128.02 | -0.09 |

fcc hollow | -128.11 | -0.17 |

The calculation was done with 340 eV plane wave cutoff, 400 eV density cutoff. One can see that the fcc hollow site is the most stable site, followed by the bridge site. In general, the energy differences between the different sites are small. Let us look at diffusion of H from one fcc site to a neighboring fcc site. The H atom will have to pass through a hcp site (transition state), which is 0.08 eV higher in energy. With this activation energy, we can use the Arrhenius equation to find rates:

k=Aexp[-E_{a}/(k_{B}T)]

with *A* as the prefactor and *E*_{a} as the activation energy.
The Boltzmann constant is *k*_{B} = 1.381 x 10^{-23} J/K.
For the different temperatures, we obtain the following rates:

*T* = 100 K: *k* = 9.3 x 10^{8}/sec

*T* = 200 K: *k* = 9.7 x 10^{10}/sec

*T* = 300 K: *k* = 4.5 x 10^{11}/sec

*T* = 500 K: *k* = 9.8 x 10^{11}/sec

For the rule of thumb for biological catalytic processes, let us
assume room temperature to be *T* = 293 K. The temperature increase
is Δ*T* = 10 K. We start from the Arrhenius equation and use that the ratio should be 2.

2exp[-

E_{a}/(k_{B}T)] = exp[-E_{a}/(k_{B}[T+ ΔT])]

E_{a}=k_{B}Tlog(2)(1 +T/ΔT)

Plugging in the numbers, one obtains

E_{a}= 0.53 eV

This is a typical activation energy in biological catalytic reactions. One can see that it is significantly larger than the diffusion activation energies we consider here.

For a proper investigation of diffusion barriers, one would take a
larger unitcell so that the coverage is not as high as it is in this
example. Probably, in this example the adsorption positions of the H
atoms are dominated by trying to reduce their repulsive dipole
interaction (they form a dipole when adsorbing to the Al(111)
surface). One would also have to scan the whole diffusion path. In
this example one could do this by keeping the *x-y* positions of the H
rigid and just relaxing the *z*-component. By changing the *x-y*
components by hand, one can scan several paths. For more complicated
problems, one would employ the Nudged Elastic Band method.

#### 1.3 Electrostatic potential

The electrostatic potential for H adsorbed in the ontop site looks as shown in this figure:

The H atom is positioned at *z* = 1.55 Å, and the dipole layer is at
*z* = 6.0 Å. The workfunction for the surface with a monolayer of
hydrogen is 5.7 eV and for the clean Al(111) surface we get 3.7 eV.
The difference between the two sides of the slab is due to the dipole
inducuced by H-adsorption: The H atom pulls electrons away from the
surface (μ = -0.078 *e*Å). The experimental value for the
workfunction of a clean Al surface is 4.25 eV - a bit far from our
calculated value of 3.7 eV, but that is to be expected for such a thin
slab. For a thicker slab we would get better agreement.

All the number can be found in Dacapos output file:

> grep DIP out.txt DIP: Z_0 [A] 6.000000 nzzero: 41 39 40 42 43 DIP: DIP: u_damp u Field_1 Field_2 Work_f1 Work_f2 V_jump E_dip DIP: [eA] [eA] [V/A] [V/A] [eV] [eV] [V] [eV] DIP: 0.011 0.011 - - - - 0.283 -0.128 DIP: 0.003 -0.142 -0.303 -0.370 - - 0.083 -0.037 DIP: -0.005 -0.180 -0.368 -0.430 - - -0.117 0.053 ... ... ... DIP: -0.077 -0.077 0.012 -0.016 5.663 3.695 -1.965 0.889 DIP: -0.077 -0.077 0.011 -0.016 5.662 3.696 -1.965 0.888 DIP: -0.077 -0.078 0.009 -0.019 5.659 3.693 -1.971 0.891 DIP: -0.077 -0.078 0.010 -0.018 5.663 3.692 -1.974 0.893