Quantum chemistry needs a modern file format

Diagram showing today's fragmented output files on the left and a single .qvf container with typed sections on the right

Quantum chemistry needs a modern file format

Diagram showing today's fragmented output files on the left and a single .qvf container with typed sections on the right
From a directory of slices to a single container

Quantum chemistry calculations in 2026 produce rich, many layered results. A single run can generate atomic structures, volumetric grids for densities and orbitals, band structures, vibrational modes, spectra, trajectories, and detailed provenance describing how the calculation was performed. In principle, all of this belongs together as one coherent dataset.

In practice, it does not.

Instead, results are fragmented across a collection of loosely related files: cube files for densities, XSF for periodic systems, XYZ for trajectories, vendor specific binaries for restart data, and a mix of text, JSON, or ad hoc logs for energies and convergence. Each file captures a slice of the calculation, but none captures the whole. There is no shared manifest, no consistent schema, and often no reliable record of provenance.

This fragmentation is not just inconvenient. It is a structural limitation that has persisted for decades.

The real problem: an aging ecosystem

QC visualization landscape: tools sorted by last release year, showing active, stalled, and effectively dead categories
The freeze line in QC visualization tools

The issue is not only the number of formats, but the state of the tools that use them.

A look at the visualization tools that quantum chemists actually reach for shows a field that is partially active and partially frozen in place. Before listing them, one distinction matters. A pure orbital or density viewer exists to look at the output of a QC calculation: orbitals, densities, electrostatic potentials, NTOs, IBOs, ELF, and so on. A molecule editor with a built in renderer exists primarily to build and edit structures, and happens to visualize cubes as a secondary capability. Avogadro 2 is the latter. So is Gabedit. So is GaussView. So is every vendor GUI tied to a specific QC engine. The distinction matters because the editor category has steady commercial backing (Schrödinger, Gaussian Inc., SCM, Wavefunction, TURBOMOLE GmbH), while the pure viewer category is overwhelmingly academic and overwhelmingly thin.

Pure orbital and density viewers, still maintained

Tool Latest release License
Multiwfn Apr 2026 (rolling), 3.8 stable (Jan 2026) Free, permissive, academic and commercial
Pegamoid Aug 2025 (last commit), 2.12.4 active 2024 to 2025 GPL
MOrbVis Active 2025, WebGPU based Open source
IQmol Jul 2024 Open source, Q-Chem affiliated
Molden / gmolden Oct 2023 Free, proprietary terms
IboView Jan 2022 Free, custom license
QMForge Feb 2020 GPL. Slow to stall.
Speck 2017 to 2020 era, sporadic MIT

Editors and full GUIs that include a renderer, still maintained

Tool Latest release License
Avogadro 2 Apr 2026 BSD 3 clause / GPL v2, open source
AMSview (AMS GUI) 2026.1 Commercial (SCM)
TmoleX 2024.1 era Commercial (TURBOMOLE, COSMOlogic)
Maestro / Jaguar GUI 2026.1 era Commercial (Schrödinger)
Chemcraft 2024, ongoing Commercial. Free non saving demo.
wxMacMolPlt Jan 2024 GPL
Spartan Spartan’24 Commercial (Wavefunction)
Vipster 1.20.0 era, active 2023 to 2024 GPL
Gabedit Jul 2021 GPL
GaussView 2019 Commercial (Gaussian, Inc.)

Web and JavaScript stack, still maintained

Tool Latest release License Notes
Jmol / JSmol Mar 2026 LGPL 2.0 Reads cube, molden, fchk
3Dmol.js Jan 2026 BSD Loads cube and VASP data, isosurfaces
Mol* (mol star) Active 2024 to 2026 MIT Successor to NGL and LiteMol
NGL Viewer Core 2.x maintained, NGLVieweR 1.4.0 (Nov 2024) MIT Primarily structural biology, reads cube
Miew (EPAM) Active MIT Mainly biomolecular
iCn3D (NCBI) Active NIH Mainly PDB, has cube support

Periodic, crystallographic, and materials oriented

Tool Latest release License
Mercury (CCDC) 2024.3.0 era Free base, full version needs CSD license
VESTA Aug 2022 Free for academic and non commercial
Olex2 1.5 era, active Free with citation
CrystalMaker 11.x era Commercial
Diamond 5.x (Crystal Impact) Commercial
XCrySDen Oct 2019 GPL v2

Biomolecular tools routinely used for orbital cubes

Tool Latest release License
PyMOL Feb 2025 Schrödinger commercial. Open source build available.
UCSF ChimeraX 2024 to 2025 Free academic, commercial fee
UCSF Chimera 1.19 (2024). Considered legacy, succeeded by ChimeraX. Free academic
VMD 1.9.3 stable (Nov 2016). 1.9.4 RC and 2.0.0 alpha Unix in 2025. Free non commercial, custom license

The graveyard

A larger group of QC visualization tools has either explicitly retired, last released a decade or more ago, or sits in maintenance limbo with sporadic patches. Several of these were dominant in their era.

Tool Last release Status
MOLEKEL 5.4 (Aug 2009) Abandoned. Code on GitHub. GPL.
gOpenMol 3.00 (2005, Windows binary 2010) Dead. Original CSC Finland page gone.
Viewmol 2.4.1 mid 2000s, last activity ~2009 Dead. GPL.
RasMol 2.7.5.1 (Jul 2009) Effectively dead, but the protein chestnut. GPL.
ArgusLab 4.0.1 mid 2000s Dead. Long promised Qt and iPad ports never appeared.
QuteMol 0.4.1 (Jun 2007) Dead. GPL with citation.
ECCE 7.0 official PNNL (Aug 2013), fork 7.3.2 beta ~2017 Effectively dead. ECL 2.0.
GaussSum 3.0.2 (2013) Effectively dead. GPL. Built on cclib.
MoleCoolQt 2.4 era, sporadic commits Largely dormant. Charge density oriented.
Luscus Source side ~2016 to 2018 Stagnant. Academic Free License v3.
MOLDA Early 2000s Dead. Site gone for years.
PyMOlyze Early 2000s Dead. Rolled into QMForge philosophy.
MOLMOL 2K.2 era (ETH Zurich) Dead. NMR oriented but used for structures.
HyperChem 8.0.10 (2011) Effectively dead. Still sold by Hypercube.
Cerius2 Last Accelrys release ~2007 Dead. Replaced by Materials Studio.
InsightII Last Accelrys release ~2005 Dead. The SGI workstation era classic.
MacMolPlt (original) Pre 7.0, mid 2000s, Carbon Long superseded by wxMacMolPlt.

What the landscape shows

The actively maintained pure orbital viewers are a short list: Molden, IboView, Multiwfn, Pegamoid (for OpenMolcas users), and MOrbVis (a recent WebGPU upstart). IQmol fits here too if classified as a viewer rather than a builder. That is essentially the full set. Everything else is either an editor with a renderer attached, a primarily biomolecular tool repurposed for cube files, or a vendor GUI tied to a specific QC engine. VMD straddles categories. The stable 1.9.3 from 2016 is what most people still run, while 1.9.4 release candidates and the 2.0.0 alpha have been in slow public development for nearly a decade.

Several widely used classics have stalled. XCrySDen last shipped in 2019, Gabedit in 2021, IboView in early 2022, VESTA in mid 2022. The truly free and open source options that a new project can realistically embed or extend are a narrow set: Avogadro 2 (BSD), Multiwfn (permissive even for commercial use), Jmol (LGPL), Gabedit (GPL), wxMacMolPlt (GPL), and Vipster (GPL). Tools that are free to download but carry restrictive or proprietary terms include Molden, VMD, IboView, and VESTA.

One format to rule them all

What this landscape calls for is not yet another viewer. It is a common file format that all of these tools, alive and yet to be written, can read.

QVF is meant to be that format. One container that consolidates what is currently spread across cube files, XSF, MOLDEN format, fchk, XYZ, ad hoc JSON output, and a long list of vendor binaries. One manifest with the structure, the volumetric fields, the bands, the spectra, the vibrations, the trajectory, and the full provenance, packaged so that any consumer can read what it understands and ignore the rest. One specification, openly licensed, that any code or viewer can implement without coordination.

The historical reason this never happened is that adding format support to a dozen aging tools was a multi year effort with no single party in a position to drive it. That argument is no longer the obstacle it was. Adding a reader or writer for a clearly specified format is now the kind of task a coding agent finishes in an afternoon, given the schema, a few example files, and the test suite that ships with the spec. Quantum chemistry is one of the fields where this leverage is most consequential. The codes are large, the maintainers are few, and the time saved is exactly the time that did not previously exist.

Many of these tools were developed by PhD students or small research groups to solve a specific problem, published alongside a paper, and then left behind as those researchers moved on or shifted focus. Over time, their dependencies age, their build systems break, and their assumptions about input formats drift out of sync with modern codes.

At the same time, quantum chemistry packages themselves evolve. Output formats change subtly, or sometimes significantly, between versions. In other cases, formats remain nominally the same while their contents shift in undocumented ways. Parsing output becomes brittle, requiring constant patching and format specific workarounds.

The result is a fragile ecosystem. Visualization tools support different and only partially overlapping formats. No single format captures all relevant data. Output parsing is error prone and version dependent. Reproducibility suffers because provenance is incomplete or lost. Moving data between tools requires manual intervention.

Even when formats are well established, they are difficult to extend. Adding new data types or metadata often breaks compatibility, so formats stagnate. The path of least resistance becomes creating yet another file type rather than improving an existing one.

The cost of text based legacy formats

A particularly visible symptom of this legacy is the continued reliance on text based formats for large numerical data.

Formats like Gaussian cube store volumetric grids as ASCII text. This was reasonable when grid sizes were small. It is no longer defensible. A 200³ grid contains around 8 million values and can exceed 180 MB as text, while the same data stored as binary float32 is about 32 MB and loads orders of magnitude faster. Larger grids routinely cross into gigabyte scale text files.

Bar chart comparing the on-disk size of a 200³ volumetric grid in four encodings: Cube ASCII 180 MB, float64 binary 64 MB, float32 binary 32 MB, float32 plus zstd 15 MB
Same data, four encodings

These formats persist not because they are efficient or well suited to modern workloads, but because they are entrenched.

What a modern format needs to do

A modern file format for quantum chemistry must reflect how the field actually works today.

It should package the entire result of a calculation, including structure, grids, spectra, trajectories, and provenance, into a single file that can be shared in one step. It must support random access so that tools can read only the data they need without loading everything. It should be machine readable, explicitly typed, and unambiguous, with clear units and schema validation.

Equally important, it must be designed to evolve. New types of data should be addable without breaking existing tools. Unknown data should be safely ignored but not lost. The format should not depend on heavy, specialized libraries, so that it can be implemented across languages and environments, including lightweight scripts and browser based tools.

In short, it should be a format that supports both stability and change.

Enter QVF

QVF (Quantum Visualization Format) is designed to meet these requirements.

A QVF file is a self contained ZIP archive. At its core is a single manifest.json that describes all data contained in the file. The rest of the archive consists of typed sections: structure, volumetric data, band structures, spectra, vibrations, trajectories, and provenance.

Large numerical arrays are stored as raw binary (typically little endian float32), making them compact and fast to load. Smaller or structured data is stored as JSON. Every section is explicitly labeled with a “kind” drawn from a controlled vocabulary, allowing tools to declare what they support and ignore what they do not.

Because ZIP provides a central directory, any section can be accessed directly without scanning the entire file. A visualization tool that only needs the structure reads a few kilobytes, even if the file also contains hundreds of megabytes of volumetric data.

The format is fully self describing. Units are explicit. Provenance is mandatory. Every binary section includes a sha256 checksum for integrity. Provenance also carries an explicit agent model trail for AI driven workflows, so a calculation produced by an agent loop records the same audit information as one run by hand at a terminal.

Cross section of a .qvf archive: manifest.json at the top followed by typed section directories for structure, volumes, bands, spectra, vibrations, trajectory, and vendor namespaces, with a consumer reading via the central directory
Inside a .qvf file

Designed for an ecosystem, not a single tool

One of the central design goals of QVF is to enable an ecosystem rather than a single implementation.

A QVF reader does not need to understand everything in the file. It reads the manifest, identifies the section types it supports, and processes only those. Unsupported sections are not errors. They are simply reported as present but unused. This allows the format to grow over time without breaking older tools.

Vendors and projects can introduce their own extensions under a namespaced convention (x_<vendor>.*) without requiring coordination. At the same time, a shared core vocabulary ensures interoperability for common data types.

Capability is declared by kind, not by version. There is no level 1 or level 2 tiering. A minimal structure viewer supports structure. A bands plotter supports structure plus bands. SemVer rules pinned to kinds keep older readers working as the registry grows: minor versions add new optional kinds, major versions are rare and deliberate.

This model reflects a key reality: no single tool will ever cover all use cases. The format must allow many tools, with different capabilities, to operate on the same file.

Matrix showing which QVF section kinds each consumer type supports: structure only viewer, bands plotter, spectra tool, reference viewer, and validator
Capability is declared by kind, not by version

Scope and practical limits

QVF is intentionally scoped as a visualization and analysis container. It stores evaluated data, including grids, spectra, and trajectories, but does not attempt to encode the full internal state of a quantum chemistry calculation, such as basis set integrals or wavefunction coefficient matrices. Including those would dramatically increase file size and blur the line between data container and simulation engine.

Similarly, extremely large time series of volumetric data is out of scope for the initial version. These constraints keep file sizes manageable, typically tens to a few hundred megabytes, while covering the vast majority of practical use cases.

What vibe-qc commits to

vibe-qc will adopt QVF as its primary visualization output format and maintain the QVF specification going forward. The writer is implemented as part of vibe-qc’s existing output infrastructure, so every converged calculation produces a .qvf file alongside the standard log output. A validation tool, qvf-validate, ships with vibe-qc and lets any other producer check that its output is conformant.

The specification is open and the license is MPL 2.0, the same license as vibe-qc itself. SemVer governance, growth of the kind registry, and the consumer contract are vibe-qc maintained responsibilities. Vendor extensions remain entirely independent of the central spec.

A reference viewer is on the roadmap. It will ship with vibe-qc and support the full v1 section vocabulary: structure with unit cell and optional symmetry overlay, isosurface rendering for the volumetric kinds with viewer suggested defaults, animated vibrational modes, IR and Raman spectra with adjustable broadening, band structure plots along the declared k path, and trajectory playback. The viewer is not written yet. The format is written first on purpose, so that the viewer can be built against a stable specification and other codes and tools can adopt the format in parallel rather than waiting on a single reference implementation.

We strongly encourage other quantum chemistry codes and visualization tools to adopt QVF. The writer side is light: a few hundred lines of code on top of any modern ZIP library. The reader side is lighter still. Parse the manifest, decompress the binary members you care about, render. Anything outside the kinds you support is something you label as present but unused and move on from.

Feedback, adoption notes, and proposals for new kinds are welcome at mpei@vibe-qc.com.

Toward a better default

The goal of QVF is not to replace every existing format overnight, but to provide a better default.

A single file that contains everything needed to inspect, share, and reproduce a calculation. A format that is fast, explicit, and extensible. A foundation that multiple tools can build on without constant reinvention.

For a field that has long relied on fragmented, aging infrastructure, that shift is overdue.



The full v0.4 QVF specification, including the JSON Schema sketch, the kind registry, and the considered-and-rejected alternatives, lives in the vibe-qc repository at docs/design_qcv_format.md.

Eight releases in nine hours: the vibe-qc v0.4.0 → v0.4.7 bug arc

Timeline of eight releases from v0.4.0 to v0.4.7 over nine hours forty-one minutes

vibe-qc tagged v0.4.0 at 14:20 on April 27. By midnight on April 28 it had tagged seven more releases. All eight share the codename “Schrödinger’s Llama.” None of the patches fixed quantum chemistry bugs. Every single one fixed something that only surfaced when a real user on a non-dev machine tried to install and run the code.

This is the story of those nine hours and forty-one minutes, and what they taught about the gap between a working release and a usable one.

Timeline of eight releases from v0.4.0 to v0.4.7 over nine hours forty-one minutes
Eight releases in 9h41m. Green: release milestones. Pink: Linux or build fixes. Blue: docs and workflow. Every patch fixed something invisible on the developer’s macOS machine.

v0.4.0: what shipped

The v0.4.0 tag represented the culmination of the development arc described in the previous post: end-to-end periodic SCF across all method and spin combinations, effective-core potentials via libecpint, Fermi-Dirac smearing, Saunders-Hillier level shifting, quadratic SCF fallback, DFT-D3(BJ) dispersion, and 25 tutorials with theory sections and verified citations. 879 tests passing, one xfailed, Sphinx building clean with -W. The GitLab CI pipeline auto-deployed the docs site on push. The release branch was public.

What the pre-flight checks cannot tell you is what happens when someone on a different OS, with a different CMake version, on a different Linux distribution, follows your quickstart from a cold start.

v0.4.1: the build system has opinions (+7h 24m)

The first stress test came from running setup_native_deps.sh on Arch Linux. Three independent failures inside an hour:

Eigen 5 detection. The check in build_libint.sh used ls a b to test for header existence. On Arch, Eigen 5 installs only at /usr/include/eigen3/ with no legacy symlink at /usr/include/Eigen/. The ls call returned non-zero, the check reported Eigen missing, and the build failed — even after a fresh pacman -S eigen. Fix: replace with pkg-config --exists eigen3 plus a [ -e ... ] short-circuit.

CMake 4.x compatibility. libxc 7.0.0, FFTW 3.3.10, and several libecpint vendored dependencies still declare cmake_minimum_required(VERSION <3.5), which CMake 4.2 on Arch rejects outright without an override. Fix: add -DCMAKE_POLICY_VERSION_MINIMUM=3.5 to every vendored-dep CMake invocation.

CMP0167 Boost warning. libint 2.13.1 uses the deprecated find_package(Boost) interface, which produces a wall of CMake dev warnings that frightened users into thinking something was broken. Fix: add -Wno-dev.

The lesson is not specific to these three bugs. It is that build systems accumulate assumptions about the developer’s environment — which headers live where, which CMake version the project owner happens to be running, which deprecated interfaces are still tolerated on the local machine — and those assumptions are invisible until someone runs the script on a different machine.

v0.4.2: “I built the docs” does not mean “I followed them” (+8h 2m)

After getting past the build, I followed my own quickstart on a clean box. Three things blocked me before the first SCF calculation ran.

The landing page pointed users to the pre-orchestrator scripts: ./scripts/build_libint.sh && ./scripts/setup_basis_library.sh. Those two scripts leave libxc, spglib, FFTW, and libecpint unbuilt. The correct entry point is ./scripts/setup_native_deps.sh, which runs everything. The fix was straightforward; the fact that the wrong path survived the docs CI is the point. CI proves the docs are warning-free. Only a real user proves the docs make sense.

The “your first calculation” snippet had no filename, no run command, and no mention that import vibeqc requires the virtualenv’s Python specifically. A reader who installed into a venv and then typed python water.py got ModuleNotFoundError and no obvious next step. The fix: explicit “save as water.py, run with .venv/bin/python water.py,” plus a tip box covering the activate pattern and the failure mode.

There was also no single page that answered “how do I run a vibe-qc calculation?” Added docs/running.md: the virtualenv Python pattern, OMP_NUM_THREADS, tee / nohup / tmux for long jobs, and a common-errors table.

v0.4.3: the update workflow needs to be one command (+8h 11m)

Existing checkouts needed a four-step manual sequence to update: git fetch, checkout the tag, pull, re-run setup_native_deps.sh, pip install, verify. That is too many steps to remember, too many places to make a mistake, and too likely to leave a user on a stale build without realising it.

scripts/update.sh replaced all of it: one command, handles tag/branch/remote refs uniformly, refuses to run on a dirty tree, prints the release banner at the end so the user knows which version they are on. Paired with docs/updating.md covering the common failure modes: banner shows old version (re-install the Python package), missing library (pass --rebuild-native-deps).

v0.4.4 + v0.4.5: two Linux loader bugs, one after the other (+8h 38m, +8h 57m)

These two patches have the same root — the mismatch between how macOS and Linux resolve shared library dependencies — but they hit at different levels of the dependency tree, which is why they needed two separate fixes.

Side-by-side diagram showing DT_RUNPATH non-transitivity causing libxc and libcerfcpp not found, versus DT_RPATH transitivity fixing both
v0.4.4 added rpath entries for all five vendored deps in vibeqc_core.so, fixing libxc.so.15. v0.4.5 addressed a deeper issue: DT_RUNPATH (CMake’s default) is non-transitive — when libecpint loads its own deps, the loader does not inherit vibeqc_core.so‘s rpath. Switching to DT_RPATH via -Wl,--disable-new-dtags plus $ORIGIN in libecpint’s own build fixed libcerfcpp.so.3.

v0.4.4 fixed ImportError: libxc.so.15: cannot open shared object file. The v0.4 vendoring commit had added libxc, spglib, FFTW, and libecpint under third_party/<dep>/install/lib/ but only listed libint and libecpint in the compiled extension’s rpath. macOS’s DYLD_LIBRARY_PATH fallback search resolved the others silently. Linux’s glibc loader does not have that fallback. Fix: a foreach loop over all five vendored deps adds an rpath entry for each.

v0.4.5 fixed libcerfcpp.so.3 not found even after v0.4.4 — same symptom, different cause. After the extension loaded libecpint successfully (because libecpint was now in the rpath), libecpint itself tried to load its own transitive dependencies (libcerfcpp, pugixml) and failed. The reason: modern CMake emits DT_RUNPATH by default, which is non-transitive. The loader searches it only for the binary’s direct dependencies. When libecpint’s loader runs, it does not inherit the extension’s rpath. Old-style DT_RPATH is transitive. Two-prong fix: -Wl,--disable-new-dtags in cpp/CMakeLists.txt forces DT_RPATH for the extension, and -DCMAKE_INSTALL_RPATH='$ORIGIN' baked into the libecpint build lets it find its own siblings via $ORIGIN.

macOS never revealed either bug because @rpath on macOS is transitive by default and DYLD has the fallback. Every month of development on macOS was silently papering over a problem that any Linux user would hit in the first thirty seconds.

v0.4.6: when you split a change across two commits, write it down (+9h 15m)

The v0.4.5 hotfix was assembled under time pressure: cherry-pick the docs-CI banner fix, ship. What got left behind was the companion commit that adds the RELEASE_CODENAMES catalogue to vibeqc.banner and threads it into the runtime banner string. The result: import vibeqc; print_banner() on v0.4.5 read Release v0.4.5 with no codename. Nine new tests in tests/test_banner_codename.py shipped with v0.4.6 to pin the catalogue contract permanently: PEP 440 dev-suffix stripping, patch-to-minor fallback, three-step lookup. The banner now reads:

Release v0.4.6 "Schrödinger's Llama"  —  Quantum chemistry for molecules and solids

The fix itself was trivial. The lesson is about process: when a logical change spans two commits and you are mid-crisis, write down which commit goes with which patch. We did not, and a 30-minute omission cost a tagged release.

v0.4.7: the relative-path trap and the Furo TOC (+9h 41m)

Two more user-visible bugs caught by reading the live v0.4.6 docs site cold.

The quickstart said: save water.py, then run .venv/bin/python water.py. A user who ran git clone, then cd ~, then followed the snippet hit .venv/bin/python: No such file or directory. The venv is in the source tree. The quickstart assumed the reader was still in it. The fix: mkdir -p ~/vibeqc-runs/<project> plus explicit absolute paths throughout — ~/path/to/vibeqc/.venv/bin/python water.py. The activate-venv tip uses the same absolute form so it works from any working directory.

Three pages (quickstart, running, updating) were rendering a visible red Furo error block inline: ERROR: Adding a table of contents.... The Furo theme fires this as an in-page assertion when a {contents} directive is present but the right-sidebar already shows the same navigation. Fix: remove the redundant {contents} blocks. The assertion was always technically correct; the sidebar was already doing the job.

Two-panel diagram: left shows what broke across five categories, right shows why each bug was invisible during development
Left: the five categories of bugs across the eight patches, with patch counts. Right: the four structural reasons each category was invisible during development on macOS.

Three patterns worth naming

First-deploy bugs are almost entirely environment, not code. The glibc loader’s strictness, DT_RUNPATH vs DT_RPATH semantics, Arch packaging conventions, CMake 4 policy changes — none of these reproduce on the dev machine. They surface within minutes of the first non-dev user trying to install. The only defence is testing on Linux before declaring “ready to ship,” and having a patch-release workflow that is fast enough to respond.

Docs need to be followed, not just built. CI proved the docs were warning-free. It did not prove they worked. v0.4.2 and v0.4.7 are both “I followed my own quickstart and hit a wall” patches. The difference between building docs and following them is the author’s accumulated context — the venv path, the working directory, the sequence of steps that is obvious to someone who wrote the code and invisible to someone who just cloned it.

The patch-release workflow has to be in place before you ship. Seven patches in nine hours is only possible because scripts/update.sh, docs/release_process.md, and the GitLab CI auto-deploy-on-tag pipeline were all documented and operational before v0.4.0 tagged. Without that infrastructure, each patch would have been a half-day affair. The workflow is part of the release.

What comes next

H2 molecular crystal band structure and density of states, HF/STO-3G
Band structure and DOS for the H₂ molecular crystal from tutorial 12: bonding band flat near the Fermi level, antibonding band showing dispersion along Γ to X. One of the tutorials that shipped in the v0.4 window alongside the installer fixes.

With the “Schrödinger’s Llama” arc closed, the next milestone is v0.5.0 “Wilson’s Otter” — analytic forces, vibrational frequencies, and the periodic gradient machinery that unlocks phonons and elastic constants. The tutorial series grows with it: 25 tutorials at v0.4.7, with the Peierls dimerisation (tutorial 17), NEB reaction paths (tutorial 19), DFT functional comparison (tutorial 15), and the band structure and DOS (tutorial 12) all shipping in the v0.4 window.

The code is at vibe-qc.com, installation instructions in the installation guide, and scripts/update.sh handles everything for existing checkouts. MPL 2.0.

vibe-qc v0.4 “Schrödinger’s Llama”: from Day 3 to first public release

vibe-qc is an open-source quantum chemistry and solid-state code — C++17 backend, Python frontend via pybind11, ASE Calculator interface — being written in the open. Day 1 shipped the molecular stack and Ewald infrastructure. Day 2 closed v0.2.0. Day 3 set the public-release bar at v0.4.0 and shipped the first convergence aid. This post covers everything from that decision to v0.4.6 “Schrödinger’s Llama” — 150 commits on main, six patch tags on release, 131 tests growing to 909, and the gap between “compiles on the developer’s machine” and “actually installs on Manjaro” navigated in full.

Three-panel summary: periodic SCF method matrix all shipped, six-patch release cascade, test count 131 to 909
Left: every cell of the periodic SCF matrix — RHF, UHF, RKS, UKS × Γ-only and multi-k — shipped by v0.4. Centre: six hotfixes in 48 hours, four driven by Linux installer bugs. Right: test count from Day 1 through v0.4.6, zero regressions across all six patches.

Closing the method matrix

The v0.4 engineering goal was to close every combination of method, spin, and k-sampling for periodic Ewald SCF. Eight cells: RHF and UHF, RKS and UKS, each at Γ-only and multi-k. Day 2 had RHF at both k-samplings. The remaining six shipped across the v0.4 arc.

Periodic UHF (Phases 15a and 15b) opened open-shell Hartree-Fock on solids — magnetic systems, spin-polarised defect cells, antiferromagnetic unit cells. The unrestricted SCF minimises the energy separately in the $\alpha$ and $\beta$ spin channels, with the constraint that both share the same Coulomb field. Periodic RKS and UKS (Phases 15c-1 through 15c-3) required a new C++ kernel, build_xc_periodic_uks, to handle open-shell libxc evaluation on a LatticeMatrixSet density. The exchange-correlation energy for the open-shell case is

$$E_{\mathrm{xc}}[\rho^\alpha, \rho^\beta] = \int \varepsilon_{\mathrm{xc}}(\rho^\alpha(\mathbf{r}), \rho^\beta(\mathbf{r}),|\nabla\rho^\alpha|, |\nabla\rho^\beta|)\,d\mathbf{r}$$

evaluated on the periodic grid with Becke partitioning extended to image atoms. With that kernel in place, the SCF matrix was closed: any combination of method and spin runs at either Γ-only or multi-k with the same dispatcher entry point. The periodic DFT tutorial, the tight-cell Becke tutorial, and the periodic SCF convergence tutorial all exercise different corners of this matrix.

ECPs, dispersion, and convergence aids

Three capability tracks landed on top of the SCF matrix.

Phase 14 added effective-core potentials via libecpint 1.0.7, vendored with pugixml and libcerf. ECP support is wired through all four molecular drivers and validated against PySCF on Zn²⁺/LANL2DZ to micro-hartree accuracy. This opens the pob-* ECP basis sets for Rb through Lu — the heavy-element range that covers metal oxides, perovskites, and lanthanoid chemistry. The ECP integrals add a non-local potential term to the core Hamiltonian:

$$\hat{V}_{\mathrm{ECP}} = \hat{U}_{\mathrm{local}} + \sum_{\ell=0}^{L} \hat{U}_\ell \hat{P}_\ell$$

where $\hat{P}_\ell$ projects onto angular momentum $\ell$ and $\hat{U}_\ell$ is the semi-local correction. The integration hit one genuinely tricky bug: the libecpint API header documentation was actively misleading about whether the function appended /xml/ to the share-dir argument. It does. The source at api.cpp:73 was the authority.

Phase D1 added DFT-D3(BJ) dispersion correction, wired through run_job and the ASE Calculator. The dispersion tutorial shows the binding-curve comparison that makes the underbinding failure of uncorrected GGAs visible in a single plot.

The convergence track (C1a through C1c) shipped three tools. C1a is Saunders-Hillier level shifting (covered in the Day 3 post), which modifies the Fock matrix before diagonalisation as $\tilde{\mathbf{F}} = \mathbf{F} + b(\mathbf{S} – \frac{1}{2}\mathbf{S}\mathbf{D}\mathbf{S})$, raising virtual eigenvalues by $b$ without changing the converged fixed point. C1b added Fermi-Dirac smearing with fractional occupations. C1c shipped a quadratic SCF fallback — a Newton step in MO space,

$$\Delta\mathbf{\kappa} = -\mathbf{H}^{-1}_{\mathrm{diag}}\,\mathbf{g}$$

where $\mathbf{g}$ is the orbital-rotation gradient and $\mathbf{H}_{\mathrm{diag}}$ is the diagonal orbital Hessian, with a trust-region cap to prevent overlarge steps. Default is off; opts.quadratic_fallback_iter = N activates after iteration $N$. The SCF convergence user guide covers all three tools in order of when to reach for them.

The bug that tutorials found

The most instructive quality event of the v0.4 arc was a 12-item bug report that came from the process of writing tutorials. The pattern: reading the API carefully enough to explain it surfaces issues that running the API does not.

The sharpest example was level-shift silently dropped in the dispatcher. Engineering had shipped C1a, added level_shift to SCFOptions, and the tests passed. But _copy_options_to_rhf and _copy_options_to_scf copied nine fields by hand and had not been updated. Setting level_shift=0.4 in a tutorial script produced no effect on the SCF trace. The fix was one line in two functions. The regression test now asserts that a dispatcher round-trip preserves all SCFOptions fields — the same class of bug cannot recur silently.

The other significant find: the basis library was not inside the wheel. basis_library/ lived at repo root with basis_library/basis/ gitignored, populated by a shell script at build time. A pip wheel does not run shell scripts. The fix was moving the whole tree to python/vibeqc/basis_library/. Ninety standard .g94 basis sets plus three solid-state pob-* sets now ship inside the wheel. No environment variable required.

Schrödinger’s Llama: the patch-release cascade

Schrödinger's Llama comic: a worried llama in a wooden box labelled macOS OK, Linux question mark, Arch triple question mark, with a radioactive pip install vibe-qc flask and a thought bubble asking whether v0.4.5 is fixed on Manjaro
The release is simultaneously “working” and “broken on Linux” — until a user opens the box. Six patch releases in 48 hours, four of them Linux installer bugs that development on macOS had never revealed.

Tagging v0.4.0 followed docs/release_process.md: bump pyproject.toml, write the changelog, tag from main, fast-forward release. The banner started reading Release v0.4.0. Then the first user on a different machine tried to install it.

v0.4.1 hit three independent Arch/CMake failures. libint’s Eigen detection used ls a b, which returns nonzero if any argument is missing — on Arch, Eigen 5 installs only at /usr/include/eigen3/Eigen/Core, not the top-level path that the detection script expected. libint’s FindBoost invocation triggered fatal CMake 4.x deprecation warnings. And all three vendored dependencies declared cmake_minimum_required(VERSION <3.5), which CMake 4.x rejects outright. None of these were vibe-qc bugs. They were upstream code aging into newer toolchains, absorbed because the install path is part of the public surface.

v0.4.4 and v0.4.5 were the genuinely interesting failures. After vendoring five native libraries, the rpath only mentioned two of them. macOS found the others via DYLD‘s fallback search. Linux does not.

Side-by-side dependency tree showing DT_RUNPATH non-transitivity causing libxc not found, versus DT_RPATH transitivity fixing it
Left: under DT_RUNPATH (the modern GNU ld default), the search path in vibeqc_core.so does not propagate to libecpint’s own dependency search — libcerfcpp.so.3 is not found. Right: -Wl,--disable-new-dtags forces DT_RPATH, which is transitive. macOS’s @rpath is transitive by default, which is why this was Linux-only.

v0.4.4 added a foreach loop over all five vendored deps — problem solved for direct dependencies. Then v0.4.5 hit libcerfcpp.so.3 not found, even though the file existed on disk. The root cause was DT_RPATH vs DT_RUNPATH in ELF binaries. Modern GNU ld defaults to emitting DT_RUNPATH, which is non-transitive: the rpath in vibeqc_core.so applies only to its direct dependencies. When libecpint loads its own deps, the loader does not inherit that rpath. Two-prong fix: -Wl,--disable-new-dtags on the Linux link of vibeqc_core.so (forcing DT_RPATH), plus -DCMAKE_INSTALL_RPATH='$ORIGIN' baked into libecpint at build time. The ELF loader’s transitivity semantics are not something most scientific software developers need to know about — until they do.

v0.4.6 closed the codename system. RELEASE_CODENAMES lives in python/vibeqc/banner.py. The docs CI container does not install vibe-qc (no C++ stack in python:3.13-slim), so importlib.metadata.version("vibe-qc") fell through to a placeholder. The fix was loading banner.py as a standalone module via importlib.util, bypassing vibeqc/__init__.py entirely. Docs and runtime now read the same dictionary without either depending on the other’s import chain. Every SCF log from v0.4.6 onwards carries Release v0.4.6 "Schrödinger's Llama" — unambiguous build provenance in every pasted traceback.

Documentation: 19 tutorials to 25

Six new tutorials shipped: natural orbitals and the idempotency diagnostic, projected density of states, periodic Bloch orbital cubes, tight-cell DFT with the periodic Becke partition, periodic SCF convergence, and symmetry-aware lattice integral storage. All 25 now carry resource callouts — peak RSS and wall-clock on a reference machine, generated by a dev script and pasted in. The quickstart was restructured from a 2000-line everything-page into four focused documents.

The CI/CD side replaced a silent nightly cron with a GitLab pipeline: sphinx-build -W --keep-going (warnings are errors), rsync to the production webroot, triggered only on the release branch. One gotcha worth documenting: GitLab File-type variables silently corrupt multi-line OpenSSH private keys. Store the key base64-encoded in a regular Variable and decode it in the YAML with base64 -d.

Going public

GitLab project visibility flipped to Public at v0.4.0. Anonymous clone lands on the release branch by default — the latest tagged release, not the dev main. The project has a communication address (mpei@vibe-qc.com), a PGP key for security reports (fingerprint CC6D 30BB DF96 F694 C615 FBDE 4CD5 65CF 26B1 E7E5, public key at vibe-qc.com/_static/pgp/mpei.asc), a SECURITY.md, and a CONTRIBUTING.md decision tree.

The number that matters

909 tests pass on a fresh install. The verification command is in installation.md. Zero regressions across six patch releases. The feature work — closing the periodic SCF matrix, ECPs, D3(BJ), three convergence tools — was the straightforward part. The work that made it real for someone on a different machine was four patches driven by Linux loader semantics, one by CMake version policy, and one by a codename dictionary that two separate CI environments needed to read from the same source. All six are in the changelog. All six are fixed permanently.

The bug arc continues: Eight releases in nine hours: the v0.4.0 to v0.4.7 story.

The code is at vibe-qc.com. MPL 2.0. Clone it, run the tests, read the 25 tutorials.

vibe-qc day 3: setting the public-release bar and the first convergence fix

Left: SCF energy vs iteration showing oscillation without level shift and smooth convergence with level shift b=0.3. Right: vibe-qc v0.4.0 public release roadmap.

vibe-qc is an open-source quantum chemistry and solid-state code — C++17 backend, Python frontend via pybind11, ASE Calculator interface — being written in the open one session at a time. Day 1 shipped the molecular stack and the first Ewald infrastructure. Day 2 closed out v0.2.0 — end-to-end 3D Ewald dispatch, bulk benchmarks, multi-k DIIS, periodic Becke partition — and opened v0.2.5 with SYM3a orbit-reduced lattice-integral storage. Day 3 was narrower by design: a release-target decision and the first piece of the next milestone. 746 tests in, 754 out.

Setting the public-release bar at v0.4.0

The first question of the day was where to draw the public-release line. Three candidates.

v0.2.0 is technically coherent — “3D periodic bulk HF/DFT with quantitative Ewald Coulomb” is a complete story, and we are 95% of the way there (the remaining piece is tightening the CRYSTAL cross-check witness bounds on LiH, NaCl, MgO, and Si). But shipping there would hand a user a code that oscillates the moment they try a metal or a narrow-gap oxide. The DIIS accelerator that works beautifully on insulators stalls on systems with near-degenerate occupied and virtual orbitals near the Fermi level. That is not a missing feature — it is a known failure mode with known remedies, and shipping before those remedies exist would damage the project’s credibility faster than it builds it.

v0.3.0 adds visualisation polish: properties, cube files, better band plots. Useful, but it does not fix the convergence problem. A code that looks nicer while still oscillating on MgO is not a more useful code.

v0.4.0 is where vibe-qc becomes “I can do real solid-state chemistry with this.” The convergence-tooling track (C1: level shifting, Fermi-Dirac smearing, second-order SCF) handles the metals and tight-gap systems. Phase 14 adds effective-core potentials via libecpint, opening the pob-* ECP basis sets for Rb through Lu — metal oxides, perovskites, lanthanoids. Phase 15 adds periodic UHF/UKS for magnetic systems and spin-polarised transition-metal compounds. That combination covers the tutorial cases that currently bounce off the wall, and it matches the capability level where researchers in solid-state chemistry will actually trust the results enough to publish. That is the bar. The full roadmap is on vibe-qc.com.

Left: SCF energy vs iteration showing oscillation without level shift and smooth convergence with level shift b=0.3. Right: vibe-qc v0.4.0 public release roadmap.
Left: the level-shift effect on a tight-gap periodic cell — unshifted SCF oscillates; $b = 0.3$ converges in 9 iterations to the same energy. Right: the v0.4.0 public-release roadmap, with today’s C1a already green.

Phase C1a: Saunders-Hillier level shift in periodic Ewald SCF

The first piece of the v0.4.0 convergence track shipped today. The Saunders-Hillier level shift (Saunders and Hillier, 1973; also Goedecker, 1996) modifies the Fock matrix before diagonalisation:

$$\tilde{\mathbf{F}} = \mathbf{F} + b\left(\mathbf{S} – \tfrac{1}{2}\mathbf{S}\mathbf{D}\mathbf{S}\right)$$

where $b$ is the shift parameter and $\mathbf{D}$ is the density matrix. The effect is to raise virtual MO eigenvalues by $b$ while leaving the occupied block untouched. This widens the effective HOMO-LUMO gap seen by the diagonaliser each iteration, suppressing the near-degenerate orbital swaps that cause oscillation in small-gap systems. The SCF fixed point is unchanged — at convergence, $\mathbf{F}\mathbf{D} = \mathbf{D}\mathbf{F}$ (in the $\mathbf{S}$-metric), so the shift term vanishes and the converged density is the same regardless of $b$.

That last point is the correctness contract that the test suite validates directly. A tight H₂ chain at $a = 8$ bohr (Γ point) and at $a = 10$ bohr with a $[2,2,2]$ k-mesh: $b = 0.3$ reproduces the $b = 0$ total energy to within $10^{-9}$ Ha. The reported mo_energies in the result object are the un-shifted physical orbital energies — the final self-consistency pass at convergence omits the shift, so the orbital spectrum is independent of the $b$ value chosen during iteration. Choosing $b$ too large just slows convergence; choosing it too small provides no benefit on the difficult cases. A value between 0.2 and 0.5 Ha covers most practical situations.

The implementation is wired into both the Γ-Ewald (run_rhf_periodic_gamma_ewald3d) and multi-k Ewald (run_rhf_periodic_multi_k_ewald3d) drivers via a new level_shift field on PeriodicRHFOptions, PeriodicSCFOptions, and PeriodicKSOptions. The default is 0.0 — the field is purely opt-in, and existing callers see bit-for-bit identical results. The SCF convergence section of the user guide covers when and how to use it. Eight new tests: field exposure on all three option structs, Γ-Ewald and multi-k inertness contracts, MO eigenvalues physical at convergence, default behaviour preserved. The molecular RHF/UHF/RKS/UKS level-shift wiring (C1a-2) needs a small parallel-code change in four C++ drivers and is deferred to a follow-up commit.

A segfault, triaged but not fixed

Honest engineering means mentioning this. A crash surfaced in the molecular gradient code: compute_gradient(Molecule, BasisSet, RHFResult)overlap_gradient_contribution → null function pointer inside __kmp_invoke_microtask. Multiple OpenMP worker threads hitting pc = 0x0 is a strong signal of either a libint2 Engine dispatch table that was not initialized when freshly-copied engines raced into .compute() from worker threads, or an out-of-bounds access into the engine pool when nested or dynamic OpenMP ran more threads than the pool was sized for.

The existing 32 gradient tests all pass cleanly, so the crash triggers on a specific input pattern — most likely a basis set with shells exceeding the prototype’s max_l, or a thread-count change between pool construction and the parallel region. The crash binary’s offsets shifted after rebuild, making direct symbolisation unhelpful. This is deferred: waiting on a reproducer that pins the exact shell configuration. It will be fixed before v0.4.0 tags, because the gradient code is load-bearing for geometry optimisation and eventually periodic forces.

What is next

Two immediate priorities. First, C1b — Fermi-Dirac smearing with fractional occupations and an electronic-entropy contribution to the free energy. This is the single biggest credibility-unblocker for v0.4.0: it covers metals, oxide surfaces, and defect cells with mid-gap states that DIIS plus level shifting still cannot handle. Second, the segfault reproducer — once we have an input that reliably triggers the crash, the fix should be straightforward.

SYM3b (kernel-level compute reduction — the $|G|$-fold wall-clock speedup from space-group symmetry, not just the memory saving from SYM3a) remains open on the v0.2.5 track and will slot in once C1b is done. The CRYSTAL cross-check pass that formally closes v0.2.0 is also still outstanding.

The full v0.4 story continues in the next post: vibe-qc v0.4 “Schrödinger’s Llama”.

The code is at vibe-qc.com. MPL 2.0.

When the tutorial is the product: how vibe-qc treats documentation as a first-class deliverable

Split view: theory section of a vibe-qc tutorial on the left, running Python code on the right

vibe-qc is an open-source quantum chemistry and solid-state code: a C++17 numerical backend, a Python frontend via pybind11, and an ASE Calculator interface that connects it to the wider atomistic simulation ecosystem. The code is being written in the open, one session at a time, with the cyclic cluster model — a method for treating solid-state defects at post-HF accuracy — as its eventual destination. But the thing that distinguishes vibe-qc from other projects at this stage is not its roadmap length. It is that the tutorials are treated as first-class deliverables, written concurrently with the code they document, and held to the same quality gate as the implementation itself.

This post is about why that choice was made and what it looks like in practice.

What a tutorial is, in this project

Every tutorial in vibe-qc follows the same structure: a working Python script, representative output, an embedded figure where the calculation produces something visual, a theory section (roughly half a page, three to six equations), a references block, and a pointer to the next tutorial in the sequence. None of these sections are optional. A tutorial without a theory section is a recipe. A tutorial without a working script is an essay. The combination is what makes it an educational artefact rather than documentation.

The theory sections are typeset with explicit equations. Since peintinger.com does not currently run a MathJax or KaTeX plugin, the equations in the docs themselves are rendered by MathJax inside the Sphinx build at vibe-qc.com. For example, the molecular Hartree-Fock tutorial works through the Roothaan equations in matrix form, $\mathbf{F}\mathbf{C} = \mathbf{S}\mathbf{C}\boldsymbol{\varepsilon}$, where $\mathbf{F}$ is the Fock matrix, $\mathbf{C}$ the MO coefficient matrix, $\mathbf{S}$ the overlap matrix, and $\boldsymbol{\varepsilon}$ the diagonal matrix of orbital energies:

Split view: theory section of a vibe-qc tutorial on the left, running Python code on the right
The dual character of a vibe-qc tutorial: theory on the left (Bloch theorem, dispersion relation, DOS formula, foundation references), running code on the right. Both sections are required; neither is optional.

The citation discipline matters as much as the equations. Foundation-era papers — Roothaan 1951, Bloch 1928, Becke 1988, Pulay 1980 — are cited with full journal-volume-page triples drawn directly from the original. Post-2010 results that are harder to verify from memory carry a [Ref: verify] marker in the draft and are only promoted to full citation after the volume and page have been confirmed against the journal. This is not the standard practice for informal software documentation. It is standard practice for a journal article, and the tutorials are being written to that standard deliberately.

The parity matrix as a dual roadmap

The roadmap page carries two distinct structures. The first is an engineering milestone list tracking Ewald phases, symmetry sub-phases, and versioned capabilities. The second is a tutorial parity matrix that pins vibe-qc against ORCA’s molecular tutorial set and CRYSTAL’s solid-state tutorial set — the two reference programs that cover the two halves of what vibe-qc is trying to do.

Table showing vibe-qc tutorial parity against ORCA and CRYSTAL reference programs
The tutorial parity matrix from vibe-qc.com/roadmap.html. Green means a tutorial covering this capability is shipped and cross-checked against the reference program. Yellow means partial. Red means the capability is on the roadmap but the tutorial does not yet exist.

The value of the parity matrix is that it answers two different questions from a single table. For a user evaluating whether vibe-qc can support their work, it gives a one-page answer: “can I reproduce the ORCA tutorial on geometry optimisation with this code?” For a developer deciding what to build next, it gives a prioritised queue: the red cells in the CRYSTAL column are where the periodic capability is missing and the tutorials cannot yet be written. The engineering roadmap and the documentation roadmap are the same document.

The framing is deliberately not competitive. The parity matrix does not claim vibe-qc exceeds ORCA or CRYSTAL at anything. It says, plainly, where vibe-qc is today relative to what those programs can teach, and what it would take to close the gap. That honesty is part of the project’s credibility.

Fourteen tutorials, all shipped today

As of the current tutorial index, fourteen tutorials are live. They run from molecular Hartree-Fock — the entry point, which teaches the Roothaan equations, basis set input, and how to read an SCF trace — through open-shell methods, geometry optimisation, vibrational frequencies, thermodynamics at finite temperature, orbital visualisation, basis set convergence, and dispersion corrections, and then into the periodic side: periodic HF on a 1D hydrogen chain, periodic KS-DFT, Madelung constants via Ewald summation, and band structure and density of states.

The band structure tutorial is worth looking at as the canonical example of what the format can carry. It delivers a full Γ-to-X k-path calculation on the hydrogen chain, plots bands and DOS in a single combined figure, explains the Bloch theorem, $\psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}}u_{n\mathbf{k}}(\mathbf{r})$, and the origin of band dispersion in terms of tight-binding bonding and antibonding combinations, cites Bloch 1928 and Monkhorst-Pack 1976, and ends with a pointer forward to the Peierls distortion. The code block is short enough to read in under a minute. The theory section is long enough to understand what the code is computing and why the result looks the way it does.

That rhythm — code, output, figure, theory, references, next — is the same across all fourteen tutorials. A student who works through them in order gets a coherent course in molecular and periodic quantum chemistry. A researcher who arrives at a specific tutorial looking for a recipe gets that too, because the code blocks are self-contained. The two use cases are not in tension; they are served by the same document structure.

The no-code-required backlog

One discipline the roadmap makes explicit is the distinction between tutorials that require new code and tutorials that require only documentation work on existing capability. Several entries in the tutorial backlog fall into the second category: the capability is already in the code, the API is stable, and the only work is writing the tutorial itself — the explanation, the example script, the figure, the theory section, the references.

This is worth naming because it is easy to defer. The code works; the tutorial can wait. What the project has committed to instead is treating those deferred tutorials as a first-class backlog item, tracked in the parity matrix, with the same visibility as unimplemented features. A shipped capability without a tutorial is a capability that most users will never find.

The user guide plays a complementary role here. It covers the reference-level detail — how to specify k-point meshes, how to configure Ewald parameters, how to read the memory budget report — that would clutter a tutorial but that a user needs once they move past the introductory examples. The split between tutorial and user guide is deliberate: tutorials teach phenomena and procedures, the user guide documents options and behaviour.

Plots as the lesson, not decoration

Every tutorial that produces a figure ships that figure embedded in the page. The figures are not screenshots pasted in after the fact. They are generated by scripts in the repository, regenerable by any user who runs the tutorial code, and updated whenever the underlying calculation changes.

This matters because the figures are doing pedagogical work that the code blocks cannot. A code block teaches syntax and API. A band structure plot teaches what band dispersion looks like, where the gap sits, how the DOS integrates the band structure into a density, and what it means for a state at the zone boundary to be antibonding. The band structure tutorial includes both panels — bands and DOS in a single combined figure — because showing them side by side is how a reader learns to interpret both simultaneously.

The same logic applies to the dispersion tutorial, which shows a binding curve with and without D3-BJ correction, making visible the underbinding failure of uncorrected GGAs on weakly-bound systems. The figure is the argument. The code block shows how to reproduce it.

The doc build as a gate

The documentation is built with sphinx-build -W --keep-going: warnings are errors, and the build does not finish if any cross-reference is broken or any code block fails to parse. This is the same discipline applied to the test suite — 746 passing tests as of the most recent commit — extended to the documentation layer. A tutorial that links to a non-existent API page, or a code block that uses a function name that was renamed, fails the build and blocks the merge.

This is not the default posture for scientific software projects. The default is to build docs as a separate step, treat failures as non-blocking, and accumulate rot over time. The cost of running docs under the same gate as code is low — a few minutes added to CI. The benefit is that the tutorials stay synchronized with the implementation without anyone having to remember to update them.

Where this fits in the larger project

vibe-qc is aiming at the cyclic cluster model: a method that treats a finite cluster of atoms embedded in a periodic Madelung field, enabling post-HF correlation calculations on solid-state defects at a cost that scales with cluster size rather than unit cell size. The cyclic cluster model is v2.0 on the roadmap. Every tutorial from the molecular HF entry point through the periodic band structure is load-bearing preparation for the moment when a researcher can open a CIF file, define a vacancy cluster, run CCSD(T) on it, and understand what the result means.

The tutorials exist not to fill a documentation requirement but because the project’s audience includes people who need to understand what they are computing, not just how to invoke the code. That is a different audience from “users who already know periodic HF/DFT and just need the API reference.” The user guide and API reference serve that second audience. The tutorials serve everyone, and they are the harder document to write well.

The code lives at vibe-qc.com. MPL 2.0.

vibe-qc day 2: closing out Ewald and laying the symmetry foundation

vibe-qc day 2: closing out Ewald and laying the symmetry foundation

vibe-qc is an open-source quantum chemistry and solid-state code with a C++17 numerical backend, a Python frontend via pybind11, and an ASE Calculator interface for geometry optimization and structure I/O. It is being written in the open, one coding session at a time, with Claude Opus as the implementation engine and me steering at the architecture level. Day 1 shipped the molecular HF/DFT/MP2 stack (v0.1.0, April 18) and then pushed into periodic territory: one-electron lattice integrals, Gamma-only and multi-k RHF, multi-k KS-DFT, and the first Ewald infrastructure through Madelung-cancellation helpers. Day 2 closed out the v0.2.0 Ewald milestone and planted the first stake for v0.2.5 space-group symmetry. Four deliverables, 26 new tests (720 to 746 total), Sphinx docs building clean.

Phase 12e-c-4: end-to-end EWALD_3D dispatch

The previous session had the Ewald machinery in place but not yet wired to the user-facing entry points. Today that changed. Both run_rhf_periodic_scf(...) and run_rhf_periodic_gamma_scf(...) now branch on options.lattice_opts.coulomb_method: pass DIRECT_TRUNCATED and you get the existing C++ driver unchanged; pass EWALD_3D and you get the new Python driver, with nuclear_repulsion_per_cell automatically rerouted through Ewald summation as well. That last part matters: an Ewald electronic energy paired with a direct-truncation nuclear repulsion would be self-consistent in neither the long-range cancellation nor the ω-dependence, and the resulting total energy would be meaningless. Electronic and nuclear sides now always agree on which summation scheme is in use.

The benchmark suite that shipped alongside covers H₂ crystal, LiH rocksalt, MgO rocksalt, and Ne FCC — 11 tests in total, checking convergence, $\omega$-invariance (energy independence of the Ewald splitting parameter), an equation-of-state scan, and the equivalence of a [1,1,1] multi-k mesh with a Gamma-only calculation. All pass.

One finding worth documenting: atom positions relative to the FFT Poisson grid origin are not numerically irrelevant. Atoms sitting at box corners, rather than centered in the cell, inflate the $\omega$-invariance residual by roughly two orders of magnitude. The physics is correct either way, but the numerical conditioning is not. This is now flagged in the docs and the benchmark suite uses centered geometries as the reference.

Phase 12e-c-4c-iii-c: multi-k Pulay DIIS

The multi-k Ewald driver shipped without DIIS. For loose cells and [1,1,1] meshes that was acceptable — pure density damping converged, slowly. For tight cells with denser k-meshes it was not. An H₂ chain at a = 10 bohr with a [2,2,2] mesh simply plateaued under damping and never converged.

The fix is Pulay DIIS extended to the Brillouin zone. The commutator error vector at each k-point is:

$$\mathbf{e}(\mathbf{k}) = \mathbf{F}(\mathbf{k})\mathbf{D}(\mathbf{k})\mathbf{S}(\mathbf{k}) – \mathbf{S}(\mathbf{k})\mathbf{D}(\mathbf{k})\mathbf{F}(\mathbf{k})$$

and the Pulay B matrix couples iterations i and j with a k-weighted inner product:

$$B_{ij} = \sum_{\mathbf{k}} w_{\mathbf{k}}\, \mathrm{Re}\,\mathrm{tr}\!\left[\mathbf{e}_i(\mathbf{k})^\dagger\, \mathbf{e}_j(\mathbf{k})\right]$$

The k-weights are the standard Monkhorst-Pack weights that already appear in the energy and density assembly, so no new parameters appear. The DIIS coefficients solve the same linear system as in the molecular case; the only change is that the scalar error measure is now a sum over the zone rather than a single matrix norm.

The numbers are clean. H₂ chain, [1,1,1]: 13 iterations under damping, 7 with DIIS. H₂ chain, [2,2,2]: previously non-converging under damping, now 4 iterations with DIIS. This is the same qualitative behavior the molecular DIIS delivered on H₂O back on day 1 (52 iterations to 9), just extended to k-space.

Phase 12f: periodic Becke partition for tight-cell DFT

The molecular Becke fuzzy-cell partition (J. Chem. Phys. 88, 2547, 1988) assigns each grid point a weight that partitions unity across atomic cells using a smooth step function of interatomic distances. In a periodic system with a tight unit cell, image atoms from neighboring cells fall within the smoothing radius of the step function, and ignoring them breaks the partition: the weights no longer sum to the cell volume.

The fix is to extend the partition denominator over home atoms plus image atoms within a user-specified radius. The implementation is wired into run_rks_periodic via two new options: PeriodicKSOptions.use_periodic_becke (boolean, default False) and becke_image_radius_bohr (float). Default behavior is unchanged from v0.1.0; the periodic extension is an explicit opt-in for users who know they are working in a tight cell.

The empirical witness is a 5-bohr cubic H₂ cell. Without periodic Becke, the molecular partition produces a total grid weight of 18,526 — wildly wrong for a cell with volume 125 bohr³. With periodic Becke, the total grid weight comes out as 126, matching V_cell = 125 bohr³ to within the grid quadrature error. The discrepancy between the two modes (18,526 vs 126) is not a subtle numerical issue; it is a categorical correctness failure of the molecular partition in tight periodic geometry. Any DFT calculation on a tight crystal without this correction is integrating a density that does not integrate to the correct electron count per cell.

Phase SYM3a: orbit-reduced storage for lattice integrals

The v0.2.5 symmetry milestone starts here. The one-electron lattice integrals $S(\mathbf{g})$, $T(\mathbf{g})$, $V(\mathbf{g})$ are indexed by a cell offset vector $\mathbf{g} = (g_1, g_2, g_3)$. For a system with space-group symmetry, many (g, atom-pair) combinations are related by point-group operations and carry the same numerical content up to a known rotation. SYM3a exploits this to compress storage.

The machinery builds on SYM2c, which already identifies atom-pair and cell-index orbits for arbitrary (including non-origin-fixed) structures like NaCl. SYM3a takes those orbits and stores one representative sub-block of shape $(n_{\mathrm{AO},a},\, n_{\mathrm{AO},b})$ per orbit, rather than a full $(n_{\mathrm{bf}} \times n_{\mathrm{bf}})$ block at every cell triple. A round-trip compress/reconstruct test on simple cubic Pm-3m He at STO-3G, with a 10 bohr cutoff, partitions 33 cell triples into 5 orbits and achieves 6.6× memory reduction. Reconstruction is exact to 10⁻¹⁶.

What SYM3a does not yet do is reduce compute: the integral kernels still evaluate every shell-pair times cell-triple combination before handing off to the compressor. That is SYM3b, the next item on the list. SYM3b is where the |G|-fold reduction hits wall-clock time rather than just peak memory, because the kernel skips shell-pair/cell-triple combinations that are equivalent under the point group instead of evaluating and then discarding them. The storage handle SYM3a provides is what SYM3b will plug into.

Where v0.2.0 stands

The v0.2.0 contract is “ω-invariant 3D periodic SCF that reduces to molecular HF in the loose-box limit and integrates the unit-cell volume correctly in tight cells.” As of today that contract is end-to-end testable with two option flags:

options.lattice_opts.coulomb_method = EWALD_3D
options.ks_opts.use_periodic_becke = True

Internal-consistency benchmarks are in place across four crystal systems. The remaining gating item for the v0.2.0 tag is a cross-check pass against published CRYSTAL reference energies on LiH, NaCl, MgO, and Si at published geometries. That is a validation exercise against an external reference rather than new implementation work, and it is next on the list.

What is next

Immediately: SYM3b, the kernel-level compute reduction that makes the |G|-fold saving from space-group symmetry show up in wall-clock time. The substrate — Wigner D-matrices (SYM1), AO permutation matrices and orbit identification (SYM2a/b), atom-pair-resolved orbits for non-origin-fixed structures (SYM2c), and now the compressed storage layer (SYM3a) — is all in place. SYM3b is the payoff.

Further out, the headline feature of the entire project remains the cyclic cluster model. The roadmap runs from v1.0 (feature-complete molecular and periodic HF/DFT/MP2) through v2.0 (HF-CCM in 3D) and on through a sequence of independently publishable steps: MP2-CCM, local-MP2-CCM, CCSD-CCM, and eventually projection-based embedding where a CCM-correlated central cluster sits inside a DFT-treated periodic environment. The capstone is multi-scale defect chemistry at chemical accuracy on systems with hundreds of atoms in the cluster. Every piece shipped in the first two days of this project is load-bearing for that goal.

The code is at vibe-qc.com. Licensed MPL 2.0.

The giants vibe-qc stands on: a guided tour of the stack

In my previous post I wrote about vibe-coding a quantum-chemical program in one day. The code lives at vibe-qc.com. That post closed with a compact acknowledgment that none of it would have been possible without decades of prior work by a lot of people. This post is where I pay that debt properly: every library we used, what it does, why it matters, and how to cite it.

If you are building anything similar, treat this as a reading list. If you are just curious what sits underneath a modern quantum chemistry code, keep scrolling. I have organized things by role rather than alphabetically, because the interesting story is the layered structure of the stack.

The numerical core

libint

This is the beating heart of vibe-qc. Every one- and two-electron Gaussian integral, and every analytic derivative of those integrals, is computed by libint. It is a C++ library whose source is partly handwritten and partly emitted by its own code generator (which is itself a substantial piece of C++). We build version 2.13.1 from source with max_am=5 (angular momentum up to g-functions) and first-order derivatives enabled. The source build takes about 10 minutes on an M1, and you only do it once per checkout.

Without libint, implementing two-electron repulsion integrals with the right recurrence relations would consume weeks on its own. With libint, it is a function call.

Valeev, E. F. Libint: A library for the evaluation of molecular integrals of many-body operators over Gaussian functions, Version 2.13.1, 2024. https://libint.valeyev.net/

libxc

libxc is the library that turns vibe-qc’s DFT layer from “LDA only” into a proper general-purpose DFT code with 500+ functionals available via aliases like "PBE", "B3LYP", or "BLYP". Under the hood it exposes a uniform C API for evaluating exchange-correlation energy densities and their derivatives with respect to density and density gradient. Hybrid functionals report their exact-exchange fraction through xc_hyb_exx_coef, which is exactly what a SCF driver needs to know.

The “right” way to expose functionals in a new DFT code is to wrap libxc, not to reimplement functionals yourself. We did the right thing here.

Lehtola, S.; Steigemann, C.; Oliveira, M. J. T.; Marques, M. A. L. Recent developments in libxc, a comprehensive library of functionals for density functional theory. SoftwareX 7, 1-5 (2018). https://doi.org/10.1016/j.softx.2017.11.002

Marques, M. A. L.; Oliveira, M. J. T.; Burnus, T. Libxc: a library of exchange and correlation functionals for density functional theory. Comput. Phys. Commun. 183, 2272-2281 (2012). https://doi.org/10.1016/j.cpc.2012.05.007

Eigen

All dense linear algebra in vibe-qc goes through Eigen: matrix products, eigendecompositions (for the Fock diagonalization), linear solves (for the DIIS coefficient equation), and the energy-weighted density matrix construction in the gradient code. Eigen is header-only, which keeps the build simple, and the expression-template design means the code reads like math.

Guennebaud, G.; Jacob, B.; et al. Eigen v3. 2010. https://eigen.tuxfamily.org

OpenMP

Included in the library list because the infrastructure is wired, although the hot loops in vibe-qc are currently serial. Adding OpenMP to the Fock build and the ERI-gradient loops is on the day-two todo list. OpenMP is the de-facto standard for shared-memory parallelism in scientific C++ and it is what we will use.

Bindings and build

pybind11

pybind11 is what makes vibe-qc feel like a Python library even though all the numerical work is in C++. The bindings expose the C++ classes (Molecule, BasisSet, RHFResult, etc.) to Python with zero-copy NumPy interop for the heavy arrays, and the native work releases the GIL via py::gil_scoped_release, so the Python frontend is not a scaling bottleneck.

Jakob, W.; Rhinelander, J.; Moldovan, D. pybind11: Seamless operability between C++11 and Python. 2017. https://github.com/pybind/pybind11

scikit-build-core

The modern Python build backend that lets pip install -e . drive a CMake build under the hood. This is the glue that makes vibe-qc feel like a normal Python package even though it compiles ~3500 lines of C++ and links against three native libraries.

CMake and Ninja

CMake generates the build, Ninja runs it. Ninja is considerably faster than make for projects with thousands of small files (libint’s code generator emits a lot of them).

Boost and GMP

These are libint’s build-time dependencies, not vibe-qc’s. Boost provides header-only utilities (MPL, Type Traits, Preprocessor) that libint’s code generator relies on. GMP provides arbitrary-precision rational arithmetic, which the generator needs to emit exact coefficients in the integral recurrence relations. You install them, libint builds, and you never think about them again.

Granlund, T.; the GMP development team. GNU MP: The GNU Multiple Precision Arithmetic Library, Edition 6.3.0, 2023. https://gmplib.org

Python runtime

NumPy

Every matrix and vector that crosses the C++/Python boundary becomes a NumPy array, usually as a zero-copy view of the underlying Eigen data. NumPy’s array protocol is the reason pybind11 can do that zero-copy handoff cleanly.

Harris, C. R. et al. Array programming with NumPy. Nature 585, 357-362 (2020). https://doi.org/10.1038/s41586-020-2649-2

ASE (Atomic Simulation Environment)

ASE deserves its own paragraph because it is the decision that turned vibe-qc from “a Python library with an SCF loop” into “a program chemists can actually use.” By subclassing ase.calculators.calculator.Calculator, vibe-qc inherits:

  • Structure file readers for XYZ, CIF, PDB, POSCAR, Gaussian inputs, and basically every format chemists touch.
  • Geometry optimizers (BFGS, FIRE, LBFGS, conjugate gradient, etc.) that just work once forces are implemented.
  • Consistent unit handling via ase.units.
  • Visualization through ase.visualize.
  • An interface that every other major code in the ecosystem (VASP, Quantum ESPRESSO, GPAW, etc.) also implements, so users who already know ASE know vibe-qc.

If I were starting another scientific code tomorrow, I would wire it into ASE on day one, same as we did here.

Larsen, A. H. et al. The atomic simulation environment, a Python library for working with atoms. J. Phys.: Condens. Matter 29, 273002 (2017). https://doi.org/10.1088/1361-648X/aa680e

Reference and validation

PySCF

PySCF is the cross-check oracle for everything vibe-qc does. Every RHF, UHF, and RKS energy, every MO coefficient, every gradient component in the 131-test suite is compared against a PySCF calculation on the same system with the same basis set. HF agreement is at machine precision (around 1e-14 Ha), DFT at grid accuracy (down to 5e-11 Ha for LDA, 1.3e-7 Ha for B3LYP on the default medium grid).

If you are writing a new quantum chemistry code in 2026, you should be validating against PySCF. It is open, well-tested, and actively maintained. There is no excuse for shipping numbers that have not been cross-checked.

Sun, Q. et al. Recent developments in the PySCF program package. J. Chem. Phys. 153, 024109 (2020). https://doi.org/10.1063/5.0006074

Sun, Q. et al. PySCF: the Python-based simulations of chemistry framework. WIREs Comput. Mol. Sci. 8, e1340 (2018). https://doi.org/10.1002/wcms.1340

pytest

Methods we implemented from the literature

These are not libraries, they are the papers whose formulas are transcribed into vibe-qc’s source code. If you look at the gradient code and wonder where the energy-weighted density matrix came from, the answer is Pople, Krishnan, Schlegel and Binkley 1979. If you look at the Becke partitioning in the grid module, the answer is Becke 1988.

Hartree-Fock and Roothaan equations

Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory. Dover, 1996.

The Szabo-Ostlund textbook is what most of us learned this from. It is the book that makes the Roothaan equations, $\mathbf{F}\mathbf{C} = \mathbf{S}\mathbf{C}\boldsymbol{\varepsilon}$, actually click.

DIIS (Direct Inversion of the Iterative Subspace)

Pulay, P. Convergence acceleration of iterative sequences. The case of SCF iteration. Chem. Phys. Lett. 73, 393-398 (1980).

Pulay, P. Improved SCF convergence acceleration. J. Comput. Chem. 3, 556-560 (1982).

The commutator error measure $\mathbf{e} = \mathbf{F}\mathbf{D}\mathbf{S} – \mathbf{S}\mathbf{D}\mathbf{F}$ drives the extrapolation. The 5-to-6x speedup we measured on H₂O (52 iterations down to 9) is pure Pulay.

Analytic HF gradients (Pople-Binkley formulation)

Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Derivative studies in Hartree-Fock and Moller-Plesset theories. Int. J. Quantum Chem. Symp. 13, 225-241 (1979).

The gradient formula from this paper is

$$\frac{dE}{dR_A} = \frac{dE_{\mathrm{nuc}}}{dR_A} + \mathrm{tr}\!\left(\mathbf{D}\frac{d\mathbf{H}}{dR_A}\right) + \frac{1}{2}\mathrm{tr}\!\left(\mathbf{D}\frac{d\mathbf{G}(\mathbf{D})}{dR_A}\right) – \mathrm{tr}\!\left(\mathbf{W}\frac{d\mathbf{S}}{dR_A}\right)$$

where $\mathbf{W}$ is the energy-weighted density matrix. If you have ever wondered why the gradient code wants $\mathrm{tr}(\mathbf{W}\cdot d\mathbf{S})$ instead of something involving orbital derivatives directly, this is why.

DFT integration grid

Treutler, O.; Ahlrichs, R. Efficient molecular numerical integration schemes. J. Chem. Phys. 102, 346-354 (1995). https://doi.org/10.1063/1.469408

The M4 radial mapping with Chebyshev-2nd-kind nodes. Treutler-Ahlrichs is the standard choice for a reason.

Becke, A. D. A multicenter numerical integration scheme for polyatomic molecules. J. Chem. Phys. 88, 2547-2553 (1988). https://doi.org/10.1063/1.454033

Becke’s fuzzy-cell partitioning with the iterated (3/2)μ − (1/2)μ³ switch function is what makes the grid work for polyatomics.

Basis sets

Pritchard, B. P.; Altarawy, D.; Didier, B.; Gibson, T. D.; Windus, T. L. New Basis Set Exchange: An Open, Up-to-date Resource for the Molecular Sciences Community. J. Chem. Inf. Model. 59, 4814-4820 (2019). https://doi.org/10.1021/acs.jcim.9b00725

libint ships a snapshot of the Basis Set Exchange collection. 90 sets including the STO-3G, 3-21G, 6-31G, 6-311G, cc-pVXZ, and def2 families, plus the JKFIT, JFIT, and CABS auxiliary sets.

vibe-qc also has infrastructure for dropping in the pob-TZVP family, which was developed by our group for periodic solid-state calculations but applies equally well on the molecular side:

Peintinger, M. F.; Oliveira, D. V.; Bredow, T. Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations. J. Comput. Chem. 34, 451-459 (2013). https://doi.org/10.1002/jcc.23153

Once the raw .g94 files are dropped into basis_library/custom/, the setup script picks them up automatically.

Exchange-correlation functionals (exposed as named aliases)

  • Slater exchange: Slater, J. C. A Simplification of the Hartree-Fock Method. Phys. Rev. 81, 385 (1951).
  • VWN5 correlation: Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations. Can. J. Phys. 58, 1200 (1980).
  • PBE: Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865 (1996).
  • B88 exchange: Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 38, 3098 (1988).
  • LYP correlation: Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 37, 785 (1988).
  • B3LYP: Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 98, 11623 (1994).

Development tooling

  • Claude Code (Anthropic). The coding agent that paired with me on this. https://www.anthropic.com/claude-code
  • GitLab (self-hosted at gitlab.peintinger.com). Version control host.
  • Homebrew. macOS package manager used for libxc, Eigen, Boost, GMP, CMake, and Ninja.

Closing thought

The list above is long, and that is the point. “Vibe-coding a quantum chemistry program in one day” is a true description of what happened, but it is a misleading one if you read it as “a quantum chemistry program was created from nothing in one day.” What actually happened is that a coding agent (pairing with someone who knows the field) threaded together the work of hundreds of authors spanning seventy-five years of physics, chemistry, numerical analysis, and software engineering.

The compressible part was the glue. The giants were already there, standing still, waiting to be stood on.

I vibe-coded a quantum-chemical program in one day

I vibe-coded a quantum-chemical program in one day

During my PhD at the Mulliken Center for Theoretical Chemistry in Bonn, I wrote a quantum-chemical code called AICCM, an object-oriented educational implementation of the Cyclic Cluster Model (CCM) at the Hartree-Fock level. If you want the formal reference, it is described in:

M. F. Peintinger, T. Bredow, The cyclic cluster model at Hartree-Fock level, Journal of Computational Chemistry 35 (11), 839-846 (2014). DOI: 10.1002/jcc.23550

That work took a good chunk of my PhD. Years of implementation, debugging, validating against CRYSTAL and other periodic codes, reading obscure papers on Wigner-Seitz integration and boundary handling of three- and four-center integrals. Real work. The kind that leaves scars.

Alongside AICCM I also developed the pob-TZVP basis sets, a consistent triple-zeta valence family tuned specifically for periodic solid-state calculations. They were built for CRYSTAL, but from the start I had AICCM in mind as a second home for them:

M. F. Peintinger, D. V. Oliveira, T. Bredow, Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations, Journal of Computational Chemistry 34 (6), 451-459 (2013). DOI: 10.1002/jcc.23153

These basis sets matter for the story later in this post, because vibe-qc already has a slot in its basis library waiting for them.

This week I listened to a podcast on vibe-coding, and the thought stuck in my head: what happens if I hand a modern coding agent a task I genuinely understand, something deeply non-trivial, and then resist the urge to help code? Not a toy CRUD app. Not a weekend script. A quantum chemistry program. The kind of thing where a single wrong sign in a density matrix silently gives you garbage for three hours of debugging before you notice.

So I created a blank git repo on my GitLab server, called it vibe-qc, and gave Claude the task. Just the task. No starter code, no reference implementation to crib from, no hand-holding on the math. I steered at the architecture level and made scoping decisions, but I did not write code. Not one line.

What I expected vs. what happened

I honestly expected it to crash out somewhere around the two-electron integrals, or to produce a plausible-looking SCF loop that silently disagreed with the reference by a few millihartree. That is the usual failure mode when someone writes this code for the first time: it runs, it looks right, and the total energy is wrong in the fourth digit.

Instead, by the end of the day we had a working program that agreed with PySCF to machine precision. It did fall into one trap I found amusing: at one point it tried to run restricted Hartree-Fock on lithium. Lithium has three electrons. You cannot do closed-shell RHF on an odd number of electrons. You need an open-shell method. We added UHF a few commits later and the problem went away, but it is telling that an AI coding agent can reproduce the exact kind of mistake a first-year grad student makes.

What got built in one day

16 commits. Here is the inventory.

Architecture

Python frontend with a C++17 numerical core, bridged by pybind11. CMake plus scikit-build-core for the build system, so pip install -e . just works. About 3500 lines of C++ for integrals, SCF drivers, gradients, and DFT machinery. libint2 v2.13.1 built from source under third_party/libint/ with max_am=5 and first-order derivatives enabled. Homebrew libxc 7.0 for the 500+ exchange-correlation functionals.

The GIL is released during native work, so the Python frontend is not a bottleneck. I had asked about this up front and was reassured that Python would not become a restriction for multi-core scaling later.

Hartree-Fock, closed and open shell

  • RHF with Hcore or SAD initial guess, symmetric orthogonalization via $\mathbf{S}^{-1/2}$, optional density damping.
  • UHF with separate alpha and beta densities, reporting <S²> as a spin-contamination diagnostic. Closed-shell UHF collapses to RHF at 1e-14.
  • Pulay DIIS convergence accelerator. On H₂O we measured 52 iterations (with damping 0.5) collapsing to 9 iterations (with DIIS). Per-spin DIIS for UHF.
  • SAD initial guess (Superposition of Atomic Densities) that runs a fractional-occupation atomic SCF per unique element. This was added specifically after UHF hit a wrong local minimum on OH in 6-31G*. The fix was immediate.
  • Analytic nuclear gradients for both RHF and UHF via the Pople-Binkley formula with an energy-weighted density matrix. Verified against PySCF at 1e-8 Ha/bohr.

Density Functional Theory

  • Numerical integration grid: Treutler-Ahlrichs M4 radial (Chebyshev-2nd-kind nodes) combined with a Gauss-Legendre-in-cos(θ) and uniform-φ angular product, partitioned across atoms by Becke fuzzy cells with the iterated switch function. The default medium grid integrates the H-1s density to 1e-10.
  • AO evaluation on the grid: both $\chi_\mu(\mathbf{r})$ and $\nabla\chi_\mu(\mathbf{r})$ at every grid point, supporting Cartesian or pure spherical shells. Analytic AO gradients verified against finite differences at 1e-10.
  • libxc wrapper accepting alias names (“LDA”, “PBE”, “B3LYP”) or explicit comma-separated XC_… IDs, with hybrid fractions detected automatically.
  • RKS SCF driver with full energy decomposition $(E_{\mathrm{core}},\, E_J,\, \alpha E_K,\, E_{\mathrm{xc}},\, E_{\mathrm{nuc}})$. Matches PySCF to 5e-11 Ha on LDA and to grid accuracy (~1e-6 Ha) on B3LYP and PBE.

ASE integration

This was the decision I am most happy about in retrospect. Rather than write a custom input parser, we wired the code into ASE (the Atomic Simulation Environment) as a proper Calculator subclass. That gave us, for free:

  • Input readers for XYZ, CIF, PDB, POSCAR, Gaussian inputs, and basically every format chemists actually use.
  • Geometry optimization via ase.optimize.BFGS. H₂O in HF/STO-3G with a distorted starting geometry converges in 7 BFGS steps to r(OH) = 0.989 Å and angle HOH = 100.0°, exactly the known literature value.
  • Implicit method routing: HF with multiplicity 1 goes to RHF with forces, HF with multiplicity > 1 goes to UHF with forces, DFT with multiplicity 1 goes to RKS (energy-only at day’s end).
  • Clean unit handling via ase.units.Bohr and ase.units.Hartree.

Basis set library

90 standard Gaussian .g94 basis sets shipped by libint are exposed automatically (STO-3G through aug-cc-pVQZ, the def2 family, plus the JKFIT, JFIT, and CABS auxiliary sets). On top of that there is a basis_library/custom/ directory where a user can drop their own .g94 files. A setup script assembles the union, with custom files overriding standard ones by name. This path is ready for dropping in my own pob-TZVP and pob-TZVP-D basis sets (the ones I developed for solid-state calculations) once I dig up the files.

Tests

131 tests across 13 files, full suite runs in about 10 seconds. Every numerical feature is cross-checked against PySCF as the reference. A representative slice of what actually passes:

  • H₂ / STO-3G at R = 1.4 bohr: E_HF = −1.116714325063 Ha, difference from PySCF = −1.87e-14.
  • H₂O / STO-3G at experimental geometry: E_HF = −74.964453863067 Ha, difference from PySCF = −2.84e-14.
  • OH doublet / STO-3G UHF: <S²> = 0.7533 (ideal 0.75, about 0.3% spin contamination).
  • H₂O / STO-3G / LDA: 5e-11 Ha vs PySCF.
  • H₂O / STO-3G / B3LYP: 1.3e-7 Ha (grid-accuracy limited).
  • CH₄ / 6-31G* / PBE: under 3e-6 Ha.

Documentation

A README with quick-start examples for both macOS and Linux, a detailed installation doc covering Homebrew, Debian/Ubuntu, Fedora/RHEL, and Conda paths, a quickstart tour, and a basis-library README. All written by the same agent, in one day, alongside the code.

What I deliberately did not add yet

UKS (unrestricted Kohn-Sham), analytic DFT gradients, ROHF, OpenMP inside the Fock and ERI-gradient loops, MPI for multi-node, post-HF correlation methods like MP2 and CCSD, and of course the actual cyclic cluster model. Every foundation for CCM is now in place: libint with periodicity hooks available, ASE’s CIF reader already wired in, and a validated HF/DFT SCF stack on finite molecules that I trust.

How I steered it

If you are thinking about trying this yourself, the interesting part was not the code, it was the scoping. A few patterns I noticed in my own behavior over the day:

Anticipatory infrastructure. Early in the morning, when the code was still just producing one-electron integrals, I made the decision to build libint from source with first-order derivatives enabled. Homebrew’s libint would have been faster to set up, but it did not include the derivative integrals. I knew we would need gradients for geometry optimization later, and I did not want to tear down and rebuild the install halfway through the day. That single 10-minute decision paid off twice: once for HF gradients, once for UHF gradients.

Validation before features. After RHF produced a plausible number, I pushed for formal pytest coverage against PySCF before we added DIIS. That caught two real bugs immediately: a Cartesian-vs-spherical d-orbital basis-purity issue, and a last-iteration MO self-consistency bug. Both were silent errors. Neither would have shown up if I had just eyeballed the total energy.

End-user surface before internals. We wired in the ASE Calculator and the logging integration before adding forces, before DFT, before SAD. The moment I caught the calculator running silently without emitting any trace, I paused the next phase and inserted a logging phase. It is much easier to debug a 50-iteration SCF when you can see the DIIS error norm falling.

Realistic descoping. My opening plan was “UHF and ROHF, then DFT.” By mid-afternoon it was “UHF only, ROHF can wait.” On the DFT grid, I accepted a Gauss-Legendre product angular grid instead of a proper Lebedev grid for the MVP. These were not compromises on correctness, they were compromises on sophistication, and they let us ship a working DFT implementation the same day.

Commit cadence. Every phase got its own focused commit, with a descriptive message, passing the full test suite. No force-pushes, no amends. 16 clean revertable checkpoints.

Standing on shoulders

None of this would have been possible in a day if the hard parts had not already been solved by other people. The integral engine is libint (Ed Valeev), the 500+ DFT functionals come from libxc (Susi Lehtola and co.), linear algebra is Eigen, the Python bindings are pybind11, the user-facing surface is ASE (Larsen et al.), and every single number in the test suite is cross-checked against PySCF (Qiming Sun et al.) as the reference oracle. On the methods side we are implementing formulas from Pulay (DIIS), Pople, Krishnan, Schlegel and Binkley (analytic HF gradients), Treutler and Ahlrichs (radial grid), and Becke (atomic partitioning). Full citations and a guided tour of the stack are in the companion post.

Watching it work was the best part

I need to say this plainly: I had a blast. I spent the day with Anthropic’s Claude Opus 4.7 as a pair-programmer, and watching the model work through this problem was genuinely thrilling. Every time I came back from a coffee or a meeting, there were new commits, clean ones, with descriptive messages, and a test suite that had grown and still passed. The SCF converged on H₂O. Then DIIS landed and the iteration count dropped from 52 to 9. Then gradients showed up and BFGS walked the water molecule to its equilibrium geometry. Then DFT arrived and B3LYP matched PySCF to grid accuracy. Each of those moments felt like a little gift.

What took me multiple years during my PhD, Opus 4.7 produced in a day. Not because vibe-qc is AICCM, it is not. AICCM implements real periodic CCM with four-center integrals at Wigner-Seitz boundaries, which is a genuinely hard problem and is not in vibe-qc yet. But the molecular HF + DFT + gradients + ASE integration layer, which in a traditional PhD timeline is itself a one- or two-year project, collapsed into a single afternoon. That is a kind of leverage I have never experienced before, and I am grinning about it.

The PhD was not wasted, quite the opposite. The reason I could steer this build effectively, and the reason I knew to demand a SAD guess the moment UHF stalled, and the reason I caught the lithium-on-RHF mistake instantly, is because I spent those years learning quantum chemistry from the inside. The expertise did not become obsolete. The expertise became the steering wheel, and it turns out steering is a deeply satisfying way to use what you know.

The really exciting part is what this unlocks. Ideas that used to be “interesting but I cannot justify six months to try it” suddenly become “let us see what happens this weekend.” The research questions I filed away after my PhD, the ones I always meant to come back to, are now within reach in a way they were not last year. That is a good day.

The next step, when I find another day, is the cyclic cluster model. That is the part where I find out whether the model can actually extend beyond well-trodden textbook ground into the corner of the literature that I personally carved out a decade ago. I am genuinely curious. I will report back.

The code lives at vibe-qc.com. If you want to reproduce the experiment yourself, the recipe is: a blank repo, a clear scoping brief, and the discipline to let the model cook.