Dear Shweta,

To evaluate the band edge positions vs NHE from VASP band structure calculations for both bulk/slab structures follow these steps:

*Using VASP

Perform a self-consistent field (SCF) calculation using INCAR, KPOINTS, POSCAR, and POTCAr files.

Run a non-self-consistent field calculation to obtain the band structure.

**Determine the Band Edges

From the EIGENVAL/PROCAR file, extract the valence band maximum and conduction band minimum.

***Align the Energy Levels to the Vacuum Level

For bulk calcs: The absolute position of the VBM and CBM needs to be referenced relative to an external standard such as the vacuum level.

For slab ones: The vacuum potential can be obtained from the planar-averaged electrostatic potential (LOCPOT).

Important is to:

£ Extract the macroscopic average of the electrostatic potential across the slab.

££Define the vacuum level as the highest potential far from the slab.

To cmpute the energy alignment:

E_(VBM/CBM}^{absolute} = E_{VBM/CBM}^{VASP} - E_{Fermi} + E_{vacuum}

****Convert to the normal Hydrogen electrode scale

£The vacuum level is referenced as 0 eV, and the NHE reference level is -4.44 eV relative to vacuum.

£Convert the band edge positions to NHE via

E_{VBM/CBM}^{NHE} = E_{VBM/CBM)}^{abs} + 4.44 eV

***** If necessary, apply GGA+U or hybrid functionals (HSE06) to improve accuracy.

PS: Now, If you want to obtain similar plots:

!!! Extract Band Edges from VASP using the DOSCAR or EIGENVAL file to determine the CBM and VBM energies.

!!! Align them to the vacuum level using the macroscopic electrostatic potential from LOCpot in slab calculations.

!!! Convert to NHE Scale using following formulae

E_{NHE] = E_{vacuum} + 4.44 eV

Meanwhile, check please following script on Python for instance:

import numpy as np

import matplotlib.pyplot as plt

from pymatgen.io.vasp import Vasprun, Locpot

# Constants

NHE_SHIFT = 4.44

FERMI_REFERENCE = 0

# Load data from VASP output files

vasprun = Vasprun("vasprun.xml") # Band structure and DOS Nfo

locpot = Locpot.from_file("LOCPOT") # Electrostatic pot data

# Extract VBM and CBM

band_gap_data = vasprun.get_band_structure()

vbm = band_gap_data.get_vbm()["energy"]

cbm = band_gap_data.get_cbm()["energy"]

# Extract the vacuum level from LOCPOT

avg_potential = np.mean(locpot.data[0]) # Macroscopic average of potential

vacuum_level = max(avg_potential) # Maximum potential far from the slab

# Compute absolute band edge positions

vbm_absolute = vbm - vasprun.efermi + vacuum_level

cbm_absolute = cbm - vasprun.efermi + vacuum_level

# Convert to NHE scale

vbm_nhe = vbm_absolute + NHE_SHIFT

cbm_nhe = cbm_absolute + NHE_SHIFT

# PrResults

print(f"VBM (vs. NHE): {vbm_nhe:.2f} eV")

print(f"CBM (vs. NHE): {cbm_nhe:.2f} eV")

# Visualise

materials = ["#Material"]

band_edges = [(vbm_nhe, cbm_nhe)]

colors = ['blue', 'red']

fig, ax = plt.subplots(figsize=(8, 5))

for i, (vbm, cbm) in enumerate(band_edges):

ax.plot([i, i], [vbm, cbm], color="black", linewidth=2)

ax.scatter([i], [vbm], color=colors[0], label="Valence Band Edge" if i == 0 else "")

ax.scatter([i], [cbm], color=colors[1], label="Conduction Band Edge" if i == 0 else "")

# Reference lines fi

ax.axhline(y=0, color="green", linestyle="--", label="H₂/H₂O Level (0 eV)")

ax.axhline(y=1.23, color="purple", linestyle="--", label="O₂/H₂O Level (~1.23 eV)")

ax.set_xticks(range(len(materials)))

ax.set_xticklabels(materials, rotation=45)

ax.set_ylabel("Energy (eV vs. NHE)")

ax.set_title("Band Edge Positions vs. NHE")

ax.legend()

ax.grid(True)

plt.show()

Hope it helps 🙂

Dear Kostadin,

Thank you so much for your detailed answer and sharing the script. It is indeed helpful.

Could you please explain more for the bulk system? What is meant by referencing "relative to an external standard such as the vacuum level"?