43 — XSF and BXSF: periodic volumetric visualisation

Visualise crystal structures, electron densities, Bloch orbitals, and Fermi surfaces with XCrySDen’s native formats — directly from vibe-qc.

You’ll learn: how to write XSF files for crystal structures, electron densities, and single Bloch orbitals; how to write BXSF files for Fermi-surface visualisation; and how to open them in moltui, VESTA, and XCrySDen.

Why: XSF is the de-facto format for periodic volumetric data. Unlike Gaussian cube — which has no concept of a lattice — XSF stores the unit-cell vectors and tells the viewer to tile the cell. The result is correct periodicity with exactly one set of voxels, no supercell hacks. BXSF extends this to reciprocal space so XCrySDen can render Fermi surfaces directly.

Prerequisites: a working vibe-qc install; optional moltui for terminal viewing.

# Install moltui as a vibe-qc extra:
.venv/bin/pip install '.[viewer]'

The three verbs: structure, volumetric, bands

vibe-qc ships four writers in the XSF family:

Writer

What it writes

Best viewer

write_xsf_structure

Crystal-only XSF (PRIMVEC + PRIMCOORD)

VESTA, XCrySDen, moltui

write_xsf_mo

One Bloch orbital on the primitive cell

VESTA, XCrySDen

write_xsf_density

Total SCF density on the primitive cell

VESTA, XCrySDen, moltui

write_bxsf

Band energies on a regular k-mesh (Fermi surface)

XCrySDen, moltui

The first three work in real space; the last one works in reciprocal (k) space. All are callable as free functions from vibeqc; write_xsf_density is also wired into run_periodic_job(write_density=True) for one-line access.

Tip

For molecular volumetric output — cube files, molden, orbital visualisation — see tutorial 11 and the volumetric data user guide.

1. XSF crystal structure

The simplest XSF file carries only the crystal — lattice vectors and atom positions — with no volumetric data. Useful as a quick visual check on your PeriodicSystem before running the SCF.

import vibeqc as vq
import numpy as np

# NaCl rocksalt primitive cell.
a = 5.640 / 0.529177210903  # Å → bohr
system = vq.PeriodicSystem(
    3,
    [[a/2, a/2, 0], [a/2, 0, a/2], [0, a/2, a/2]],  # fcc primitive
    [vq.Atom(11, [0.0, 0.0, 0.0]),                    # Na at origin
     vq.Atom(17, [a/2, a/2, a/2])],                   # Cl at body-centre
)

vq.write_xsf_structure("nacl_primitive.xsf", system)
# → nacl_primitive.xsf — open in VESTA, XCrySDen, or moltui

The file contains a CRYSTAL block with PRIMVEC (lattice vectors in ångström) and PRIMCOORD (atom positions in ångström). vibe-qc converts from bohr automatically; XSF’s native unit is ångström.

For the one-call workflow, run_periodic_job writes a {stem}.xsf alongside the other output files:

basis = vq.BasisSet(system.unit_cell_molecule(), "sto-3g")

vq.run_periodic_job(
    system, basis,
    method="RHF",
    output="nacl",
    write_xsf_structure_file=True,   # default: True
    write_density=False,
)
# → nacl.out, nacl.system, nacl.molden, nacl.xsf (.structure.xsf
#   when write_density=True), nacl.POSCAR, nacl.cif

When write_density=True (see next section), the density lands in {stem}.xsf and the structure-only XSF is routed to {stem}.structure.xsf to avoid a suffix collision.

2. XSF density

Pass write_density=True to run_periodic_job and the job emits a {stem}.xsf with the SCF total electron density sampled on the primitive cell:

vq.run_periodic_job(
    system, basis,
    method="RHF",
    output="nacl",
    write_density=True,
    density_spacing_bohr=0.2,        # default
)
# → nacl.xsf — CRYSTAL block + DATAGRID_3D_density

The density ρ(r) = Σ D_{μν} χ_μ(r) χ_ν(r) is evaluated on a regular fractional-coordinate grid spanning one primitive cell. The BEGIN_DATAGRID_3D_density block holds the voxel values; VESTA renders them as a scalar isosurface, moltui as a raycast isosurface.

When both write_density=True and write_xsf_structure_file=True, the density lands in {stem}.xsf and the structure-only file is routed to {stem}.structure.xsf to avoid a suffix collision.

For programmatic access, use the standalone writer:

from vibeqc import write_xsf_density

# After a converged SCF, D_real is the real-space density.
write_xsf_density("nacl_density.xsf", system, basis, result.density_real)

3. XSF Bloch orbitals

To visualise a single crystalline (Bloch) orbital — ψ_{n,k}(r) — on the primitive cell, use write_xsf_mo. This is the periodic counterpart of write_cube_mo for molecules.

# Assumes a Γ-only SCF converged above.  At Γ the coefficients are
# real (up to an arbitrary sign) and k = (0, 0, 0).

homo = system.n_electrons() // 2 - 1

vq.write_xsf_mo(
    "nacl_homo.xsf", system, basis,
    np.asarray(result.mo_coeffs),       # C(k) — at Γ this is real
    np.array([0.0, 0.0, 0.0]),          # k = Γ
    band_index=homo,
    spacing_bohr=0.2,
    component="real",                    # Re ψ — the standard "MO picture"
)

write_xsf_mo constructs a primitive-cell grid, evaluates the Bloch orbital

\[ \psi_{n,k}(\mathbf{r}) = \sum_{\mathbf{T}} e^{i\mathbf{k} \cdot \mathbf{T}} \sum_\mu C_{\mu n}(\mathbf{k}) \, \chi_\mu(\mathbf{r} - \mathbf{T}) \]

at every voxel, and writes the result as a DATAGRID_3D block. The component= keyword selects what scalar to write:

Component

Scalar

Use when…

"real"

Re ψ

Orbitals at Γ or time-reversal-invariant k (standard MO picture)

"imag"

Im ψ

Gauge-paired view at general k

"abs"

|ψ|

Gauge-invariant magnitude

"density"

|ψ|²

Cross-k comparison; no phase ambiguity

For a full walk-through of Bloch-orbital evaluation — the evaluate_bloch_orbital kernel, the PrimitiveCellGrid factory, the cube-file workaround for molecular viewers, and the gauge- invariance considerations — see tutorial 22: periodic orbital cubes.

4. BXSF Fermi surfaces

BXSF (“band XSF”) stores band energies sampled on a regular 3D k-mesh across the Brillouin zone. XCrySDen / moltui render an isosurface at the Fermi energy — the Fermi surface — showing which bands cross the Fermi level and where in k-space.

What you need

write_bxsf expects:

  • A PeriodicSystem (for the reciprocal-lattice vectors),

  • A 4D energies array of shape (n_kx, n_ky, n_kz, n_bands) in Hartree (the writer converts to eV automatically),

  • An optional e_fermi in Hartree.

The k-mesh must be a regular fractional-coordinate grid over [0, 1) along each reciprocal-lattice direction — exactly what np.linspace(0, 1, nk, endpoint=False) produces.

End-to-end: NaCl rocksalt at pob-TZVP

This example runs a multi-k Hcore diagonalisation on a dense k-mesh and writes a BXSF file. NaCl is a wide-gap insulator, so no bands cross the Fermi level — the resulting “Fermi surface” will be empty. The workflow is identical for metals; just substitute a metallic system (e.g. an H chain at half-filling from tutorial 4) and use SCF-converged Fock matrices for production-quality surfaces.

import numpy as np
import vibeqc as vq
from vibeqc import (
    Atom, BasisSet, PeriodicSystem, LatticeSumOptions,
    compute_overlap_lattice, compute_kinetic_lattice,
    compute_nuclear_lattice, bloch_sum, diagonalize_bloch,
    write_bxsf,
)

# ---- 1. Build NaCl rocksalt -----------------------------------------
a = 5.640 / 0.529177210903            # Å → bohr
system = vq.PeriodicSystem(
    3,
    [[a/2, a/2, 0], [a/2, 0, a/2], [0, a/2, a/2]],
    [Atom(11, [0.0, 0.0, 0.0]),
     Atom(17, [a/2, a/2, a/2])],
)
basis = BasisSet(system.unit_cell_molecule(), "pob-tzvp")

# ---- 2. Build real-space lattice matrices (Hcore) -------------------
opts = LatticeSumOptions()
opts.cutoff_bohr = 12.0
opts.nuclear_cutoff_bohr = 25.0

S_lat = compute_overlap_lattice(basis, system, opts)
T_lat = compute_kinetic_lattice(basis, system, opts)
V_lat = compute_nuclear_lattice(basis, system, opts)

# ---- 3. Run a multi-k SCF for the Fermi energy ----------------------
# (For a production-quality BXSF you'd use the converged Fock
#  matrices here; this example uses Hcore to keep the recipe short.)
km_scf = vq.KPoints.monkhorst_pack(system, [4, 4, 4], symmetry=True)
opts_scf = vq.PeriodicSCFOptions()
opts_scf.lattice_opts.cutoff_bohr = 12.0
opts_scf.conv_tol_energy = 1e-6

result = vq.run_rhf_periodic_scf(system, basis, km_scf, opts_scf)

# Fermi energy: midpoint between the HOMO maximum and LUMO minimum
# across all k-points in the SCF mesh.
n_occ = system.n_electrons() // 2
homo = max(float(eps[n_occ - 1]) for eps in result.mo_energies)
lumo = min(float(eps[n_occ])     for eps in result.mo_energies)
e_fermi = (homo + lumo) / 2.0

# ---- 4. Sample band energies on a dense regular k-mesh --------------
nk = 16                                # 16×16×16 = 4096 k-points
kx = np.linspace(0, 1, nk, endpoint=False)
ky = np.linspace(0, 1, nk, endpoint=False)
kz = np.linspace(0, 1, nk, endpoint=False)
KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing="ij")

n_bands = basis.nbasis
energies = np.zeros((nk, nk, nk, n_bands))

for i in range(nk):
    for j in range(nk):
        for k in range(nk):
            # Convert fractional k → Cartesian k (bohr⁻¹).
            rec_lat = system.reciprocal_lattice()  # (3, 3), columns = b_i
            k_cart = KX[i, j, k] * rec_lat[:, 0] \
                   + KY[i, j, k] * rec_lat[:, 1] \
                   + KZ[i, j, k] * rec_lat[:, 2]
            H_k = bloch_sum(T_lat, k_cart) + bloch_sum(V_lat, k_cart)
            S_k = bloch_sum(S_lat, k_cart)
            sol = diagonalize_bloch(H_k, S_k)
            energies[i, j, k, :] = sol.eigenvalues

# ---- 5. Write BXSF --------------------------------------------------
write_bxsf("nacl_bands.bxsf", system, energies, e_fermi=e_fermi)
# → nacl_bands.bxsf — open in XCrySDen or moltui

The loop evaluates H(k) = T(k) + V(k) at every k-point of the dense 16³ mesh and diagonalises. For a converged SCF, replace the Hcore matrices with the real-space Fock lattice set from the SCF result; the sampling loop is identical.

Tip

For metals, the Fermi surface is non-empty. Substitute a metallic system — an H₂ chain at half-filling, or bulk Al — and the BXSF file will contain one or more bands that cross E_F. XCrySDen / moltui will render a surface where those bands equal E_F, exactly the textbook Fermi-surface picture.

Performance notes

  • A 16³ dense mesh diagonalising n_bands at each k-point takes O(nk³ · n_bands³) operations. For pob-TZVP on NaCl (42 basis functions, 2 atoms), a 16³ mesh completes in under a minute on a laptop.

  • File sizes scale as nk³ × n_bands × 6 bytes (ASCII float). 16³ × 42 bands ≈ 11 MB.

  • For production Fermi surfaces, 32³ or 48³ meshes are common; file sizes reach 100–300 MB.

5. Viewing in moltui

moltui reads XSF and BXSF natively. Launch it on any of the files written above:

moltui nacl_primitive.xsf          # crystal structure
moltui nacl.xsf                    # density isosurface
moltui nacl_homo.xsf               # Bloch orbital isosurface
moltui nacl_bands.bxsf             # Fermi surface / band isosurfaces

Keyboard controls (from the geometry / isosurface view):

Key

Action

h j k l / arrow keys

Rotate the structure

+ / -

Zoom in / out

b

Toggle periodic replication (XSF only)

m

Cycle view mode (geometry / orbitals / normal modes)

o

Toggle isosurface visibility

n / p

Navigate grids / bands (when multiple are present)

Tab / Shift-Tab

Cycle through bands (BXSF)

i / I

Decrease / increase isovalue

V

Open visual panel (isovalue, style, lighting)

e

Export PNG screenshot

q

Quit

Tip

Periodic replication (key b) tiles the primitive cell across the display. This is particularly useful for XSF structure and density files — you see the full crystal, not just one cell. The replication respects the lattice vectors, so off-diagonal cells (monoclinic, triclinic) tile correctly.

For more on moltui — including molecular-orbital viewing, normal-mode animation, and trajectory playback — see tutorial 27: viewing with moltui.

6. Viewing in VESTA / XCrySDen

VESTA

VESTA (free, cross-platform GUI) reads XSF files directly:

# macOS:
open -a VESTA nacl.xsf

# Linux:
vesta nacl.xsf

  • Structure: File → Open → select the .xsf. VESTA displays atoms + bonds with the unit-cell wireframe. Right-click → Style to switch to polyhedral / ball-and-stick / space-filling.

  • Volumetric data: VESTA automatically detects the DATAGRID_3D block. Edit → Volumetric Data → adjust the isosurface level and colour map.

  • Periodic replication: Objects → Boundary → set the number of replicated cells along each axis.

XCrySDen

XCrySDen (free, Linux/macOS) is the reference viewer for the format:

# XSF structure + volumetric data:
xcrysden --xsf nacl.xsf

# BXSF Fermi surface:
xcrysden --bxsf nacl_bands.bxsf

XCrySDen’s --bxsf mode opens the Brillouin-zone wireframe and renders an isosurface for every band that crosses the Fermi level. The wireframe is the reciprocal-lattice bounding box; the surfaces inside it are the Fermi sheets. Use Tools → k-path to overlay a band-structure path on the same display.

Tip

moltui vs VESTA vs XCrySDen. moltui is the terminal-native triage tool — fast, runs over SSH, no GUI. VESTA is the standard crystallographic viewer — publication-quality figures, polyhedral rendering, diffraction patterns. XCrySDen is the specialist for Fermi surfaces and band-structure overlay; its BXSF reader is the reference implementation. Pick the right viewer for the job.

Theory

XSF format conventions

The XSF format (Kokalj 2003) was designed for periodic crystals. A file consists of a CRYSTAL block containing:

  • PRIMVEC — three row vectors (in ångström) defining the primitive-cell lattice. vibe-qc converts from its internal bohr representation at write time.

  • PRIMCOORD — atom count, a format flag (1 for fractional coordinates), and one line per atom (Z  x  y  z in ångström).

  • DATAGRID_3D_ (optional) — a 3D scalar grid spanning a chosen real-space volume. The volume is specified by an origin (ångström) and three span vectors (ångström); for a primitive-cell density these are the lattice vectors. Voxel (0, 0, 0) is at the origin; voxel (n_a − 1, n_b − 1, n_c − 1) is one step short of the opposite corner. The voxel at the far corner — which coincides with the periodic image of (0, 0, 0) — is not stored. XSF / VESTA treat it as implicit, and including it introduces a spurious stripe at every cell boundary (the “off-by-one trap”).

The grid data runs with the first index varying fastest, then the second, then the third — the Fortran-order convention that vibe-qc honours via np.transpose(data, (2, 1, 0)).ravel().

BXSF format conventions

BXSF (“band XSF”) uses the same block-structured syntax but in reciprocal space:

  • The spanning vectors are the reciprocal-lattice vectors \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\) in 1/Å.

  • The origin is \((0, 0, 0)\) — the Γ point.

  • Band energies are in eV (vibe-qc converts from Hartree internally using 1 Ha = 27.211386 eV).

  • Each band is one BAND: block; the band index counts from 1.

  • The Fermi Energy in the BEGIN_INFO block is in eV.

XCrySDen renders the band structure as \(N_\text{bands}\) isosurfaces inside the BZ wireframe, one per band, at the Fermi energy. Bands that never reach \(E_F\) produce no visible isosurface; bands that cross \(E_F\) produce Fermi sheets — the textbook Fermi surface.

Real-space grids (XSF) vs k-space grids (BXSF)

XSF

BXSF

Domain

Real space \(\mathbf{r} \in\) unit cell

Reciprocal space \(\mathbf{k} \in\) BZ

Grid

Primitive-cell fractional coordinates

Monkhorst–Pack fractional coordinates

Scalar

ρ(r), ψ_{n,k}(r), or any f(r)

ε_n(k) — band energy

Isosurface meaning

Shape of the density / orbital

Shape of the Fermi surface

Span vectors

Lattice \(\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3\)

Reciprocal lattice \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\)

Units

Å for coordinates, a.u. for data

1/Å for reciprocal vectors, eV for energies

The real-space grid tells you where the electrons are; the k-space grid tells you which electrons cross the Fermi level and at what crystal momentum. Together they give the complete one-electron picture of a periodic system.

Fermi surfaces and the BZ

A Fermi surface is the set of k-points where

\[ \varepsilon_n(\mathbf{k}) = E_F \]

for some band \(n\). In an insulator, no band satisfies this condition — the occupied bands all lie below \(E_F\) and the virtual bands all lie above it, so the BXSF isosurface is empty. In a metal, at least one band crosses \(E_F\), and the resulting isosurface encloses the occupied region of k-space — the Fermi sea. The shape of the Fermi surface determines transport properties (conductivity, thermopower, magnetoresistance) and is the central object of metallic electronic structure.

The mesh density nk controls the smoothness of the isosurface. For publication-quality figures, 32³ or denser is recommended; 16³ is enough for orientation and triage.

References

  • Kokalj, A. “Computer graphics and graphical user interfaces as tools in simulations of matter at the atomic scale,” Comput. Mater. Sci. 28, 155 (2003). The XSF format specification. Updated at http://www.xcrysden.org/doc/XSF.html.

  • Bloch, F. “Über die Quantenmechanik der Elektronen in Kristallgittern,” Z. Phys. 52, 555 (1928). Bloch’s theorem — the foundation of crystalline orbitals and band theory.

  • Ashcroft, N. W.; Mermin, N. D. Solid State Physics. Holt, Rinehart and Winston (1976), Chapters 8–10 and 15. The standard textbook on Bloch states, band structures, and Fermi surfaces. Chapter 15 is the canonical introduction to Fermi-surface topology and measurement.

Next