Your first band structure in VASP: a walkthrough for chemical engineers

When I started DFT for my undergraduate thesis, I hit a wall. Chemical engineering taught me transport phenomena, thermodynamics and reactor design. Quantum mechanics was one short course. VASP tutorials assumed I already knew what a Brillouin zone was, why the band gap came out wrong, and which of the four input files I had broken this time.

This is the guide I wanted then. It takes one simple material, silicon, from an empty folder to a band structure plot and a density of states. It explains each setting in plain language, and it ends with the changes you need for a real research material.

The quantum mechanics you actually need

You do not need to derive anything to run a good calculation, but you do need four ideas.

Electrons in a crystal behave as waves. In a periodic crystal, each electron state can be labelled by a wave vector $\mathbf{k}$, which says how the wave's phase changes from one unit cell to the next. All the distinct $\mathbf{k}$ values fit inside one cell of reciprocal space, called the first Brillouin zone.

A band structure is energy plotted against k. For each $\mathbf{k}$ there is a ladder of allowed energies. Follow one rung of the ladder as $\mathbf{k}$ moves along a path through the Brillouin zone and you trace a band. The gap between the highest filled band (valence band maximum, VBM) and the lowest empty one (conduction band minimum, CBM) is the band gap.

DFT replaces many electrons with their density. Solving for every electron at once is impossible. Density functional theory solves for the electron density $n(\mathbf{r})$ instead, with one approximate term, the exchange–correlation functional $E_{xc}[n]$, standing in for everything quantum-mechanically hard. PBE is the usual first choice. It is fast and gets structures right, but it underestimates band gaps, often by 30–50%. Remember this before you compare your number with an experiment.

The calculation is an iterative loop. Guess a density, compute the potential it creates, solve for new electron states, build a new density from them, and repeat until nothing changes. This self-consistent field (SCF) loop works like converging a recycle stream in Aspen or DWSIM: same idea, same failure modes.

Two numerical settings control accuracy, and both work like mesh resolution in CFD:

• ENCUT, the plane-wave cutoff energy, sets how finely the wavefunctions are resolved in space.
• The k-point mesh sets how finely the Brillouin zone is sampled.

Too coarse and the answer is wrong; too fine and you wait for nothing. You test both until the total energy stops changing (below about 1 meV per atom).

The four input files

Every VASP run reads four files from the folder it runs in.

For silicon, the primitive cell of the diamond structure has two atoms:

textSi (diamond, primitive FCC cell)
5.43
  0.0  0.5  0.5
  0.5  0.0  0.5
  0.5  0.5  0.0
Si
2
Direct
  0.00 0.00 0.00
  0.25 0.25 0.25

Line 2 is a scale factor (the lattice constant in Å), lines 3–5 are the lattice vectors in units of it, then the element names, how many atoms of each, and their positions as fractions of the lattice vectors.

The POTCAR comes from the PAW library that ships with your VASP licence. For several elements you concatenate their files in the same order as the element line in POSCAR. Getting that order wrong gives no error, just a nonsense result.

The workflow: relax, SCF, bands, DOS

A band structure is not one calculation but four, each in its own folder:

1. Relax the structure to its lowest-energy geometry.
2. SCF: converge the charge density on a uniform k-mesh and save it.
3. Bands: read that density and compute energies along a high-symmetry path.
4. DOS: read the same density and compute energies on a dense uniform mesh.

Steps 3 and 4 are non-self-consistent: they do not update the density, they only evaluate states in it. That is why they can use k-points that would be useless for converging the density.

Step 1: relax the structure

INCAR:

textSYSTEM = Si relax
PREC   = Accurate
ENCUT  = 400        ! eV; about 1.3 x the largest ENMAX in POTCAR when the cell can change
EDIFF  = 1E-6       ! electronic convergence, eV
ISMEAR = 0          ! Gaussian smearing: the safe choice for semiconductors
SIGMA  = 0.05
IBRION = 2          ! conjugate-gradient ionic relaxation
ISIF   = 3          ! relax atoms, cell shape and cell volume
NSW    = 50         ! maximum ionic steps
EDIFFG = -0.01      ! stop when all forces are below 0.01 eV/Å

KPOINTS, a Γ-centred 8 × 8 × 8 mesh:

textAutomatic mesh
0
Gamma
8 8 8
0 0 0

When it finishes, CONTCAR holds the relaxed structure. Copy it to POSCAR and run the relaxation once more. When the cell volume changes, the plane-wave basis set no longer fits it perfectly (an error called Pulay stress), and a second run from the new cell fixes that.

Step 2: the self-consistent run

Copy the relaxed POSCAR, the POTCAR and the KPOINTS into a new folder, then use:

textSYSTEM = Si scf
PREC   = Accurate
ENCUT  = 400
EDIFF  = 1E-6
ISMEAR = 0
SIGMA  = 0.05
IBRION = -1         ! no ionic moves
NSW    = 0
LCHARG = .TRUE.     ! write CHGCAR, the converged charge density
LWAVE  = .FALSE.

When it finishes, read off the Fermi energy. You need this value later, not the one from the band run:

bashgrep E-fermi OUTCAR

Step 3: the band structure

The bands are computed along straight lines between high-symmetry points of the Brillouin zone. For a face-centred cubic cell like silicon's, the standard path visits Γ, X, W, K and L. In a new folder, copy POSCAR, POTCAR and the SCF run's CHGCAR, then write a KPOINTS in line mode:

textSi band path (FCC): G-X-W-K-G-L
40              ! points along each segment
Line-mode
Reciprocal
0.000 0.000 0.000 ! \Gamma
0.500 0.000 0.500 ! X

0.500 0.000 0.500 ! X
0.500 0.250 0.750 ! W

0.500 0.250 0.750 ! W
0.375 0.375 0.750 ! K

0.375 0.375 0.750 ! K
0.000 0.000 0.000 ! \Gamma

0.000 0.000 0.000 ! \Gamma
0.500 0.500 0.500 ! L

You do not have to look these points up by hand. SeeK-path and VASPKIT (task 303 for bulk crystals) generate the path for any structure. One trap: SeeK-path may return a standardised cell. If it does, use that cell for every step, because the coordinates are only valid for it.

INCAR for the band run:

textSYSTEM = Si bands
PREC   = Accurate
ENCUT  = 400        ! must match the SCF run
ISMEAR = 0
SIGMA  = 0.05
ICHARG = 11         ! read CHGCAR and keep it fixed
NBANDS = 16         ! enough empty bands to see the conduction band clearly
LORBIT = 11         ! also write orbital projections
IBRION = -1
NSW    = 0
LCHARG = .FALSE.
LWAVE  = .FALSE.

Step 4: the density of states

The DOS counts how many states sit at each energy, so it needs a dense, uniform mesh, not a path. Copy the same POSCAR, POTCAR and CHGCAR into a fourth folder, use a 16 × 16 × 16 Γ-centred KPOINTS, and change the band INCAR to:

textICHARG = 11
ISMEAR = -5         ! tetrahedron method: sharp, accurate DOS
NEDOS  = 2001       ! energy grid points
LORBIT = 11         ! DOS projected on each atom and orbital (s, p, d)
NBANDS = 16

Plot and read it

pymatgen reads vasprun.xml and does the bookkeeping:

pythonfrom pymatgen.io.vasp import Vasprun
from pymatgen.electronic_structure.plotter import BSPlotter

efermi = Vasprun("scf/vasprun.xml").efermi          # from the SCF run, not the band run
bs = Vasprun("bands/vasprun.xml").get_band_structure(
    kpoints_filename="bands/KPOINTS", line_mode=True, efermi=efermi)

print(bs.get_band_gap())        # energy (eV), direct or indirect, and where the transition is
ax = BSPlotter(bs).get_plot(ylim=(-6, 6))
ax.figure.savefig("si_bands.png", dpi=300, bbox_inches="tight")

For silicon you should see the valence band maximum at Γ and the conduction band minimum partway along Γ–X: an indirect gap. PBE gives about 0.6 eV against the experimental 1.12 eV. That is the underestimate from the first section, not a broken calculation. If you need accurate gaps, the next steps are hybrid functionals such as HSE06, or the cheaper mBJ meta-GGA.

On the DOS, the projections tell you which atoms and orbitals make up each band edge. That is often the chemically interesting part: whether the gap opens between oxygen 2p and metal d states, for example.

Moving to a real material

Silicon is forgiving. A transition-metal oxide such as the NdCuO₃ perovskite in my thesis needs more care:

• Magnetism. Set ISPIN = 2 and give starting moments with MAGMOM. Antiferromagnetic orderings usually need a supercell, so that up and down spins sit on different atoms. Compare the total energies of the orderings to find the ground state.
• Localised d and f electrons. Plain PBE spreads them out too much. DFT+U (LDAU = .TRUE. with LDAUTYPE, LDAUL, LDAUU, LDAUJ) adds an on-site correction. Take U from the literature for similar compounds, and check how much your conclusions depend on it.
• Rare-earth 4f electrons. POTCARs such as Nd_3 freeze the 4f electrons into the core, which keeps the calculation stable and is often a sensible first approximation.
• Consistency. Every step must use the same ENCUT, ISPIN, MAGMOM and DFT+U settings. A band run with different settings from its SCF run is reading a density that does not belong to it.

Mistakes that cost me time

• POTCAR elements in a different order from POSCAR.
• Skipping the ENCUT and k-point convergence tests, then trusting the third decimal place.
• Using the band run's Fermi energy. Non-self-consistent runs do not determine it reliably.
• Setting ISMEAR = -5 in the band run. The tetrahedron method needs a uniform mesh and fails in line mode.
• Relaxing the cell with a low ENCUT and no second run.
• Reporting a PBE gap as if it were the experimental one.

If something here is unclear or wrong, or you get stuck on your own system, email me. I learnt most of this by asking, and I am happy to return the favour.