Skip to content

Calculating Band Structure in VASP

Jump to bottom
roeeasher edited this page Jun 18, 2026 · 4 revisions

How to Calculate a Band Structure in VASP Using the HSE06 Hybrid Functional

Calculating band structures with hybrid functionals like HSE06 in VASP is computationally demanding but provides significantly more accurate band gaps compared to standard GGA functionals. Because hybrid calculations evaluate exact exchange, the standard non-self-consistent approach used for PBE band structures requires a different workflow.

This tutorial walks you through the process, starting with the crucial first step: generating the converged ground-state charge density and wavefunctions.

Prerequisites and Best Practices

Before starting the Self-Consistent Field (SCF) calculation, ensure you have addressed the following:

  • Relaxed Structure: Your input structure (POSCAR) must be fully relaxed beforehand. Ideally, this relaxation should also be done using the HSE06 functional to ensure the lattice parameters and atomic positions are consistent with the electronic structure you are about to calculate.
  • Primitive Unit Cell: Always use the primitive unit cell rather than a conventional cell. Hybrid functional calculations scale poorly with the number of atoms and k-points. Reducing your system to the smallest possible primitive cell is essential to keep the computational cost manageable.

Step 1: The HSE06 SCF Calculation

The goal of this first step is to obtain a highly accurate, converged ground state. This run will generate the CHGCAR and WAVECAR files, which are absolutely required for the subsequent band structure calculation.

Example Input Files: Bulk Silicon

To make this concrete, let us look at the structural and k-point inputs for bulk silicon.

The POSCAR File

Here is the POSCAR for the primitive cell of silicon. Notice that it contains only 2 atoms. Using this minimal representation is crucial for hybrid calculations to keep the computational cost manageable.

Si2
1.0
  -0.0000000000000002    2.7218511864697263    2.7218511864697263
   2.7218511864697263    0.0000000000000000    2.7218511864697263
   2.7218511864697263    2.7218511864697263    0.0000000000000003
Si
2
direct
   0.2500000000000002    0.2500000000000001    0.2500000000000000 Si
   0.0000000000000000    0.0000000000000000    0.0000000000000000 Si

The KPOINTS File

For the SCF step, a uniform grid is required to accurately sample the Brillouin zone and converge the charge density. An 8x8x8 Gamma-centered mesh is an excellent choice for the primitive cell of bulk silicon.

Si
0
Gamma
8 8 8
0 0 0

The INCAR File

Below is the INCAR file optimized for an HSE06 static calculation.

# HSE06 Functional Parameters
LHFCALC = True
GGA = Pe
AEXX = 0.25
HFSCREEN = 0.2

# Electronic Convergence and Basis Set
ENCUT = 520
PREC = Accurate
PRECFOCK = Normal
LASPH = True
LREAL = False
EDIFF = 1e-06
NELM = 200
ALGO = All

# Smearing and Spin
ISPIN = 1
ISMEAR = 0
SIGMA = 0.03

# Ionic Relaxation (Static Calculation)
IBRION = -1
NSW = 0

# Output and Convergence Aids
ICHARG = 2
LCHARG = True
LWAVE = True
LELF = True
LAECHG = True
LMAXMIX = 4

# Potential Outputs
LVTOT = True
LVHAR = True

Breakdown of Important INCAR Parameters

Understanding the parameters in the INCAR file is crucial for a successful hybrid calculation. Here is a detailed breakdown of the key flags used in this SCF step:

1. Defining the HSE06 Functional

  • LHFCALC = True: Activates the Hartree-Fock (exact exchange) routine.
  • GGA = Pe: Specifies the PBE generalized gradient approximation for the correlation part.
  • AEXX = 0.25: Sets the fraction of exact exchange to 25%, which is standard for HSE06.
  • HFSCREEN = 0.2: Defines the range-separation (screening) parameter to 0.2 Å^-1.

2. Precision and Basis Set Settings

  • ENCUT = 520: The plane-wave cutoff energy. Always check your specific POTCAR recommendations, but 520 eV is generally a safe, high-quality default.
  • PREC = Accurate: Ensures a fine FFT grid is used to avoid aliasing errors, which is critical for accurate forces and stresses.
  • PRECFOCK = Normal: Controls the FFT grid for the exact exchange calculation. Normal provides an excellent balance of accuracy and computational cost.
  • LASPH = True: Includes non-spherical contributions within the PAW spheres. This is absolutely essential for accurate band gaps in hybrid calculations, especially for d-block and f-block elements.
  • LREAL = False: Forces projection in reciprocal space. For small primitive cells, this is faster and more accurate.

3. Electronic Convergence and Minimization

  • ALGO = All: Selects the electronic minimization algorithm. Standard Davidson (Normal) often struggles with exact exchange. All (or Damped) is highly recommended for hybrid functionals to ensure stable convergence.
  • EDIFF = 1e-06: Sets a tight electronic convergence criterion (1e-6 eV). Since this SCF run provides the foundation for the entire band structure, a strictly converged ground state is critical.
  • NELM = 200: Sets a high maximum number of electronic self-consistency steps, allowing the HSE06 calculation enough iterations to converge.

4. Smearing Settings

  • ISMEAR = 0: Uses Gaussian smearing, which is the correct and safest choice for semiconductors and insulators like bulk silicon.
  • SIGMA = 0.03: Sets a small smearing width to prevent unphysical entropy contributions to the total energy while aiding electronic convergence.

5. Crucial Outputs for the Band Structure Step

  • LWAVE = True & LCHARG = True: Writes the exact wavefunctions (WAVECAR) and charge density (CHGCAR). Both files are required to initialize the subsequent non-self-consistent band structure run.
  • LMAXMIX = 4: Writes and mixes the charge density up to l=4 (d-orbitals). If working with f-block elements, increase this to 6. This parameter significantly improves convergence in subsequent steps that read the CHGCAR.

The POTCAR File

POTCAR file could be obtained from the shared BMD directory, or wherever you get your POTCARs from.

Step 2: The Non-SCF Band Structure Calculation

With the ground state fully converged, the next step is to calculate the eigenvalues along the high-symmetry paths of the Brillouin zone. Because hybrid functionals require the full set of occupied wavefunctions from the uniform grid to calculate the exact exchange potential, we cannot simply replace the uniform grid with a line-mode grid like we do in standard PBE calculations.

Instead, we must keep the uniform grid and explicitly append the high-symmetry path k-points with a weight of exactly zero.

Setting Up the Calculation Directory

Create a new directory for this band structure calculation. You must copy the following files from your successful Step 1 (SCF) directory into this new folder:

  • CHGCAR: The converged charge density.
  • WAVECAR: The converged wavefunctions.
  • POSCAR: Your primitive unit cell.
  • POTCAR: The pseudopotentials.
  • IBZKPT: The irreducible Brillouin zone k-points generated during the SCF run (you will need this for the Python script below).

The INCAR File

Create a new INCAR file for this step. While it looks similar to the SCF INCAR, there are a few critical modifications required to read the previous state and extract the band projections.

# HSE06 Functional Parameters
LHFCALC = True
GGA = Pe
AEXX = 0.25
HFSCREEN = 0.2

# Electronic Convergence and Basis Set
ENCUT = 520
PREC = Accurate
PRECFOCK = Normal
LASPH = True
LREAL = False
EDIFF = 1e-06
NELM = 200
NELMIN = 5
ALGO = Normal

# Smearing and Spin
ISPIN = 1
ISMEAR = 0
SIGMA = 0.03

# Ionic Relaxation (Static Calculation)
IBRION = -1
NSW = 0

# Output and Read Flags
ISTART = 1
ICHARG = 1
LCHARG = True
LWAVE = True
LELF = True
LAECHG = True
LMAXMIX = 4
LORBIT = 11

# Potential Outputs
LVTOT = True
LVHAR = True

Key Differences: Non-SCF vs. SCF INCAR

When moving from the SCF ground state calculation to the non-SCF band structure step, several parameters must be changed to ensure VASP reads the pre-converged data and diagonalizes the Hamiltonian correctly.

  • ISTART = 1 and ICHARG = 1: These are the most critical changes.
    • ISTART = 1 instructs VASP to read the converged wavefunctions from the WAVECAR generated in Step 1.
    • ICHARG = 1 instructs VASP to read the charge density from the CHGCAR and extrapolate it. This locks in your converged ground state.
  • ALGO = Normal: The electronic minimization algorithm is switched from All (or Damped) back to Normal (the standard Davidson block Davidson algorithm). Because you are starting from highly converged wavefunctions, the standard and more stable Normal algorithm is perfectly sufficient and preferred for this exact exchange diagonalization.
  • NELMIN = 5: This parameter is added to force a minimum of 5 electronic steps. This guarantees that the wavefunctions for the newly added zero-weight path k-points are sufficiently converged within the fixed potential.
  • LORBIT = 11: Added to compute the spherical harmonic projections of the wavefunctions. This generates the PROCAR and DOSCAR files, which you will need if you intend to plot orbital-projected (fat) band structures.

Generating the Zero-Weight KPOINTS File

s mentioned, hybrid band structures require a combined KPOINTS file containing both the uniform grid (with their standard weights) and the path k-points (with zero weights). This prevents the newly added path points from altering the Fermi level and the exchange-correlation potential. You can automate this process using pymatgen. Below is a Python script (generate_kpoints.py) that reads your POSCAR and the IBZKPT file from the SCF run, generates the high-symmetry path, and outputs the correct KPOINTS format.

from pymatgen.core import Structure
from pymatgen.symmetry.bandstructure import HighSymmKpath
from pymatgen.io.vasp.inputs import Kpoints

# 1. Load the primitive structure and generate the path
struct = Structure.from_file("POSCAR")
kpath = HighSymmKpath(struct)

# 2. Get the explicit k-points along the path (40 points per segment)
path_kpts, path_labels = kpath.get_kpoints(line_density=40, coords_are_cartesian=False)

# 3. Read the uniform grid from the IBZKPT file copied from the SCF run
ibz = Kpoints.from_file("IBZKPT")
new_kpts = list(ibz.kpts)
new_weights = list(ibz.kpts_weights)

# 4. Append the path k-points with exactly 0.0 weight
for kpt in path_kpts:
    new_kpts.append(kpt)
    new_weights.append(0.0)

# 5. Write the combined list to the new KPOINTS file
hse_kpoints = Kpoints(comment="HSE Zero-Weight Path",
        num_kpts=len(new_kpts),
        style=Kpoints.supported_modes.Reciprocal,
        kpts=new_kpts,
        kpts_weights=new_weights)
hse_kpoints.write_file("KPOINTS")

Plotting the Band Structure

After your calculations finish, it is time to plot the band structure!

Here is a python script that plots the band structure and labels the high-symmetry points in the Brillouin zone:

import numpy as np
import matplotlib.pyplot as plt
from pymatgen.io.vasp.outputs import Vasprun
from pymatgen.electronic_structure.core import Spin
from pymatgen.core.structure import Structure
from pymatgen.symmetry.bandstructure import HighSymmKpath

# 1. Load the raw data
vr = Vasprun("vasprun.xml", parse_projected_eigen=True)
bs = vr.get_band_structure(line_mode=True)

kpts_cart = np.array([k.cart_coords for k in bs.kpoints])
energies_up = bs.bands[Spin.up]

is_spin_polarized = Spin.down in bs.bands
if is_spin_polarized:
    energies_down = bs.bands[Spin.down]

# Load POSCAR to get high symmetry points and their fractional coordinates
try:
    structure = Structure.from_file("POSCAR")
except FileNotFoundError:
    print("POSCAR not found. Falling back to the final structure from vasprun.xml.")
    structure = vr.final_structure

# Generate standard high symmetry paths based on the structure's spacegroup
kpath_generator = HighSymmKpath(structure)
symm_kpoints = kpath_generator.kpath['kpoints']

def find_kpoint_label(frac_coord, symm_dict, tol=1e-3):
    """Finds the closest high symmetry point label for a given fractional coordinate."""
    for label, coord in symm_dict.items():
        diff = frac_coord - coord
        # Wrap coordinates to account for periodic boundary conditions in reciprocal space
        diff -= np.round(diff)
        if np.linalg.norm(diff) < tol:
            # Format label for matplotlib LaTeX rendering
            if label.upper() == 'GAMMA' or label == r'\Gamma':
                return r"$\Gamma$"
            elif '_' in label:
                # Formats subscripts properly, e.g., B_1 to $B_{1}$
                parts = label.split('_')
                return f"${parts[0]}_{{{parts[1]}}}$"
            else:
                return label
    return "?"

# 2. Dynamically scan the path for nodes based on distance anomalies
x_dist = np.zeros(len(kpts_cart))
segment_starts = [0]

# tick_info stores tuples of (x_position, left_kpt_index, right_kpt_index)
# This allows us to handle jump nodes that have two different symmetry labels
tick_info = [(0.0, 0, 0)] 

for i in range(1, len(kpts_cart)):
    dist = np.linalg.norm(kpts_cart[i] - kpts_cart[i-1])
    
    if dist < 1e-4:
        # Distance is 0: Duplicate point (Continuous Node)
        x_dist[i] = x_dist[i-1]
        segment_starts.append(i)
        tick_info.append((x_dist[i], i-1, i))
        
    elif dist > 0.15:
        # Distance is large: Brillouin Zone Jump (Jump Node)
        x_dist[i] = x_dist[i-1]
        segment_starts.append(i)
        tick_info.append((x_dist[i], i-1, i))
        
    else:
        # Normal step along the path
        x_dist[i] = x_dist[i-1] + dist

# The end of the path is always the final node
segment_starts.append(len(kpts_cart))
tick_info.append((x_dist[-1], len(kpts_cart) - 1, len(kpts_cart) - 1))

print(f"Total zero-weight points processed: {len(kpts_cart)}")
print(f"Nodes dynamically detected: {len(tick_info)}")

# 3. Apply the LaTeX labels dynamically based on POSCAR symmetry
custom_labels = []
for pos, idx1, idx2 in tick_info:
    label1 = find_kpoint_label(bs.kpoints[idx1].frac_coords, symm_kpoints)
    label2 = find_kpoint_label(bs.kpoints[idx2].frac_coords, symm_kpoints)
    
    if idx1 == idx2 or label1 == label2:
        custom_labels.append(label1)
    else:
        # Construct compound labels for Brillouin zone jumps (e.g., X | Q)
        custom_labels.append(f"{label1} $\\mid$ {label2}")

tick_positions = [info[0] for info in tick_info]

# 4. Render the plot
fig, ax = plt.subplots(figsize=(12, 8))
num_bands = energies_up.shape[0]

for b in range(num_bands):
    for s in range(len(segment_starts)-1):
        start = segment_starts[s]
        end = segment_starts[s+1]
        
        ax.plot(x_dist[start:end], energies_up[b, start:end] - bs.efermi, 
                color='#1f77b4', linestyle='-', linewidth=1.5)
        
        if is_spin_polarized:
            ax.plot(x_dist[start:end], energies_down[b, start:end] - bs.efermi, 
                    color='#1f77b4', linestyle='--', linewidth=1.5)

# Format axes
for pos in tick_positions:
    ax.axvline(x=pos, color='black', linewidth=0.5)

ax.set_xticks(tick_positions)
ax.set_xticklabels(custom_labels, fontsize=16)
ax.set_ylabel(r"$E - E_f$ (eV)", fontsize=20)
ax.set_xlim(0, max(x_dist))
ax.set_ylim(-6, 6) 

# Legend
from matplotlib.lines import Line2D
custom_lines = [Line2D([0], [0], color='#1f77b4', linestyle='-', linewidth=1.5)]
legend_labels = ['Spin Up']
if is_spin_polarized:
    custom_lines.append(Line2D([0], [0], color='#1f77b4', linestyle='--', linewidth=1.5))
    legend_labels.append('Spin Down')
ax.legend(custom_lines, legend_labels, loc='lower left', fontsize=12)

fig.savefig("HSE06_BandStructure.png", dpi=300, bbox_inches="tight")
print("Success! Plot generated using POSCAR derived high symmetry labels.")

Clone this wiki locally