Volumetric data: cube, XSF, and BXSF files

vibe-qc writes the interchange formats every solid-state / molecular viewer reads:

  • Gaussian cube (.cube) — molecular electron densities, molecular orbitals, and ESP. Read by VMD, Avogadro, PyMOL, ChimeraX, and moltui.

  • XCrySDen XSF (.xsf) — periodic crystal structures and volumetric data on a primitive-cell grid. Read by VESTA, XCrySDen, and moltui (since v0.8.x+ moltui).

  • XCrySDen BXSF (.bxsf) — band energies on a 3D k-mesh for Fermi-surface plots. Read by XCrySDen and moltui.

Quick-reference: which writer for which viewer

Viewer

Format

vibe-qc writer

VESTA, XCrySDen

XSF

write_xsf_structure, write_xsf_volume, write_xsf_density, write_xsf_mo, run_periodic_job(... write_density=True)

XCrySDen (Fermi surface)

BXSF

write_bxsf

VMD, Avogadro, PyMOL, ChimeraX

Cube

write_cube_density, write_cube_mo, write_cube_mos, run_job(... write_cube=...)

moltui (terminal)

XSF, BXSF, Cube

Any of the above — moltui reads all three

What you sample

A cube file is a uniformly-spaced 3D scalar field. For an electron density,

\[ \rho(\mathbf{r}) \;=\; \sum_{\mu\nu} D_{\mu\nu}\, \chi_\mu(\mathbf{r})\, \chi_\nu(\mathbf{r}), \]

vibe-qc evaluates \(\rho\) on every voxel and writes the values in Bohr\(^{-3}\) in the Gaussian-95 cube convention (header, atom block, then a flat float stream with the third lattice index running fastest). For a single molecular orbital,

\[ \phi_i(\mathbf{r}) \;=\; \sum_\mu C_{\mu i}\, \chi_\mu(\mathbf{r}), \]

the same evaluator is invoked, but the values written are \(\phi_i(\mathbf{r})\) (signed amplitude) rather than \(|\phi_i|^2\), so isosurface viewers (Avogadro, VMD, ChimeraX) can render both lobes of a \(\pi\) orbital. To render a density-of-states-like view of a single MO, square the file post-hoc in your viewer.

The file format itself is documented in Gaussian’s cube specification; XSF is described in the XCrySDen XSF Format reference and the binary BXSF variant in the same document § “BXSF file format”.

Molecular: cube files

import vibeqc as vq

mol = vq.Molecule.from_xyz("h2o.xyz")
basis = vq.BasisSet(mol, "6-31g*")
res = vq.run_rhf(mol, basis)

# Density: ρ(r) = Σ D_μν χ_μ χ_ν
vq.write_cube_density("rho.cube", res.density, basis, mol)

# A single MO (zero-based index)
homo = mol.n_electrons() // 2 - 1
vq.write_cube_mo("homo.cube", res.mo_coeffs, homo, basis, mol)

# Multiple MOs in one file (multi-volume cube; viewer flips between them)
vq.write_cube_mos("mos.cube", res.mo_coeffs, [homo - 1, homo, homo + 1, homo + 2],
                   basis, mol)

The grid is built automatically: an axis-aligned bounding box around the molecule plus 4 bohr of padding, with 0.2 bohr cubic voxels (~10⁶ voxels for a small molecule). Override with spacing= and padding=, or pass a fully custom :class:vibeqc.CubeGrid:

grid = vq.make_uniform_grid(mol, spacing=0.1, padding=6.0)
vq.write_cube_density("rho-fine.cube", res.density, basis, mol, grid=grid)

For UHF / UKS pass D = D_alpha + D_beta (the total density). A typical sanity check: integrating ρ over the box should return the electron count to ~1% on the default grid, exact in the limit.

Periodic: XSF files

vibe-qc offers three paths for XSF output, from simplest to most flexible:

1. Automatic via run_periodic_job

The simplest route. Pass write_density=True and the job writes {stem}.density.xsf alongside the other output files:

from vibeqc import Atom, PeriodicSystem, BasisSet, run_periodic_job

# NaCl rocksalt, sto-3g
a = 5.640 / 0.529177210903  # Å → bohr
system = PeriodicSystem(
    3,
    [[a, 0, 0], [0, a, 0], [0, 0, a]],
    [Atom(11, [0, 0, 0]), Atom(17, [a/2, a/2, a/2])],
)
basis = BasisSet(system.unit_cell_molecule(), "sto-3g")

run_periodic_job(
    system, basis,
    method="RHF",
    output="nacl",
    write_density=True,
    density_spacing_bohr=0.2,
)
# → nacl.out, nacl.system, nacl.molden, nacl.POSCAR,
#   nacl.xsf, nacl.cif, nacl.density.xsf

The density grid spans the primitive cell with density_spacing_bohr voxel spacing (default 0.2 bohr). The lattice vectors are written in Ångström; VESTA / XCrySDen / moltui will tile the cell automatically.

2. Structure-only via write_xsf_structure

vq.write_xsf_structure("crystal.xsf", system)

Writes the CRYSTAL block only (PRIMVEC + PRIMCOORD) — no volumetric data. Useful for quick visual sanity checks on a PeriodicSystem before running anything expensive. Opens in VESTA / XCrySDen / moltui as the bare unit cell.

3. Custom volumetric data via write_xsf_volume

import numpy as np

# Evaluate your scalar field on a primitive-cell grid …
rho = evaluate_density_on_grid(...)  # shape (n1, n2, n3), in bohr

vq.write_xsf_volume(
    "rho.xsf", system,
    data=rho, name="density",
    origin=np.zeros(3),          # Bohr — position of voxel (0,0,0)
    span=system.lattice.T,        # Bohr — grid spanning vectors
)

The span= rows define how far the grid extends — for a primitive-cell density, pass the lattice vectors. For a supercell, pass multiples thereof. The origin= is the position of voxel (0, 0, 0) in Bohr; use np.zeros(3) for a cell-aligned grid.

The most common XSF grid shape is the primitive-cell orbital evaluated via tutorial 22:

vq.write_xsf_mo("lih_bonding.xsf", system, basis, C, k, band_index=1,
                spacing_bohr=0.15)

4. XSF density writer

vq.write_xsf_density("density.xsf", system, basis, D,
                     spacing_bohr=0.2)

Shorthand for evaluating the total density D on a primitive-cell grid and writing it as a DATAGRID_3D_density block. Same grid conventions as write_xsf_mo.

Periodic: BXSF files (Fermi surfaces)

BXSF stores band energies sampled on a regular Monkhorst–Pack k-mesh in the Brillouin zone. Each band is one volumetric block; XCrySDen (or moltui) renders an isosurface at the Fermi energy to show the Fermi surface.

End-to-end: NaCl Fermi surface

import numpy as np
from vibeqc import Atom, PeriodicSystem, monkhorst_pack, run_rhf_periodic_scf, write_bxsf
from vibeqc import PeriodicSCFOptions, BasisSet

# Setup
a = 5.640 / 0.529177210903
system = PeriodicSystem(
    3,
    [[a, 0, 0], [0, a, 0], [0, 0, a]],
    [Atom(11, [0, 0, 0]), Atom(17, [a/2, a/2, a/2])],
)
basis = BasisSet(system.unit_cell_molecule(), "pob-tzvp")

# SCF on a 4×4×4 k-mesh
km = monkhorst_pack(system, [4, 4, 4])
opts = PeriodicSCFOptions()
opts.lattice_opts.cutoff_bohr = 12.0
opts.conv_tol_energy = 1e-6
result = run_rhf_periodic_scf(system, basis, km, opts)

# The result gives band energies per k-point.  Sample them on a
# denser mesh for smooth Fermi-surface rendering:
nk = 16
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")

# For each k-point in the dense mesh, evaluate band energies by
# diagonalising H(k).  Store as (nk, nk, nk, n_bands) in Hartree.
energies = np.zeros((nk, nk, nk, result.n_bands))
for i in range(nk):
    for j in range(nk):
        for k in range(nk):
            kpt = np.array([KX[i,j,k], KY[i,j,k], KZ[i,j,k]])
            H_k, S_k = build_fock_k(system, basis, result, kpt, opts)
            eigs = vq.diagonalize_bloch(H_k, S_k)
            energies[i, j, k, :] = eigs.eigenvalues

# e_fermi = midpoint between HOMO and LUMO at the converged k-mesh
e_fermi = result.fermi_energy

write_bxsf("nacl_bands.bxsf", system, energies, e_fermi=e_fermi)

Open nacl_bands.bxsf in XCrySDen (xcrysden --bxsf nacl_bands.bxsf) or moltui. Both render a Fermi-surface isosurface for each band at E = E_Fermi, with the Brillouin-zone wireframe as the bounding box. Cycle between bands to see which ones cross the Fermi level.

Units and conversions

  • BXSF energies are conventionally in eV. vibe-qc’s internal energies are in Hartree; write_bxsf converts automatically.

  • The reciprocal lattice spanning vectors are in 1/Å.

  • The k-mesh origin is (0, 0, 0); the three spanning vectors are the reciprocal lattice vectors \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3\).

  • Grid data runs with the first k-index varying fastest (same Fortran-order convention as XSF DATAGRID_3D).

Viewing in moltui

moltui reads XSF and BXSF natively alongside its existing cube support:

# Install (one-time)
./scripts/install_optional_tools.sh moltui

# View a crystal structure
moltui nacl.xsf

# View a density isosurface inside the unit cell
moltui nacl.density.xsf

# Cycle through Fermi surfaces (bands crossing E_F)
moltui nacl_bands.bxsf

Keyboard controls in moltui:

Key

Action

b

Toggle periodic replication

m

Cycle view mode (geometry / orbitals / normal modes)

o

Toggle isosurface visibility

o/p

Navigate grids/bands (when multiple are present)

V

Open visual panel (isovalue, style, lighting)

e

Export PNG screenshot

Set the XSF viewer as your default with the moltui command; for desktop-quality rendering, pipe the same file to VESTA or XCrySDen.

End-to-end example: HOMO / LUMO of water

A reproducible water HOMO/LUMO cube generation, ready to copy:

# input-water-cubes.py
from pathlib import Path
import vibeqc as vq

HERE = Path(__file__).parent

mol = vq.Molecule([
    vq.Atom(8, [ 0.0,  0.00,  0.00]),   # bohr
    vq.Atom(1, [ 0.0,  1.43, -0.98]),
    vq.Atom(1, [ 0.0, -1.43, -0.98]),
])
basis = vq.BasisSet(mol, "6-31g*")
result = vq.run_rhf(mol, basis)

# Closed-shell H2O has 10 electrons → 5 occupied MOs.
homo = mol.n_electrons() // 2 - 1            # MO index 4
lumo = homo + 1                              # MO index 5

vq.write_cube_density("rho.cube", result.density, basis, mol)
vq.write_cube_mo("homo.cube", result.mo_coeffs, homo, basis, mol)
vq.write_cube_mo("lumo.cube", result.mo_coeffs, lumo, basis, mol)

Drag the three files into VMD / Avogadro / ChimeraX. rho.cube shows the total electron density (canonical isovalue 0.05 e/bohr\(^3\) for a textbook bond-density picture); homo.cube and lumo.cube render with isovalue ±0.05 a.u. to show the lobes (red = negative amplitude, blue = positive — viewer-dependent).

Sanity check: \(\int \rho \,d^3r\) should integrate to 10 (vibe-qc’s test suite verifies this to ~1% on the default grid). \(\int |\phi_i|^2 \,d^3r\) should integrate to 1 for any occupied MO.

Validation

tests/test_cube.py checks:

  • the writer header parses back as expected (atom count, voxel vectors, grid shape, multi-value flag),

  • ρ(r) integrates to the electron count on a fine grid,

  • single-MO files have ⟨φ|φ⟩ ≈ 1.

tests/test_xsf.py checks:

  • ångström conversion is correct for lattice and atom positions,

  • XSF DATAGRID traversal order (first index runs fastest, then j, then k) — required by VESTA/XCrySDen,

  • BXSF conventions: 1/ångström reciprocal lattice and eV energies.