This tutorial will demonstrate how to compute the work function of a metal; our chosen example is silver.

The basic idea is simple: we create a slab model of the Ag [100] surface, and by setting the appropriate boundary conditions in the multi-grid Poisson solver in ATK, to normalize the effective potential to zero far away from the surface, the work function becomes equal to the Fermi level.

It should be pointed out, already from the beginning, that the precise value of the work function is rather sensitive to several numerical parameters. The proper approach for obtaining a converged value is to carefully study the influence of each of these parameters. This tutorial will show how to compute the work function itself, for a specific parameter choice.

It will be assumed that you have some experience with ATK, on the level of having computed some of the basic tutorials and performed e.g. a band structure calculation or similar. Thus we will not go into every step in full detail.

Setting up the geometry

1. Start by picking up silver from the Database . Send the structure to the Builder .
2. In the Builder, cleave the crystal along the (default) [100] direction, by choosing Transform>Cleave from the menu.
   
3. Repeat the resulting structure 7 times along Z. (This is one of the parameters that may need to be changed to ensure convergence.)
   
4. We need to add vacuum on both sides of the slab. We will leave 5 Å on the left side, so that no basis functions extend to the edge of the cell (the basis set radius for silver with default parameters is 7 Bohr), and 20 Å on the right hand side, to ensure that the effective has enough distance to decay. (If this is enough or not is another convergence parameter to study.)
   
   Therefore, extend the unit cell by 25 Å by changing the Z-component of the C vector to 53.5999 Å.
5. Then, select all atoms (Ctrl+A) and translate them by 5 Å along the Z axis.
   
6. To finalize the geometry setup, select the atom with the largest Z coordinate (31.5571 Å) and turn it into a ghost atom (menu Transform>Ghost).
   A ghost atom in ATK is an "atom" with a basis set but without any pseudopotential core or charge. The basis set orbitals that can be populated in order to host a finite electron density in a region where there are no real atoms. The alternative to having a ghost atom on the surface is to extend the range of the basis functions. This setting applies to all atoms, however, and makes the calculation take longer, especially if the system is not as simple as in this example.
7. In result, we will have 13 layers of silver in the slab, while the ghost atom on the surface extends the range of the electron density into the vacuum so that it has time to decay.
8. Send the ready structure to the Script Generator .

Defining the parameters of the calculation

1. In the Script Generator, insert a New Calculator by double-clicking it in the "Blocks" panel, as well as an Analysis>MolecularEnergySpectrum, which will help us compute the Fermi level. Also set a suitable file name for the NetCDF file (with full path, if you intend to run the script via the Job Manager).
   
2. Next, double-click the New Calculator in the Script panel and set the following parameters:
   • For a good description of the electron density, let us use 9x9 k-points in the surface plane. (Whether this is a high enough value - or too high, perhaps - for this or any other material, requires a systematic investigation.)
     
   • For silver, the SingleZetaPolarized basis set is most likely sufficient, so we can save some calculation time this way. (Other materials may not be described well enough by this basis set; this has to be carefully checked in each case.)
     
   • In general it is worthwhile considering performing the calculation also for GGA, to check the difference, but we will not do that here.
   • Finally - the most important setting: select the Multi-grid Poisson solver, and set the Neumann boundary condition on the left C face, and Dirichlet on the right. To achieve this, you must click the "chain" icon (marked in green in the picture) to unlink the settings for the two respective surfaces.
     The Dirichlet boundary condition means that we normalize the effective potential to zero on the right cell boundary in C, far away from the surface. We may not choose the same value (or any other constant value) on the left boundary as well, and obviously we also cannot use periodic boundary conditions. We can, however demand the gradient of the potential to vanish on the left of the system (Neumann condition).

Computing the work function

There is one more detail to take care of, and for that you need to edit the script manually.

Therefore, send the script from the Script Generator to the Editor .

Add the following line to the end of the script:

nlprint("Work function = %s eV" % numpy.abs(molecular_energy_spectrum.fermiLevel().inUnitsOf(eV)))

Since we have normalized the effective potential to be zero far away from the surface, the work function equals the Fermi level in the system, and this line extracts it from the molecular energy spectrum and prints it out to the log file. We use the ATK function nlprint which only prints on the master node, if you run in parallel, while a standard Python print statement would be executed on all nodes.

Save the script!

Then, run it by sending it to the Job Manager , or by transferring it to a computational server. On a typical machine it should take less than 5 minutes to finish, even in serial.

Results

The computed work function, printed by the script at the end of the log file, is 4.45 eV, which is indeed close to the experimental value, which ranges from 4.22 to 4.81 eV [1,2], depending on the method used to measure it.

Just to demonstrate that we can also get accurate results for other metals, see the comparisons with experiments in the table below (all values in eV):

* Experimental values are from Ref. 3, except where otherwise indicated.
** Fortuitously, the calculated work function for gold with SingleZetaPolarized is 5.28 eV.
‡ The 110 surface is known to undergo reconstruction; the computed value shown here does not take that into account.
† All calculations were performed using DoubleZetaPolarized basis set and 9x9 k-points, and other parameters as in the example above, for (100) surfacese except where otherwise indicated. Note that this mean that all values are not necessarily converged in these parameters. The surfaces are also not relaxed.

Comments

• Above we have commented in a few places on the importance of carefully checking convergence in certain numerical as well as geometrical parameters. The work function is quite sensitive to all of them, and so, for convenience, let us summarize the most important ones:
  • Number of slab layers
    (for silver, tests show that 13 layers, as we used above, is sufficient, but any smaller number gives different results)
  • Extension of vacuum, on both sides
    (20 Å on the right surface, as used in this example, is on the safe side; at least for silver [100], even just 10 Å is enough)
  • K-point sampling
    (9x9 gives reasonable results, but as always with k-point sampling the convergence in this parameter is not a simple monotonous function)
  • Basis set
    (we need the polarization orbitals to get p-shaped functions for the transition metals, and hence we choose SingleZetaPolarized for silver, but in other parts of the periodic table other rules may apply, which make e.g. the DoubleZeta basis set a good work horse; DoubleZetaPolarized is always a safe choice, although it demands a bit more computation time)
  • Exchange-correlation functional
    ("simple" metals are generally well described by LDA, but this is certainly not true for all materials)
• The script can be run in parallel (MPI) for some performance benefit, due to the large number of k-points which is one of the parameters ATK parallelizes over.
• The work function depends on which face you are looking at.
• For silver, the population of the ghost atom is rather small, but not zero (0.078 electrons). In other cases, like Ba, the ghost atom population is almost 0.3 electrons, which are taken from the outermost "real" Ba atom, indicating that indeed the kinetic energy can be lowered substantially if we let the electron cloud protrude from the surface a bit more.

References

[1] A. W. Dweydari and C. H. B. Mee, physica status solidi (a) 27, 223 (1975)
[2] M. Chelvayohan and C. H. B. Mee, J. Phys. C: Solid State Phys. 15, 2305 (1982)
[3] CRC Handbook on Chemistry and Physics, pages  12-114 (2008)
[4] J. K. Grepstad, P. O. Gartland and B. J. Slagsvold, Surface Science 57, 348-362 (1976)
[5] J. C. Rivière, Appl. Phys. Lett. 8, 172 (1966)