Tutorial 01: Electronic Structure#
In this tutorial you will run the complete PAOFLOW post-processing workflow on bulk silicon: build a PAO Hamiltonian from atomic projections, interpolate the electronic band structure, and compute the density of states.
What you will learn
How Kohn–Sham wavefunctions are projected onto pseudo-atomic orbital (PAO) bases
How a real-space PAO Hamiltonian \(H(\mathbf{R})\) is built by Fourier transform
How interpolation produces smooth bands on arbitrarily dense k-grids
How the density of states is computed from interpolated eigenvalues
Tip
Try it locally: click here to download this tutorial as a notebook: tutorial01.ipynb.
The required .save directory assets can be obtained from: PAOFLOW Releases.
Required input#
The .save directory must contain:
data-file-schema.xml— crystal structure, k-point grid, and Kohn–Sham eigenvaluesatomic_proj.xml— projection amplitudes and overlap matrices
If atomic_proj.xml is not available, PAOFLOW can construct the projections internally from the pseudopotential data instead. That alternative is introduced in Section 2.
1. Initialise PAOFLOW#
PAOFLOW.PAOFLOW(...) reads data-file-schema.xml to load the crystal structure, k-point grid, and Kohn–Sham eigenvalues, then prepares all internal data structures. All output files are written to outputdir.
# Adjust these paths if needed
savedir = '../../examples/qe_examples/example01/silicon.save' # Path to the QE save directory
outputdir = 'output' # Path to the output directory for PAOFLOW results
from PAOFLOW import PAOFLOW
paoflow = PAOFLOW.PAOFLOW(
savedir=savedir,
outputdir=outputdir,
smearing='gauss',
npool=1,
verbose=True,
)
2. Read atomic projections#
PAOFLOW represents the electronic structure in a pseudo-atomic orbital (PAO) basis. The basis functions are localised orbitals derived from isolated-atom Kohn-Sham solutions for each angular momentum channel in the pseudopotential.
For silicon (\(3s^2\,3p^2\) valence configuration), this gives one \(s\) and three \(p\) orbitals per atom: eight PAOs for the two-atom unit cell.
The projection amplitudes quantify how much of Kohn-Sham band \(n\) at k-point \(\mathbf{k}\) is captured by PAO \(\alpha\):
The overlap matrix between PAOs on different sites is used in the subsequent Lowdin orthogonalisation:
paoflow.read_atomic_proj_QE()
If one does not want to run projwfc.x, PAOFLOW can instead construct the basis internally from the pseudopotential data and the radial functions stored in BASIS/. The rest of the workflow is unchanged.
For silicon, configuration='standard' is a practical starting point: it keeps the valence shells and adds extra radial functions that improve the representation of the low conduction bands. A larger extended basis is available when a wider conduction window is needed. See the Internal Basis Workflow for the full set of options.
# Alternative to read_atomic_proj_QE(): build the projections internally
# paoflow.projections(
# basispath='../../BASIS/',
# configuration='standard',
# )
3. Check projectability#
The projectability of band \(n\) at k-point \(\mathbf{k}\) measures how completely the PAO basis captures that band:
Values near 1 indicate a good representation; values well below 1 correspond to high-energy, plane-wave-like states not captured by localised orbitals.
Bands with \(p_{n\mathbf{k}}\) below the threshold (default 0.95) are excluded from Hamiltonian construction. For silicon, the four valence bands and four lowest conduction bands are well within the PAO window.
paoflow.projectability()
4. Build the PAO Hamiltonian#
With the projection amplitudes in hand, the k-space PAO Hamiltonian is assembled after Lowdin orthogonalisation of the PAO basis using \(S(\mathbf{k})\):
A discrete Fourier transform then maps this to real space:
The real-space Hamiltonian \(H(\mathbf{R})\) decays rapidly with distance; only a handful of Wigner-Seitz shells carry significant weight. This compact, first-principles tight-binding model is the basis for all subsequent interpolation: evaluating bands at any new k-point costs only a fast Fourier sum over these shells followed by an \(8\times8\) matrix diagonalisation.
paoflow.pao_hamiltonian()
5. Compute the band structure#
bands() evaluates the interpolated Hamiltonian from the real-space PAO matrix elements:
The calculation is performed along the standard FCC high-symmetry path (\(\Gamma\)-\(X\)-\(W\)-\(K\)-\(\Gamma\)-\(L\)-\(U\)-\(W\)-\(L\)-\(K\), selected by ibrav=2) at 2000 k-points.
Output:
bands_0.dat - contains the band energies along the k-path defined by ibrav=2
kpath_points.txt - contains the high-symmetry k-points and the number of points between them
paoflow.bands(ibrav=2, nk=2000)
Tip
Instead of using the default high-symmetry path selected by ibrav, you can explicitly define your own path with band_path and high_sym_points. This writes bands_custom_0.dat in outputdir, allowing you to keep both the default and custom-path results.
custom_path = 'L-G-X-W-K-G'
custom_points = {
'G': [0.0, 0.0, 0.0],
'L': [0.5, 0.5, 0.5],
'X': [0.5, 0.5, 0.0],
'W': [0.5, 0.25, 0.75],
'K': [0.375, 0.375, 0.75],
}
paoflow.bands(
ibrav=2,
band_path=custom_path,
high_sym_points=custom_points,
nk=2000,
fname='bands_custom',
)
One can then visualize the silicon band structure using bands_0.dat as shown below.
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
bands_file = Path(outputdir) / 'bands_0.dat'
kpath_file = Path(outputdir) / 'kpath_points.txt'
arr = np.loadtxt(bands_file)
k = arr[:, 0]
energies = arr[:, 1:]
labels = []
steps = []
with open(kpath_file) as f:
for line in f:
parts = line.split()
if len(parts) != 2:
break
labels.append(parts[0].replace('gG', r'$\Gamma$'))
steps.append(int(parts[1]))
# cumulative tick locations
tick_idx = np.cumsum([0] + steps[:-1])
# keep only valid tick positions
mask = tick_idx < len(k)
tick_idx = tick_idx[mask]
labels = np.array(labels)[mask]
tick_pos = k[tick_idx]
fig, ax = plt.subplots(figsize=(8, 6))
for i in range(energies.shape[1]):
ax.plot(k, energies[:, i], 'k-', lw=1.2)
ax.axhline(0, color='gray', ls='--', lw=0.5)
for x in tick_pos:
ax.axvline(x, color='gray', lw=0.5)
ax.set_xticks(tick_pos)
ax.set_xticklabels(labels)
ax.set_xlim(k[0], k[-1])
ax.set_xlabel('k-path')
ax.set_ylabel('Energy (eV)')
plt.tight_layout()
plt.show()
The silicon band structure displays the characteristic indirect gap: the valence band maximum sits at \(\Gamma\) while the conduction band minimum lies between \(\Gamma\) and \(X\). PBE-GGA underestimates the gap (~0.6 eV calculated vs 1.1 eV experimental) — a well-known limitation of the approximate exchange-correlation functional. One can see the low projectability bands that have been excluded from the PAO Hamiltonian in the conduction band region above ~7.4 eV.
# Example: custom path for silicon (FCC primitive coordinates)
custom_path = 'L-G-X-W-K-G'
custom_points = {
'G': [0.0, 0.0, 0.0],
'L': [0.5, 0.5, 0.5],
'X': [0.5, 0.5, 0.0],
'W': [0.5, 0.25, 0.75],
'K': [0.375, 0.375, 0.75],
}
paoflow.bands(
ibrav=2,
band_path=custom_path,
high_sym_points=custom_points,
nk=2000,
fname='bands_custom',
)
6. Interpolate the Brillouin-zone grid#
interpolated_hamiltonian() doubles the k-grid density by Fourier interpolation (e.g. 12x12x12 -> 24x24x24), providing a finer BZ sampling for accurate DOS integrations.
pao_eigh() diagonalises the interpolated Hamiltonian at every k-point on the denser grid:
The resulting eigenvalues and eigenvectors are stored for downstream transport and DOS calculations.
paoflow.interpolated_hamiltonian()
paoflow.pao_eigh()
7. Momentum matrix elements and adaptive smearing#
gradient_and_momenta() computes the k-gradient of the Hamiltonian:
From this, PAOFLOW evaluates the momentum matrix elements and group velocities used in transport:
adaptive_smearing() assigns each (band, k-point) a Gaussian broadening width proportional to:
Flat bands receive more broadening than dispersive ones, improving BZ integration accuracy without over-smearing sharp features.
paoflow.gradient_and_momenta()
paoflow.adaptive_smearing()
8. Compute and plot the density of states#
The electronic density of states is evaluated with adaptive Gaussian broadening at 1000 energy points spanning the full valence band and the lower conduction edge:
The silicon DOS shows a clear gap separating the valence and conduction manifolds. Van Hove singularities at band edges appear as steps or peaks. The total integrated weight under the valence DOS equals the number of valence electrons per cell (8 for Si).
Output:
dosdk_0.dat - contains the DOS as a function of energy
*_pdosdk_*.dat - contains the projected DOS onto each PAO as a function of energy. The postfix indicates the spin channel (0 for non-spin-polarised calculations). The prefix indicates the PAO index (e.g. 0 for the first PAO, 1 for the second, etc.).
paoflow.dos(emin=-12.0, emax=2.2, ne=1000)
# ----- total DOS -----
dos_file = Path(outputdir) / 'dosdk_0.dat'
dos_data = np.loadtxt(dos_file)
# ----- PDOS files -----
pdos_files = sorted(Path(outputdir).glob('*_pdosdk_*'))
fig, (ax_dos, ax_pdos) = plt.subplots(
2, 1, figsize=(8, 6), sharex=True, gridspec_kw={'height_ratios': [1, 2]}
)
# Total DOS
ax_dos.plot(dos_data[:, 0], dos_data[:, 1], lw=1.5, color='k', label='Total DOS')
ax_dos.axvline(0, color='gray', ls='--', lw=0.5)
ax_dos.set_ylabel('DOS')
ax_dos.legend()
# PDOS
for f in pdos_files:
data = np.loadtxt(f)
label = f.name.replace('_pdosdk_', ' ')
ax_pdos.plot(data[:, 0], data[:, 1], lw=1.2, label=label)
ax_pdos.axvline(0, color='gray', ls='--', lw=0.5)
ax_pdos.set_xlabel('Energy (eV)')
ax_pdos.set_ylabel('PDOS')
ax_pdos.legend(fontsize=8, ncol=2)
plt.tight_layout()
plt.show()
9. Further reading#
For a deeper explanation of concepts such as projectability, interpolation of Hamiltonian, please refer to the following papers: