TREXIO for CI wavefunctions and periodic restarts

After the molecular exchange tutorial, two further checks matter: keeping CI coefficients attached to their actual orbital basis, and keeping a periodic density attached to its k-point sampling. Both examples below retain the default QVF output and request TREXIO explicitly. Run them from a source checkout using an environment with vibe-qc[trexio] installed.

LiH with a frozen core and an active space

This example uses LiH at a separation of 3 bohr, STO-3G, and CASSCF with two active electrons in two active orbitals. One spatial core orbital remains doubly occupied. This is a small interchange model, not a basis-converged LiH benchmark.

export OMP_NUM_THREADS=1
python examples/trexio/correlated_wavefunction.py --method casscf \
    --output-dir ~/vibeqc-runs/trexio/lih-casscf
python examples/trexio/correlated_wavefunction.py --method casscf --backend text \
    --output-dir ~/vibeqc-runs/trexio/lih-casscf-text

Download the complete correlated_wavefunction.py. Its calculation is:

molecule = Molecule([Atom(3, [0., 0., 0.]), Atom(1, [0., 0., 3.])])
options = {} if args.method == "fci" else {"active_space": (2, 2)}
result = run_job(molecule, basis="sto-3g", method=args.method,
                 output=stem, trexio=True, trexio_backend=args.backend,
                 progress=False, verbose=False, **options)
suffix = ".trexio.h5" if args.backend == "hdf5" else ".trexio"
path = Path(str(stem) + suffix)
data = read_trexio(path)

The automatic export contains the determinant expansion and optimized CASSCF orbitals together. Replacing those orbitals with the original RHF coefficients changes the represented wavefunction even if the determinant coefficients are unchanged. Frozen core occupation must also be present in every determinant.

The script checks normalization of the CI vector, spin-resolved RDM traces, the four-electron AO density, and the frozen-core bits:

coefficients = np.asarray(data.fields["determinant_coefficient"])
norm = float(np.vdot(coefficients, coefficients).real)
np.testing.assert_allclose(norm, 1., rtol=0, atol=1e-10)
assert data.mo_type == "CI"
assert "mo_energy" not in data.fields  # no invented canonical SCF energies
np.testing.assert_allclose(data.energy, result.energy, rtol=0, atol=1e-10)
for spin, count in (("up", data.n_up), ("dn", data.n_dn)):
    np.testing.assert_allclose(np.trace(data.fields["rdm_1e_" + spin]), count,
                               rtol=0, atol=1e-9)
density = data.density_matrices()[0]  # uses the full correlated 1-RDM
overlap = compute_overlap(data.basis_set())
np.testing.assert_allclose(np.trace(density @ overlap), 4., rtol=0, atol=1e-9)
if args.method != "fci":
    assert data.mo_class.count("Core") == 1
    assert data.mo_class.count("Active") == 2
    # Six spatial orbitals fit in one 64-bit word per spin. Core MO 0
    # must be occupied in BOTH words of every determinant.
    assert all(int(word) & 1 for det in data.fields["determinant_list"] for word in det)

For this small example there are six spatial MOs, so one 64-bit word per spin suffices. Do not copy its one-word bit check into a general determinant reader without handling larger orbital spaces.

The total energy should be about -7.881214343110 Ha. The output ends in PASS casscf, with four determinants and CI norm close to one. Orbital energies are absent because canonical SCF eigenvalues would not describe this optimized correlated orbital set. The density comes from the full stored 1-RDM, including its off-diagonal elements.

Change --method to casci for a fixed-orbital active-space calculation or fci for all four electrons in all six orbitals. FCI does not use the active_space argument and is only inexpensive here because the model is tiny. The script accepts both backends for every method.

Independent CI energy reconstruction

With PySCF installed in the reference interpreter:

python examples/regression/runner_trexio_pyscf.py \
    ~/vibeqc-runs/trexio/lih-casscf/lih-casscf.trexio.h5

The checker forms a Hamiltonian from its own integrals and contracts it with the stored CI vector. Its result must report "wavefunction": "CI", "verdict": "pass", and an energy difference below 1e-8 Ha. It does not rerun a separate CASSCF optimization. That distinction lets it test whether the file describes the original state.

The job exporter selects root zero. Use the low-level root= option for a different computed state; see CI and CASSCF wavefunctions. A CI density can seed an SCF calculation, but READ is not a continuation of the determinant solver or CASSCF optimizer. Coupled-cluster amplitudes and higher RDMs require explicitly supplied fields rather than automatic export.

A periodic two-k-point wavefunction

The periodic example uses one He atom in a cubic cell of side 7 bohr, RHF/STO-3G, and a (2, 1, 1) k-point mesh. It deliberately uses a small GDF cutoff and disables automatic cutoff convergence to keep the IO exercise small. SCF convergence is still required. Its energy is not a converged prediction for solid helium.

python examples/trexio/periodic_restart.py \
    --output-dir ~/vibeqc-runs/trexio/helium
python examples/trexio/periodic_restart.py --backend text \
    --output-dir ~/vibeqc-runs/trexio/helium-text

Download periodic_restart.py. The same settings are used for the original calculation and the restart:

system = PeriodicSystem(3, np.eye(3) * 7., [Atom(2, [0., 0., 0.])])
basis = BasisSet(system.unit_cell_molecule(), "sto-3g")
settings = dict(method="RHF", jk_method="gdf", kpoints=(2, 1, 1),
                aux_basis="def2-svp-jk", gdf_method="rsgdf", rsgdf_ke_cutoff=12.,
                convergence="off", max_iter=60, conv_tol_energy=1e-8,
                progress=False, verbose=0,
                write_molden_file=False, write_population_file=False)
result = run_periodic_job(system, basis, output=args.output_dir / "helium",
                          trexio=True, trexio_backend=args.backend, **settings)
suffix = ".trexio.h5" if args.backend == "hdf5" else ".trexio"
path = args.output_dir / ("helium" + suffix)
data = read_trexio(path)

Keep the sampling with the orbitals

assert result.converged
assert data.periodic
blocks = data.mo_blocks()
assert [(block.k_point, block.spin) for block in blocks] == [(0, 0), (1, 0)]
np.testing.assert_allclose(data.lattice, system.lattice, rtol=0, atol=1e-12)
np.testing.assert_allclose(sum(data.kpoint_weights), 1., rtol=0, atol=1e-12)
electrons = sum(data.kpoint_weights[block.k_point] * sum(block.occupations)
                for block in blocks)
np.testing.assert_allclose(electrons, 2., rtol=0, atol=1e-10)
for block, original in zip(blocks, result.mo_coeffs):
    np.testing.assert_allclose(block.coefficients, original, rtol=0, atol=1e-12)

There are two restricted orbital blocks, identified by k_point 0 and 1. The per-cell electron count is \(N=\sum_k w_k\sum_i f_{ik}=2\), not the unweighted sum over both blocks. For an unrestricted calculation, include both spin blocks at each k point. The reader’s densities are also unweighted; apply the weights when integrating a cell quantity.

data.lattice uses column vectors in bohr. data.kpoints gives Cartesian inverse-bohr vectors, while raw pbc_k_point holds reduced coordinates. Complex coefficients remain complex even when this particularly small case happens to have real-valued orbitals. A molecular overlap matrix is not the Bloch overlap at a nonzero k point; the molecular overlap check from the first tutorial must not be reused here.

READ with the same target mesh

restarted = run_periodic_job(data.periodic_system(), data.basis_set(),
                             output=args.output_dir / "helium-read",
                             initial_guess="read", read_from=path, **settings)
assert restarted.converged
error = abs(restarted.energy - result.energy)
assert error < 1e-8, error

The output should end in PASS periodic, a weighted electron count of 2, and a READ energy difference below 1e-8 Ha. The second run creates helium-read.qvf and the ordinary job files. The first run’s TREXIO file remains the explicit restart source.

The target lattice, basis, k points and weights must be compatible with the source; READ can map a different ordering of the same mesh, but it does not interpolate a density onto a new mesh. Multi-k AO integral arrays are omitted because the corresponding TREXIO datasets have no k-index dimension. Foreign one- or two-dimensional files need an explicit dim= when calling periodic_system() unless they carry vibe-qc’s dimension metadata.

The molecular PySCF checker rejects periodic files. This example verifies preservation and restart consistency; it does not establish independent periodic integral agreement. For physical convergence, continue with periodic SCF convergence and Bloch orbitals and k points.