Volumetric data: cube, XSF, and BXSF files

vibe-qc’s native visualization format is QVF (.qvf): a single archive carrying the structure, density, orbitals, and basis together, opened in vibe-view. run_job and run_periodic_job write it by default (output_qvf=False disables it), and it is the recommended way to visualize vibe-qc results: see output files, the QVF tech spec, and The QVF file format, end to end. For periodic QVF consumers, the downloadable chi-ccm-b-qvf fixture bundle provides validated 1D and 3D finite-BvK archives with density, orbital, replicated-cell, and Wannier-centre overlay captures.

This page covers the interchange formats, for when an external solid-state or molecular viewer is what you need:

  • 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 (K. Momma, F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011), doi:10.1107/S0021889811038970), XCrySDen (A. Kokalj, J. Mol. Graphics Modell. 17, 176 (1999), doi:10.1016/S1093-3263(99)00028-5), 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

vibe-view (native)

QVF

high-level runners by default; keep output_qvf=True explicit when required

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 the orbital density \(|\phi_i|^2\) 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 ASCII Band XSF (BXSF) syntax in the same document under “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:

# doc-audit: skip - continues the molecule, basis, and result constructed above
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}.xsf alongside the other output files. If structure output is also enabled, the structure-only sibling is {stem}.structure.xsf:

import numpy as np
import vibeqc as vq

# NaCl rocksalt primitive cell, sto-3g.
a = 5.640 / 0.529177210903  # angstrom to bohr
a1 = np.array([0.0, a / 2, a / 2])
a2 = np.array([a / 2, 0.0, a / 2])
a3 = np.array([a / 2, a / 2, 0.0])
lattice = np.column_stack((a1, a2, a3))
cl_position = 0.5 * (a1 + a2 + a3)

system = vq.PeriodicSystem(
    3,
    lattice,
    [vq.Atom(11, [0.0, 0.0, 0.0]), vq.Atom(17, cl_position)],
)
basis = vq.BasisSet(system.unit_cell_molecule(), "sto-3g")

vq.run_periodic_job(
    system,
    basis,
    method="RHF",
    output="nacl",
    write_density=True,
    density_spacing_bohr=0.2,
)
# The density is nacl.xsf; the structure-only sibling is nacl.structure.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.

The same spacing sets the density volume in the .qvf archive. That volume is written only if the grid sum over the primitive cell, sum(rho) * dV, equals the cell’s electron count (excluding ECP cores) within 0.005 electrons. A uniform grid integrates a periodic function with an error made of the function’s Fourier components at the grid’s reciprocal vectors (aliasing; Trefethen and Weideman, SIAM Rev. 56, 385 (2014), Sec. 2). The all-electron cores of second-row and heavier atoms still have large components there at 0.2 bohr, so the sum can be off by whole electrons in either direction. Ne/STO-3G in a 7-bohr cube gives 11.98 instead of 10 electrons at 0.2 bohr and 10.00002 at 0.05 bohr, with the same SCF energy. When the check fails, the archive is still written without the volume, a warning names the integral, and the file’s outcome records the error under volume.density. Decrease density_spacing_bohr until the volume is written. Heavier cores need finer grids, and the number of grid points grows with the cube of the refinement.

2. Structure-only via write_xsf_structure

# doc-audit: skip - continues the periodic system constructed above
vq.write_xsf_structure("crystal.xsf", system)

Writes the structure 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.

The block keyword follows system.dim: CRYSTAL for bulk, SLAB for a 2D system, POLYMER for a 1D one. For dim < 3 the padded PRIMVEC rows are bookkeeping, not cell edges; see XSF format conventions.

3. Custom volumetric data via write_xsf_volume

# doc-audit: skip - schematic custom field evaluated by user code
import numpy as np

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

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 Periodic orbital cubes and XSF files:

# doc-audit: skip - continues explicit periodic lattice objects from the linked workflow
vq.write_xsf_mo("lih_bonding.xsf", system, basis, C, k, band_index=1,
                spacing_bohr=0.15)

4. XSF density writer

# doc-audit: skip - continues an explicit lattice-resolved density workflow
vq.write_xsf_density("density.xsf", system, basis, density_real,
                     spacing_bohr=0.2)

This lower-level writer requires density_real to be an explicit real-space LatticeMatrixSet, such as the output of vq.real_space_density_from_kpoints. A generic periodic result does not promise a density_real attribute. Use the automatic runner route above unless your workflow already owns the lattice-resolved density.

Periodic: BXSF files (regular k-space band grids)

BXSF is the ASCII Band XSF format. It stores band energies sampled on a regular three-dimensional reciprocal-space grid. Each band is one ASCII BAND: block. XCrySDen renders a constant-energy isosurface, normally at the Fermi energy. Only a self-consistent metallic spectrum with a physically defined chemical potential produces a scientific Fermi surface.

End-to-end: NaCl Hcore format demonstration

The complete recipe below uses public vibe-qc APIs to build an Hcore band grid for a true two-atom rocksalt primitive cell. Hcore contains only kinetic and nuclear-attraction terms. The file is useful for testing BXSF layout and viewer operation, but it is not an RHF or DFT band structure and must not be interpreted as a NaCl Fermi surface.

import numpy as np
import vibeqc as vq

# 1. True NaCl rocksalt primitive cell. Lattice vectors are columns.
a = 5.640 / 0.529177210903  # angstrom to bohr
a1 = np.array([0.0, a / 2, a / 2])
a2 = np.array([a / 2, 0.0, a / 2])
a3 = np.array([a / 2, a / 2, 0.0])
lattice = np.column_stack((a1, a2, a3))
cl_position = 0.5 * (a1 + a2 + a3)

system = vq.PeriodicSystem(
    3,
    lattice,
    [vq.Atom(11, [0.0, 0.0, 0.0]), vq.Atom(17, cl_position)],
)
basis = vq.BasisSet(system.unit_cell_molecule(), "sto-3g")

# 2. Build the real-space Hcore and overlap lattice blocks once.
opts = vq.LatticeSumOptions()
opts.cutoff_bohr = 12.0
opts.nuclear_cutoff_bohr = 25.0
overlap_real = vq.compute_overlap_lattice(basis, system, opts)
kinetic_real = vq.compute_kinetic_lattice(basis, system, opts)
nuclear_real = vq.compute_nuclear_lattice(basis, system, opts)

# 3. Sample [0, 1) in fractional reciprocal coordinates.
nk = 8
mesh = (nk, nk, nk)
energies = np.empty((*mesh, basis.nbasis), dtype=float)
reciprocal = np.asarray(system.reciprocal_lattice())

for index in np.ndindex(mesh):
    k_fractional = np.asarray(index, dtype=float) / np.asarray(mesh)
    k_cartesian = reciprocal @ k_fractional
    overlap_k = vq.bloch_sum(overlap_real, k_cartesian)
    hcore_k = (
        vq.bloch_sum(kinetic_real, k_cartesian)
        + vq.bloch_sum(nuclear_real, k_cartesian)
    )
    overlap_k = 0.5 * (overlap_k + overlap_k.conj().T)
    hcore_k = 0.5 * (hcore_k + hcore_k.conj().T)
    solution = vq.diagonalize_bloch(hcore_k, overlap_k)
    energies[index] = np.asarray(solution.energies, dtype=float)

# 4. Put the display reference midway between the Hcore band edges.
n_occupied = system.n_electrons() // 2
valence_maximum = float(energies[..., n_occupied - 1].max())
conduction_minimum = float(energies[..., n_occupied].min())
e_reference = 0.5 * (valence_maximum + conduction_minimum)
hcore_edge_separation = conduction_minimum - valence_maximum

if hcore_edge_separation <= 0.0:
    print("Hcore manifolds overlap; displayed surfaces are nonphysical")

vq.write_bxsf(
    "nacl_hcore_bands.bxsf",
    system,
    energies,
    e_fermi=e_reference,
)

NaCl is physically an insulator, but that fact does not guarantee that this non-SCF Hcore model has a positive indirect gap. The chosen reference is only the midpoint between its occupied and virtual band edges. If those manifolds overlap, the viewer may draw an Hcore isosurface; it is still not a physical Fermi surface. For a production metallic surface, obtain the dense-grid eigenvalues and chemical potential from a route-specific converged SCF workflow. Generic periodic result objects do not expose a universal real-space Fock lattice set.

Units and conversions

  • BXSF is plain ASCII, not a binary container.

  • 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.xsf

# Inspect the Hcore band grid (NaCl has no physical Fermi surface here)
moltui nacl_hcore_bands.bxsf

Keyboard controls in moltui:

Key

Action

b

Toggle periodic replication

m

Cycle view mode (geometry / orbitals / normal modes)

o

Toggle isosurface visibility

n / 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.