COOP/COHP: bonding analysis for periodic systems¶

COOP (Crystal Orbital Overlap Population) and COHP (Crystal Orbital Hamilton Population) are energy-resolved projections of the overlap or Hamiltonian matrix onto atom-pair subspaces. They are the standard tools for bonding/antibonding analysis of periodic systems – the solid-state analogues of Mayer bond orders, resolved by energy rather than summed into a single scalar.

Feature status

vibeqc.coop_cohp.compute_coop_cohp() ships COOP and COHP for restricted and spin-polarized (UHF/UKS) periodic calculations. Both CPU (C++ kernels) and Python fallback paths are available. Output is a COOPCOHPResult with energy-resolved curves and integrated ICOOP / -ICOHP values. The Lobster convention (-COHP) is used so bonding states appear positive in both COOP and COHP.

What COOP and COHP are¶

At each k-point, for band n and atom pair (A, B):

COOP_AB(E) = Σ_k w_k Σ_n [ Σ_{μ in A, ν in B} S_{μν}(k) * C*_{μn}(k) * C_{νn}(k) ] * δ(E - ε_n(k))
COHP_AB(E) = Σ_k w_k Σ_n [ Σ_{μ in A, ν in B} H_{μν}(k) * C*_{μn}(k) * C_{νn}(k) ] * δ(E - ε_n(k))

where the delta function is broadened by a Gaussian of width sigma. The energy integral up to E_F gives ICOOP / -ICOHP – the bond strength.

Negative COHP = bonding
Positive COHP = antibonding
Following the Lobster convention, -COHP is returned so bonding states appear positive (same sign as COOP)

Quick start¶

from vibeqc.coop_cohp import compute_coop_cohp
from vibeqc import run_rhf_periodic_multi_k_ewald3d, KPoints

# Run a periodic SCF
system, basis, kmesh = ...  # your PeriodicSystem, BasisSet, KPoints
result = run_rhf_periodic_multi_k_ewald3d(system, basis, kmesh)

# Compute COOP/COHP from the SCF result
cohp_result = compute_coop_cohp(
    F_terms=result.fock_terms,
    S_real=result.overlap_real,
    system=system,
    basis=basis,
    kmesh=result.kmesh,
    H_terms=result.hcore_terms,  # optional, omit for COOP-only
    n_electrons_per_cell=system.n_electrons_per_cell(),
    sigma=0.01,
)

# cohp_result.energies  -- energy grid (Ha)
# cohp_result.coop      -- COOP(E), shape (n_pairs, n_e)
# cohp_result.cohp      -- -COHP(E), same shape (or None if no H_terms)
# cohp_result.integrated_coop  -- ICOOP, (n_pairs,)
# cohp_result.integrated_cohp  -- -ICOHP, (n_pairs,)
# cohp_result.pairs     -- list of {i, j, symbol_i, symbol_j, distance_ang}

Spin-polarized COOP/COHP¶

Pass F_terms_beta (the beta-spin Fock terms) to compute separate alpha and beta projections:

result = compute_coop_cohp(
    F_terms=result_uhf.fock_terms,
    F_terms_beta=result_uhf.fock_terms_beta,
    S_real=result_uhf.overlap_real,
    system=system,
    basis=basis,
    kmesh=kmesh,
    H_terms=result_uhf.hcore_terms,
    n_electrons_per_cell=system.n_electrons_per_cell(),
)
# result.coop.shape  -- (2, n_pairs, n_e)  [alpha, beta]
# result.cohp.shape  -- (2, n_pairs, n_e)

AO-pair grouping¶

The ao_pairs_per_atom_pair() helper partitions AO indices into atom-pair subspaces. Pairs farther apart than pair_distance_cutoff (default 8.0 bohr) are excluded:

from vibeqc.coop_cohp import ao_pairs_per_atom_pair

pairs = ao_pairs_per_atom_pair(system, basis, pair_distance_cutoff=8.0)
# {(atom_a, atom_b): (ao_indices_a, ao_indices_b), ...}

Where to look in the code¶

Layer	File
Python API + compute	`python/vibeqc/coop_cohp.py`
C++ weight kernels	`_coop_weights_k`, `_cohp_weights_k` in `cpp/src/bindings.cpp`
Gaussian broadening	`_gaussian_broaden_projected` in `cpp/src/bindings.cpp`
Citation route	`python/vibeqc/output/citations/database.toml` (routes: `coop`, `cohp`)
Tests	`tests/test_coop_cohp.py`

References¶

The computation follows the standard formulation described by:

COOP: Hughbanks & Hoffmann, J. Am. Chem. Soc. 105, 3528 (1983).
COHP: Dronskowski & Blochl, J. Phys. Chem. 97, 8617 (1993).
Modern plane-wave/LCAO bridge: Deringer, Tchougreeff & Dronskowski, J. Phys. Chem. A 115, 5461 (2011).
LOBSTER 3.0 (canonical implementation reference): Maintz, Deringer, Tchougreeff & Dronskowski, J. Comput. Chem. 37, 1030 (2016).

Citations are automatically included in the output manifest when COOP/COHP is computed; see docs/user_guide/citations.md.

Runner integration¶

The easiest way to get COOP/COHP is through the periodic runner:

from vibeqc import run_periodic_job

result = run_periodic_job(
    system, basis, method="RKS", functional="pbe",
    kpoints=[2, 2, 2], output_qvf=True, coop_cohp=True,
    dos_kmesh=[8, 8, 8],  # optional, controls DOS/COOP k-mesh
)

This computes COOP/COHP on the same DOS k-mesh and embeds dos.coop and dos.cohp sections in the .qvf archive. For UHF/UKS, the beta Fock channel is auto-detected and n_spin=2 sections are written.

Plotting¶

from vibeqc.plot import coop_figure, cohp_figure, bands_cohp_figure
from vibeqc import band_structure

bs = band_structure(system, basis, kpath)

# Standalone COOP/COHP plots
coop_figure(cohp_result)
cohp_figure(cohp_result)

# Combo: band structure + COHP side-by-side
bands_cohp_figure(bs, cohp_result, max_pairs=8)

CLI¶

# Pair table with ICOOP/-ICOHP values
vibeqc coop output.qvf

# Machine-readable JSON
vibeqc coop --json output.qvf

# Mayer bond order table
vibeqc mayer output.qvf
vibeqc mayer --json --threshold 0.1 output.qvf

Periodic Mayer bond orders¶

The k-space generalization of Mayer bond orders is available as a scalar companion to the energy-resolved COOP/COHP:

from vibeqc.coop_cohp import periodic_mayer_bond_orders

M = periodic_mayer_bond_orders(
    [T_real, V_real, F2e_real], S_real, system, basis, kmesh,
    n_electrons_per_cell=n_elec,
    F_terms_beta=fock_beta,  # optional, for unrestricted
)
# M.shape == (n_atoms, n_atoms), symmetric, B_AB on off-diagonals