Cohesive energies and the energy-engine abstraction

The cohesive (atomization) energy of a solid is the energy needed to pull it apart into free atoms. It is the quantity the solid-state basis-optimisation work optimises against, and computing one is not a single calculation but a sequence: relax the cell, take a single point, correct for zero-point motion, compute the free atoms, and assemble.

vibe_basis.pipeline runs that sequence. vibe_basis.engine is the seam that decides which program runs each step.

Note

These live in the co-located vibe-basis driver, so they need vibe-qc’s basisopt extra. See The external-program path. The plan they implement is vibe-basis/ROADMAP.md.

The engine abstraction

Every engine implements one protocol, EnergyEngine, so which program computes a number is configuration rather than structure:

Role

Engine

Where it lives

Primary, Gaussian

VibeQcEngine

vibeqc.basis_optimization.calculators

Primary, plane wave

vibe-qc GPW / GAPW

not wired yet

Fallback, Gaussian

Crystal23Engine

vibe_basis.engines.crystal23

Fallback, plane wave

GpawEngine

vibe_basis.engines.gpaw (declared, not built)

vibe-qc is the engine. The external ones are fallbacks and cross-checks. They exist to unblock work whose vibe-qc capability has not landed, and to give an independent second opinion: agreement between two implementations is evidence in a way that one implementation agreeing with itself is not.

The primary engine lives outside the driver package while the fallbacks live inside it. That inversion is deliberate. VibeQcEngine needs vibe-qc’s compiled core; Crystal23Engine needs only subprocess and file I/O. Keeping the external engines dependency-free is what makes pip install vibe-basis a ten-second install for a collaborator who has their own CRYSTAL and no interest in building vibe-qc.

Every number carries its provenance

An engine returns an EngineEnergy, not a bare float:

from vibe_basis.engines import Crystal23Engine
from vibe_basis.transports.local import LocalTransport
from vibe_basis.io.structures import STRUCTURES

engine = Crystal23Engine(transport=LocalTransport())
result = engine.crystal_energy("<inline basis text>", STRUCTURES["MgO"], method="pbe")

print(result.ok)              # did it converge and complete?
print(result.energy)          # Hartree per unit cell
print(result.engine)          # "crystal23"
print(result.engine_version)  # what the output banner actually said
print(result.failure_mode)    # None on success

“MgO is -275.4776 Ha” is not a result. “CRYSTAL23 computed -275.4776 Ha for MgO at RHF/pob-TZVP” is. A campaign that swaps engines as easily as this one does has to record which one produced each point, or it cannot be audited afterwards.

If you only want a number, evaluate_crystal and evaluate_atom return float | None. They are derived from the rich methods, not a parallel implementation, so the two views cannot disagree.

VibeQcEngine defaults are wiring defaults, not converged settings

Two of them will bite a first campaign if taken at face value:

  • Gamma point only. The periodic single point runs at a single k-point unless you say otherwise. That is enough to check the wiring and far too coarse for a production number. Pass a converged mesh through periodic_kwargs.

  • BIPOLE, not GDF. GDF needs an auxiliary fitting basis, and a freshly emitted candidate basis has none registered. Borrowing another basis’s JKFIT set would lay a fitting error on top of the basis-incompleteness signal the objective exists to measure, so the route that needs no auxiliary basis is the default. It is also the route the relaxation drives. Override through periodic_kwargs if you have a fitting set you trust for the basis in hand.

The free-atom reference: aspherical and spin-polarised

Half of a cohesive energy is the free atoms, and an atom reference is only comparable to an external one computed under the same convention. Crystal23Engine emits SYMMREMO, UHF (or SPIN in the DFT block) and SPINLOCK for an open-shell atom, so the solution is aspherical and spin-polarised; VibeQcEngine runs the spin-polarised ground-state multiplicity. Both record which convention they used in EngineEnergy.detail["atom_reference"].

This is worth checking rather than assuming, because the failure is silent: a spin-restricted atom SCF converges perfectly well and simply returns the wrong energy. In the plane-wave reference work, spin-restricting the bromine reference alone moved KBr’s atomization energy by 37 kJ/mol, which is far larger than the basis-set effects a campaign is trying to resolve.

Counterpoise: matching the published pob protocol

Every free-atom deck in the pob cohesive-energy reference set uses CRYSTAL’s ATOMBSSE, meaning the atom is computed in the ghost basis of its own crystal rather than bare. The published numbers are therefore counterpoise corrected, and a cohesive energy assembled from bare atoms is not the same quantity.

# doc-audit: skip
pipeline = CohesivePipeline(engine, method="r2scan", counterpoise=True)
result = pipeline.run(basis_text, structure)
assert result.counterpoise_applied

The difference between the two is the basis-set superposition error, which for an incomplete basis is exactly what a basis-optimisation campaign is trying to reduce. Leaving it out does not cancel when you compare against a corrected reference.

Two consequences:

  • It is off by default, and CohesiveResult.counterpoise_applied records which you got. An engine that cannot do it (currently VibeQcEngine) produces a recorded skipped stage rather than quietly falling back to bare atoms.

  • Counterpoise atoms are not shared across systems. A bare oxygen is the same calculation for MgO and CaO, so the cache computes it once. A counterpoise oxygen is computed in its own host’s ghost basis, so the two are genuinely different runs and the cache keys them by host. That is a real cost: it removes the saving that makes a campaign affordable, so use it for reference comparison rather than for every point of an optimisation.

nstar and rmax control how far the ghost shell extends. The reference set uses per-system converged values (LiF 30 / 10.0, MgO 10 / 5.0), so treat Crystal23Engine’s defaults as a starting point to converge rather than a setting to trust.

Failures are returned, not raised

A failed evaluation comes back with ok=False and a failure_mode. It does not raise. Basis optimisation walks into infeasible regions routinely: an exponent collapses, the overlap goes singular, the queue drops a job. A derivative-free optimiser handles that by scoring the point inf and moving on, so an exception would abort a campaign that is minutes per evaluation and days long.

One consequence worth internalising: energy_if_ok() returns None for a failed evaluation even when .energy holds a number. A non-converged SCF’s last-cycle energy is kept for diagnosis, and silently consuming it as if it were converged is exactly the failure that method exists to prevent.

The pipeline

from vibe_basis.cache import EnergyCache
from vibe_basis.pipeline import CohesivePipeline
from vibe_basis.engines import Crystal23Engine
from vibe_basis.transports.local import LocalTransport
from vibe_basis.io.structures import STRUCTURES

pipeline = CohesivePipeline(
    Crystal23Engine(transport=LocalTransport()),
    method="pbe",            # bulk single point and the free atoms
    relax_method="pbe",      # geometry relaxation, if the engine can
    relax=True,
    zero_point=True,
    cache=EnergyCache("runs/campaign.jsonl"),
)

result = pipeline.run("<inline basis text>", STRUCTURES["MgO"])
print(result.summary())

Stages, in order:

  1. Relax the cell and internal coordinates at relax_method, so the basis is judged at its own equilibrium geometry rather than at one some other basis produced.

  2. Single point at method, on the relaxed geometry. The usual pattern is a cheap GGA relaxation followed by a hybrid single point.

  3. Zero-point energy of the solid.

  4. Isolated atoms, one per distinct element.

  5. Assemble.

The arithmetic, per formula unit:

E_coh  =  SUM_atoms E_atom  -  E_bulk_per_formula_unit  -  ZPE_per_fu

positive for a bound solid, matching the sign convention of vibeqc.atomization.AtomizationResult. Zero-point motion reduces the cohesive energy: the solid’s true ground state sits ZPE above its static-lattice minimum, while free atoms have no vibrations to correct.

Formula units are counted from the cell’s stoichiometry, as the greatest common divisor of the per-element counts. Rocksalt MgO with 4 Mg and 4 O is 4 formula units; fluorite CaF2 with 4 Ca and 8 F is also 4. Pass n_formula_units= explicitly for a cell whose formula is not the reduced one, such as a defect supercell.

Read zero_point_applied before comparing to anything

# doc-audit: skip
if not result.zero_point_applied:
    ...  # this is a STATIC-LATTICE cohesive energy

A static-lattice cohesive energy and a ZPE-corrected one are different quantities, and for these solids they differ by tens of kJ/mol, which is the same size as the basis-set effects being measured. Comparing one to a reference computed as the other produces a discrepancy that looks like physics and is not.

So CohesiveResult states which it is, and the pipeline never quietly substitutes a default. If an engine cannot relax, or cannot do phonons, that stage is recorded as skipped, with the engine named in the reason:

# doc-audit: skip
for stage in result.stages:
    print(stage.name, stage.status, stage.reason)
# relax        skipped  engine 'crystal23' has no relaxer
# single_point ok       None
# zero_point   skipped  engine 'crystal23' has no phonons
# atom_Z8      ok       None
# atom_Z12     cached   None

Returning the input geometry for a missing relaxation would look like a converged relaxation; returning zero for a missing ZPE would look like a solid with no vibrations. Both would produce a number that looks finished and is not.

A failed zero-point stage is survivable, and downgrades the result to static-lattice with zero_point_applied=False. A failed relaxation is not: an unconverged geometry is not a geometry, so the run stops before spending an SCF on it.

Caching and resume

EnergyCache is content-addressed over an append-only JSONL journal. Pass one cache across every system in a campaign:

# doc-audit: skip
cache = EnergyCache("runs/campaign.jsonl")
for name in ("MgO", "CaO", "BaO"):
    CohesivePipeline(engine, method="pbe", cache=cache).run(basis_text, STRUCTURES[name])

Free-atom entries are keyed without a system identifier, so the oxygen computed for MgO is reused for CaO and BaO. That is the win that matters: atoms are the expensive part. In the plane-wave reference work a single free fluorine atom cost 100 to 240 seconds while the solid was quick, so recomputing oxygen once per compound would dominate a campaign.

Because the journal is append-only and is the cache, a campaign killed mid-run resumes rather than re-spending the queue time. A journal truncated by kill -9 costs the one torn entry, not the file.

Keys derive from declared configuration only, never from observed output. In particular the engine’s reported version is deliberately not part of the key: Crystal23Engine cannot report one until its first evaluation, so keying on it would change the key halfway through a run and silently split the cache. The observed version is recorded in the value instead, for audit.

Which CRYSTAL produced this?

Crystal23Engine targets CRYSTAL23, and the parsed result records the major version read off the output banner:

from vibe_basis.backends.crystal import parse_output_file

parsed = parse_output_file("mgo.out")
print(parsed.crystal_version)        # 23, or None if no banner in the text
print(parsed.version_matches(23))    # True / False / None

version_matches is three-valued on purpose. None means the banner was not in the text, which is normal for a truncated output or a tail slice, and is not a failure: an output with no visible banner but a converged total energy is still a usable energy. False means the banner said something else, and the energy is still parsed, but attributing it to CRYSTAL23 would be false.

Enforcement is a policy decision, so it is opt-in. Set Crystal23Engine(..., strict_version=True) for a production campaign, where a stray binary on one host quietly contributing a few points is the real hazard. Leave it off for a parity study that legitimately wants CRYSTAL14 output.

See also

  • Molecular basis-set optimisation, the in-process analytic-gradient path.

  • Basis sets, the bundled library.

  • vibe-basis/ROADMAP.md, the milestones this implements.

  • vibe-basis/TUTORIAL.md, a worked walkthrough.

  • examples/basisset_dev/cohesive_pipeline_demo.py, runnable with no external program and no native build.