12567
Comment:

← Revision 6 as of 20101020 09:11:46 ⇥
12567
converted to 1.6 markup

No differences found! 
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 xyplane 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 zdirection.
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_{N1} by substituting N by N  1 in the equation above. By calculating NE_{N1} and (N  1)E_{N} and subtracting them one can show
σ = (NE_{N1}  (N  1)E_{N}) / (2A)
Using a kpoint 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 ycoordinates can be calculated using the scaled positions and the unit vectors. For the zcoordinate 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 AlH bonds to be the sum of the covalent radii d_{AlH} = r_{Al} + r_{H}. Then the zcoordinate 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 kpoint 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 = A exp[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}T log(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 xy positions of the H rigid and just relaxing the zcomponent. By changing the xy 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 Hadsorption: 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