Workfunction of hcp (0001) surfaces

In this notebook, we will show how to calculate the workfunction of selected hcp(0001) surfaces using VASP. Please keep in mind that the parameters used here give no converged results. They have been chosen to demonstrate the workflow using inexpensive calculations. For converged results, parameters such as lattice parameters, plane-wave energy cutoffs, reciprocal space sampling or the need to perform spin polarized calculations have to be carefully chosen

import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
import pandas as pd
import time
from pyiron import Project
pr = Project("hcp_workfunction")

Calculating the Workfunction of Mg(0001)

Structure creation

We use the create_surface() function which uses the ASE surface generator to build our surface slab structure

 # Now we set-up the Mg (0001) surface
a = 3.1919
c = 5.1852

# Vacuum region to break the periodicity along the z-axis
vac = 10
size = (2, 2, 4)
Mg_0001 = pr.create_surface("Mg",

Using selective dynamics

We use selective dynamics to restrict relaxation to the surface atoms (first and last Mg layers). We use the advanced array indexing options available in the NumPy package (see here) to detect which atoms are at the surface and then freeze the rest

# Initially freeze all the atoms
Mg_0001.add_tag(selective_dynamics=[False, False, False])

# Find which atoms are at the surface
# (based on the z-coordinate)
pos_z = Mg_0001.positions[:, 2]
z_min, z_max = np.min(pos_z), np.max(pos_z)
eps = 1e-4
relax_indices = np.argwhere(((pos_z - eps) > z_min)
                            & ((pos_z + eps) < z_max ))
relax_indices  = relax_indices.flatten()

# Now allow these atoms to relax

Mg_0001.selective_dynamics[relax_indices] = [True, True, True]

Setup and execution

To automate the calculation we define a function that has as input the project object, structure, job_name, Fermi smearing width, the type of k-point sampling and the plane-wave energy cutoff

def get_ham(proj, basis, name, sigma=0.1, mesh="GP", encut=350):
    ham = proj.create_job(pr.job_type.Vasp, name)
    # Setting fermi-smearing
    ham.set_occupancy_smearing(smearing="fermi", width=sigma)
    # Ionic_minimization
    ham.structure = basis
    if mesh == "GP":
        # Only the Gamma point
    elif len(mesh) == 3:
    return ham
ham_vasp = get_ham(proj=pr,

Submitting to the queue (optional)

If you use a cluster installation of pyiron, you can send the created jobs to the cluster by specifying the name of the queue and the number of cores

# queue = ham_vasp.server.list_queues()[-1]
# ham_vasp.server.queue = queue
# ham_vasp.server.cores = 20

Choosing an appropriate executable


Since this example uses the \(\Gamma\) point only, we can use the VASP Gamma-only version. If you use more k-points choose an appropriate executable

ham_vasp.executable.version = "5.4_gamma"


The job is ready for execution


Post processing

To analyze the results we ensure that the job is finished (the if statement in the first line). We then compute the work function by subtracting the Fermi-level from the vacuum level

\(\Phi = V_{vac} - \epsilon_F\)

if ham_vasp.status.finished:
    # Get the electrostatic potential
    epot = ham_vasp.get_electrostatic_potential()

    # Compute the lateral average along the z-axis (ind=2)
    epot_z = epot.get_average_along_axis(ind=2)

    # Get the final relaxed structure from the simulation
    struct = ham_vasp.get_structure(iteration_step=-1)
    r = np.linalg.norm(struct.cell[2])
    z = np.linspace(0, r, len(epot_z))

    # Computing the vacuum-level
    vac_level = np.max(epot_z)

    # Get the electronic structure
    es = ham_vasp.get_electronic_structure()
    print("wf:", vac_level - es.efermi)
    plt.plot(z, epot_z - vac_level)
    plt.xlim(0, r)
    plt.axhline(es.efermi - vac_level,
    plt.xlabel("z ($\AA$)")
    plt.ylabel("V - V$_{vac}$");
wf: 3.37343565133

Looping over a series of hcp(0001) surfaces

We now repeat the workflow for a set of hcp metals (the chosen lattice parameters are approximate). Note that if you use the same naming convention, pyiron detects that a job with the same name exists (“Mg_0001”) and loads the output from this calculation rather than launch a new job with the same name.

hcp_dict = {"Zn": {"a":2.6649, "c": 4.9468},
            "Mg": {"a": 3.1919, "c": 5.1852},
            "Co": {"a": 2.5071 , "c": 4.0695},
            "Ru": {"a": 2.7059 , "c": 4.2815}}
vac = 10
size = (2, 2, 4)
for element, lattice_parameters in hcp_dict.items():
    surf = pr.create_surface(element,
                             orthogonal=True, vacuum=vac)
    surf.add_tag(selective_dynamics=[False, False, False])
    pos_z = surf.positions[:, 2]
    z_min, z_max = np.min(pos_z), np.max(pos_z)
    eps = 1e-4
    relax_indices = np.argwhere(((pos_z - eps) > z_min)
                                & ((pos_z + eps) < z_max ))
    relax_indices  = relax_indices.flatten()
    surf.selective_dynamics[relax_indices] = [True, True, True]
    job_name = "{}_0001".format(element)
    ham = get_ham(pr, surf,
    #ham.server.cores = 20
    #ham.server.queue = queue
    ham.executable.version = '5.4_gamma'

Loading and analyzing

Now we iterate over all jobs in this project and calculate the workfunction. We also time how long the cell takes to execute

t1 = time.time()
for ham in pr.iter_jobs():
    if ham.status.finished:
            final_struct = ham.get_structure(iteration_step=-1)
            elec_structure = ham.get_electronic_structure()
            e_Fermi = elec_structure.efermi
            epot = ham.get_electrostatic_potential()
            epot_z = epot.get_average_along_axis(ind=2)
            vacuum_level = np.max(epot_z)
            wf = vacuum_level - e_Fermi
            element = final_struct.get_majority_species()[-1]
            hcp_dict[element]["work_func"] = wf
t2 = time.time()
print("time: {}s".format(t2-t1))
time: 9.250723838806152s

Compiling data in a table using pandas

df = pd.DataFrame(hcp_dict).T
df = df.rename(columns={'a': 'a [A]',
                        'c': 'c [A]',
                        'work_func': 'wf [eV]'})
    a [A]  c [A]  wf [eV]
Co  2.507  4.069    5.569
Mg  3.192  5.185    3.373
Ru  2.706  4.282    5.305
Zn  2.665  4.947    3.603
[ ]: