Work Function

This tutorial explains how to calculate and analyze slab work functions with VQkit for both VASP and Quantum ESPRESSO.

For a surface or slab calculation, VQkit uses:

\[\Phi = V_{\mathrm{vacuum}} - E_F\]

where:

\(V_{\mathrm{vacuum}}\) is the vacuum electrostatic potential;
\(E_F\) is the Fermi energy.

For an asymmetric slab, the two sides of the slab can have different vacuum levels and therefore different work functions.

VQkit reports:

lower-side vacuum level;
upper-side vacuum level;
lower-side work function;
upper-side work function;
mean work function;
vacuum-window flatness.

Open the work-function workflow

From the VQkit main menu choose:

[2] Post-processing

Then choose:

[3] Electronic structure

and select:

Work function

VQkit detects whether the current directory contains VASP or Quantum ESPRESSO work-function data.

VASP workflow

Required files

For VASP, VQkit requires:

LOCPOT

and a Fermi-energy source from either:

OUTCAR

or:

DOSCAR

The preferred Fermi-energy source is OUTCAR.

If it is unavailable, VQkit attempts to read the Fermi energy from DOSCAR.

Generating LOCPOT

A VASP slab calculation must write the electrostatic potential.

VQkit reports an error if LOCPOT is absent and reminds you to run an appropriate VASP static calculation with a suitable potential-output option, such as:

LVTOT = .TRUE.

or:

LVHAR = .TRUE.

depending on the potential definition required for your workflow.

The choice between total local potential and ionic + Hartree potential should be scientifically consistent with the work-function convention you intend to use.

Choose the slab-normal axis

When VQkit detects a VASP LOCPOT, it asks:

Slab-normal axis [z]:

For a conventional slab with vacuum along z, press Enter.

Allowed axes are:

x
y
z

Choose the axis perpendicular to the slab surface and vacuum region.

Planar averaging of LOCPOT

VQkit reads the three-dimensional LOCPOT potential grid and computes a planar average along the selected slab-normal axis.

For a z-normal slab, the planar average is conceptually:

\[\overline{V}(z) = \frac{1}{A} \int_A V(x,y,z)\,dx\,dy\]

The resulting one-dimensional potential profile is used to estimate the vacuum level on both sides of the slab.

Vacuum regions

VQkit uses edge windows on both ends of the planar potential profile.

By default, approximately 10% of the profile on each side is used.

The lower-side vacuum level is calculated as the mean potential in the first edge window.

The upper-side vacuum level is calculated as the mean potential in the final edge window.

VQkit also calculates the standard deviation of the potential inside each vacuum window.

A small standard deviation indicates a flatter and more reliable vacuum plateau.

Work-function calculation

Lower side

\[\Phi_{\mathrm{lower}} = V_{\mathrm{vac,lower}} - E_F\]

Upper side

\[\Phi_{\mathrm{upper}} = V_{\mathrm{vac,upper}} - E_F\]

Mean of both sides

VQkit also reports:

\[\Phi_{\mathrm{mean}} = \frac{ \Phi_{\mathrm{lower}} + \Phi_{\mathrm{upper}} }{2}\]

For a symmetric slab, the two side values should normally be similar.

For an asymmetric slab, they may differ significantly.

VASP output report

A typical VQkit report contains:

Code
Fermi energy
Vacuum level, lower
Vacuum level, upper
Work function, lower
Work function, upper
Mean of two sides
Vacuum-window std (L/U)
Source
Data
Plot

Example:

Fermi energy            : 5.120000 eV
Vacuum level, lower     : 9.850000 eV
Vacuum level, upper     : 10.120000 eV
Work function, lower    : 4.730000 eV
Work function, upper    : 5.000000 eV
Mean of two sides       : 4.865000 eV

Generated files

VQkit writes:

VQkit_WORK_FUNCTION.dat
VQkit_WORK_FUNCTION.png
VQkit_WORK_FUNCTION.pdf

The data file contains the planar potential profile and metadata including:

code;
axis;
Fermi energy;
lower vacuum level;
upper vacuum level.

The numerical table contains:

coordinate
planar_potential_eV

Work-function plot

The generated plot shows:

planar-averaged potential;
Fermi level;
lower-side vacuum level;
upper-side vacuum level when different.

The x-axis is the coordinate along the chosen slab-normal direction.

The y-axis is potential energy in eV.

Vacuum-flatness warning

VQkit evaluates the standard deviation of each vacuum window.

If either vacuum-window standard deviation exceeds approximately:

0.10 eV

VQkit displays a warning.

A non-flat vacuum region may indicate problems such as:

insufficient vacuum thickness;
a strong slab dipole;
inappropriate slab-normal axis;
incomplete electrostatic convergence;
need for dipole correction;
a vacuum window that still overlaps the slab potential.

Example: VASP slab

Assume your completed slab directory contains:

LOCPOT
OUTCAR
POSCAR
INCAR
KPOINTS

and vacuum is along z.

Step 1 — Start VQkit

Run VQkit in the slab directory.

Step 2 — Open work function

Navigate to:

Post-processing
Electronic structure
Work function

Step 3 — Confirm the slab axis

At:

Slab-normal axis [z]:

press Enter.

Step 4 — Inspect the report

Check:

Fermi energy;
lower vacuum level;
upper vacuum level;
lower work function;
upper work function;
vacuum-window standard deviations.

Step 5 — Inspect the potential profile

Open:

VQkit_WORK_FUNCTION.png

Verify that the potential reaches a clear plateau in the vacuum region.

Step 6 — Check asymmetry

If the lower and upper work functions differ substantially, determine whether the slab is physically asymmetric.

Possible reasons include:

different terminations;
adsorption on only one surface;
an interfacial dipole;
an applied field;
incomplete dipole correction.

Quantum ESPRESSO workflow

For QE, VQkit uses a planar-averaged electrostatic potential profile.

A preferred pre-existing profile filename is:

VQkit_QE_PLANAR_POTENTIAL.dat

VQkit also recognizes filenames such as:

average.dat
avg.dat
planar_average.dat
planar_potential.dat
potential.avg

If a supported planar potential already exists, VQkit analyzes it directly.

QE Fermi energy

VQkit searches local pw.x outputs for the Fermi energy.

It can recognize common QE output lines such as:

the Fermi energy is ...

and also occupied-level information where appropriate.

Preferred output filenames include:

scf.out
pw.out
nscf.out

VQkit also scans other .out and .pwo files in the directory.

No QE planar profile found

If no planar-average potential exists, VQkit prepares the QE post-processing workflow automatically.

It creates:

VQkit_QE_WORKFUNCTION_PP.in
VQkit_QE_WORKFUNCTION_AVERAGE.in
VQkit_QE_WORKFUNCTION_STEPS.txt

The workflow uses:

pp.x to generate electrostatic-potential data;
average.x to calculate the one-dimensional planar average.

QE pp.x input

The generated pp.x input has the form:

&INPUTPP
  prefix = '...'
  outdir = '...'
  filplot = 'VQkit_QE_VBARE_VH.dat'
  plot_num = 11
/

VQkit attempts to determine prefix and outdir from the completed QE calculation.

QE saved data

QE post-processing requires access to:

<prefix>.save

VQkit checks several possible locations for the saved-data directory, including locations such as:

current directory
tmp/
out/
scratch/
nearby tmp/
nearby scf/tmp/

If the calculation was copied from an HPC system but the <prefix>.save directory was not copied, automatic local post-processing cannot run.

In that case, either:

run pp.x and average.x on the HPC system; or
copy the complete <prefix>.save directory to the local machine.

QE average.x input

VQkit also creates:

VQkit_QE_WORKFUNCTION_AVERAGE.in

The default planar-averaging direction is:

idir = 3

which corresponds to the third lattice direction.

The generated workflow uses a default of approximately:

npt = 500
awin = 1.0

These values should be reviewed when needed for the specific slab and QE setup.

Automatic QE execution

If all of the following are available:

accessible <prefix>.save data;
pp.x in PATH;
average.x in PATH,

VQkit attempts to run the post-processing automatically.

It first runs the equivalent of:

pp.x < VQkit_QE_WORKFUNCTION_PP.in > VQkit_QE_WORKFUNCTION_PP.out

Then:

average.x < VQkit_QE_WORKFUNCTION_AVERAGE.in > VQkit_QE_WORKFUNCTION_AVERAGE.out

If average.x successfully creates:

avg.dat

VQkit copies it to:

VQkit_QE_PLANAR_POTENTIAL.dat

and immediately proceeds with work-function analysis.

Manual QE execution

If automatic execution is not possible, follow the instructions in:

VQkit_QE_WORKFUNCTION_STEPS.txt

The manual sequence is:

Step 1

Run:

pp.x < VQkit_QE_WORKFUNCTION_PP.in > VQkit_QE_WORKFUNCTION_PP.out

Step 2

Run:

average.x < VQkit_QE_WORKFUNCTION_AVERAGE.in > VQkit_QE_WORKFUNCTION_AVERAGE.out

Step 3

After average.x produces:

avg.dat

rename or copy it to:

VQkit_QE_PLANAR_POTENTIAL.dat

Step 4

Run the VQkit Work Function option again.

QE unit handling

VQkit reads the numerical planar-average file and determines whether the potential is already expressed in eV.

If the potential is identified as Rydberg, VQkit converts it using:

\[1\ \mathrm{Ry} = 13.605693122994\ \mathrm{eV}\]

The final work function is reported in eV.

QE work-function calculation

Once the profile is available, the calculation is the same basic relation:

\[\Phi = V_{\mathrm{vacuum}} - E_F\]

VQkit evaluates both edge vacuum windows and reports:

lower vacuum level;
upper vacuum level;
lower work function;
upper work function;
mean work function.

Example: QE slab

Assume a completed QE SCF directory contains:

scf.in
scf.out
tmp/

and:

tmp/<prefix>.save/

is accessible.

Step 1 — Start VQkit

Start VQkit in the SCF directory.

Step 2 — Open Work Function

Navigate to:

Post-processing
Electronic structure
Work function

Step 3 — Let VQkit search for a profile

If no planar potential exists, VQkit generates the pp.x and average.x inputs.

Step 4 — Automatic route

If pp.x, average.x and the saved QE data are available, VQkit attempts to generate the profile automatically.

Step 5 — Manual route

If automatic execution is unavailable, use:

VQkit_QE_WORKFUNCTION_STEPS.txt

to run the two QE post-processing commands manually.

Step 6 — Analyze the resulting profile

Once:

VQkit_QE_PLANAR_POTENTIAL.dat

exists, run Work Function again.

Step 7 — Inspect the result

Check:

VQkit_WORK_FUNCTION.dat
VQkit_WORK_FUNCTION.png
VQkit_WORK_FUNCTION.pdf

and verify the vacuum plateau.

Symmetric versus asymmetric slabs

Symmetric slab

For an ideal symmetric slab:

\[\Phi_{\mathrm{lower}} \approx \Phi_{\mathrm{upper}}\]

A large difference may indicate:

numerical noise;
inadequate vacuum;
an unintended structural asymmetry;
an electrostatic dipole.

Asymmetric slab

For an asymmetric slab, different work functions on the two sides can be physically meaningful.

Examples include:

different surface terminations;
one-sided adsorption;
supported films;
heterointerfaces;
polar slabs.

In such cases, report the side-specific work function rather than relying only on the simple average.

Common problems

VASP LOCPOT missing

If VQkit reports that LOCPOT is missing, rerun an appropriate static calculation that writes the electrostatic potential.

Fermi energy not found

For VASP, keep:

OUTCAR

or:

DOSCAR

For QE, keep a pw.x output containing the Fermi energy.

Wrong slab-normal axis

If the selected axis is not perpendicular to the slab, the planar potential will not show the correct vacuum plateau.

Choose x, y or z according to the vacuum direction.

Vacuum region is not flat

If the vacuum-window standard deviation is large, inspect:

vacuum thickness;
dipole correction;
slab asymmetry;
electrostatic convergence;
selected averaging direction.

QE saved data unavailable

The text output alone is not sufficient for pp.x.

The complete:

<prefix>.save

directory is required.

QE pp.x not found

If VQkit reports:

pp.x was not found on PATH

load the appropriate QE module or add the QE binary directory to PATH.

QE average.x not found

If average.x is not available locally, run it on the system where your QE installation provides it.

Scientific cautions

Vacuum thickness

The vacuum region should be sufficiently large for the planar potential to reach a stable plateau.

A small vacuum region can produce an unreliable vacuum level.

Dipole correction

Asymmetric slabs can generate an artificial periodic dipole.

Use the appropriate dipole-correction strategy for the electronic-structure code and slab geometry when required.

Two-sided reporting

For asymmetric surfaces, report the lower and upper work functions separately when they correspond to physically distinct surfaces.

Energy reference

The Fermi energy and vacuum potential must come from mutually consistent electronic-structure data.

Potential choice

Different potential definitions can lead to different numerical profiles.

State clearly which VASP or QE potential definition and post-processing procedure was used.

Convergence

The calculated work function should be checked for convergence with respect to:

slab thickness;
vacuum thickness;
k-point mesh;
plane-wave cutoff;
electronic convergence;
dipole treatment.

VQkit automates the numerical extraction and plotting, but the physical quality of the work function still depends on the underlying slab calculation.