COOP and COHP bonding analysis

Crystal Orbital Overlap Population (COOP) and Crystal Orbital Hamilton Population (COHP) resolve periodic bonding interactions by energy and atom pair. They complement a band or PDOS plot:

  • PDOS asks which AO groups contribute to a state.

  • COOP asks whether the overlap population between two atoms is bonding or antibonding.

  • COHP weights that interaction by a Hamiltonian matrix element.

vibe-qc returns and plots the Lobster-style negative COHP, written -COHP, so positive values indicate bonding in both the COOP and displayed COHP curves.

COOP and COHP are production on the GDF routes. The analysis projects the same spin Hamiltonian that produced the analysed eigenpairs, the accepted SCF Fock matrix per spin channel, and refuses any other operator: an H_terms whose Bloch sum differs from the analysed Fock raises ValueError instead of being labelled COHP (python/vibeqc/coop_cohp.py). On the BIPOLE, AICCM and full-grid external-XC routes coop_cohp=True is not computed; the request fails before SCF with the reason stated (NotImplementedError: the generic property path would rebuild a fixed-cutoff surrogate operator that is not the converged Hamiltonian of that route), and a QVF written on such a route records property_hamiltonian = "omitted (...)" in the .system manifest rather than shipping curves of a different operator. The label vocabulary is defined in Methods: theory and status.

Definitions and sign convention

For atom pair \(A,B\), band \(n\), and wave vector \(\mathbf{k}\), define the pair-resolved overlap and chosen-Hamiltonian populations

\[ P^{S}_{AB,n\mathbf{k}}= \operatorname{Re}\sum_{\mu\in A}\sum_{\nu\in B} S_{\mu\nu}(\mathbf{k}) C^*_{\mu n}(\mathbf{k})C_{\nu n}(\mathbf{k}), \]
\[ P^{H}_{AB,n\mathbf{k}}= \operatorname{Re}\sum_{\mu\in A}\sum_{\nu\in B} H_{\mu\nu}(\mathbf{k}) C^*_{\mu n}(\mathbf{k})C_{\nu n}(\mathbf{k}). \]

After Gaussian broadening,

\[ \mathrm{COOP}_{AB}(E)= \sum_\mathbf{k}w_\mathbf{k}\sum_n P^{S}_{AB,n\mathbf{k}}\,\delta_\sigma(E-\varepsilon_{n\mathbf{k}}), \]
\[ \mathrm{-COHP}_{AB}(E)= -\sum_\mathbf{k}w_\mathbf{k}\sum_n P^{H}_{AB,n\mathbf{k}}\,\delta_\sigma(E-\varepsilon_{n\mathbf{k}}). \]

For standard COHP, \(\mathbf{H}\) is the same spin Hamiltonian used to generate \(\mathbf{C}\) and \(\varepsilon\) (Dronskowski, Richard and Blöchl, Peter E. (1993). *Crystal orbital Hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations.* Journal of Physical Chemistry, **97**(33), 8617–8624. doi:[10.1021/j100135a014](https://doi.org/10.1021/j100135a014), Eqs. 4 to 12, which rewrite the band energy \(\sum_j f_j\varepsilon_j\) exactly as \(\int^{E_F} \mathrm{COHP}(E)\,\mathrm dE\); the identity holds only for the operator that produced the eigenpairs). A projection of any other operator must be identified as an explicit operator projection, rather than standard COHP. COOP is the overlap-weighted analogue of Hughbanks and Hoffmann (Hughbanks, Timothy and Hoffmann, Roald (1983). *Chains of trans-edge-sharing molybdenum octahedra: metal-metal bonding in extended systems.* Journal of the American Chemical Society, **105**(11), 3528–3537. doi:[10.1021/ja00349a027](https://doi.org/10.1021/ja00349a027)); the AO projection follows Deringer, Tchougreeff and Dronskowski (Deringer, Volker L., Tchougréeff, Andrei L., and Dronskowski, Richard (2011). *Crystal orbital Hamilton population (COHP) analysis as projected from plane-wave basis sets.* Journal of Physical Chemistry A, **115**(21), 5461–5466. doi:[10.1021/jp202489s](https://doi.org/10.1021/jp202489s)) and the sign and plotting convention that of LOBSTER (Maintz, Stefan, Deringer, Volker L., et al. (2016). *LOBSTER: A tool to extract chemical bonding from plane-wave based DFT.* Journal of Computational Chemistry, **37**(11), 1030–1035. doi:[10.1002/jcc.24300](https://doi.org/10.1002/jcc.24300)). The returned COHP array already includes the leading minus sign. Interpret the displayed convention as:

Curve value

COOP

Returned -COHP

positive

bonding overlap population

bonding Hamilton population

negative

antibonding overlap population

antibonding Hamilton population

The integrated values through the reference energy are ICOOP and -ICOHP. They summarize the occupied contribution for each pair. Always state the sign convention because some programs plot raw COHP instead.

Result shape

The lower-level COOPCOHPResult contains:

Attribute

Restricted shape

Unrestricted shape

Meaning

energies

(n_energy,)

(n_energy,)

Hartree energy grid

coop

(n_pairs, n_energy)

(2, n_pairs, n_energy)

alpha/beta split when unrestricted

cohp

same as coop, or None

same as coop, or None

returned -COHP; absent without H_terms

integrated_coop

(n_pairs,)

(2, n_pairs)

ICOOP through the reference energy

integrated_cohp

(n_pairs,), or None

(2, n_pairs), or None

returned -ICOHP

pairs

list

list

atom indices, symbols, and distance in Angstrom

The high-level workflow serializes these arrays into QVF rather than attaching real-space Fock terms to the generic SCF result object.

Pair selection

ao_pairs_per_atom_pair(system, basis, pair_distance_cutoff=8.0) creates AO index lists only for distinct atom pairs \(A<B\) whose raw home-cell coordinate distance is within the cutoff. The cutoff is in bohr. The high-level runner does not currently expose this cutoff.

This selector does not apply a minimum-image convention, enumerate image-resolved bonds, or include a bond from an atom to one of its own periodic images. A one-atom primitive cell therefore produces no pairs. Bloch-summed weights include lattice translations, but the pair label and distance still describe the raw home-cell atoms. Inspect that limitation before assigning a curve to a bond crossing the cell boundary.

Advanced explicit-block API

compute_coop_cohp is also available when a workflow already owns the real-space lattice blocks. Its positional arguments are Fock terms, overlap, system, basis, and k mesh. Set include_cohp=True to project the analyzed Hamiltonian in each spin channel:

import vibeqc as vq


def analyze_lattice_blocks(
    fock_terms,
    overlap_real,
    system,
    basis,
    kmesh,
):
    return vq.compute_coop_cohp(
        fock_terms,
        overlap_real,
        system,
        basis,
        kmesh,
        include_cohp=True,
        pair_distance_cutoff=8.0,
        sigma=0.01,
        n_electrons_per_cell=system.n_electrons(),
    )

The compatibility H_terms input also enables COHP, but is accepted only when its Bloch sum equals the analyzed Hamiltonian for every spin. A different operator is refused. The low-level cohp_weights_k kernel remains available for explicitly named arbitrary-operator projections.

For an unrestricted operator, pass the beta-spin lattice blocks as F_terms_beta=. Do not copy older examples that access result.fock_terms, result.overlap_real, or result.hcore_terms: generic periodic results do not promise those fields.

periodic_mayer_bond_orders uses the same explicit-block model for a scalar bond-order matrix; the next section gives its definition. The production runner computes its own consistent version for QVF output.

Periodic Mayer bond orders

The Mayer bond order (Mayer, I. (1983). *Charge, bond order and valence in the ab initio SCF theory.* Chemical Physics Letters, **97**(3), 270–274. doi:[10.1016/0009-2614(83)80005-0](https://doi.org/10.1016/0009-2614(83)80005-0)) is the energy-integrated, scalar companion of COOP: where COOP resolves the overlap population between two atoms by energy, the Mayer index sums the occupied contribution into one number per atom pair. In a periodic system the molecular definition is averaged over the k mesh,

\[ B_{AB} = \sum_{\mathbf k} w_{\mathbf k}\sum_{\mu\in A}\sum_{\nu\in B} \big[\mathbf P(\mathbf k)\,\mathbf S(\mathbf k)\big]_{\mu\nu}\, \big[\mathbf P(\mathbf k)\,\mathbf S(\mathbf k)\big]_{\nu\mu}, \qquad \mathbf P(\mathbf k) = \mathbf C(\mathbf k)\,\mathbf n\,\mathbf C(\mathbf k)^\dagger, \]

with \(\mathbf n\) the diagonal occupation matrix (2 per occupied band in the restricted case). For an unrestricted density the \(\alpha\) and \(\beta\) contributions are formed with occupation 1 each and summed with the standard molecular factor of 2, so a closed shell reproduces the restricted value. The result is a symmetric (n_atoms, n_atoms) matrix with zero diagonal; it depends on the basis and on the Mulliken-type partition exactly as the molecular index does, so compare values only across runs with the same basis and route.

periodic_mayer_bond_orders(F_terms, S_real, system, basis, kmesh, *, n_electrons_per_cell=None, F_terms_beta=None) in vibeqc.coop_cohp is the explicit-block API: it Bloch-sums the Fock terms and the overlap at every mesh point, diagonalises, fills the lowest n_electrons_per_cell / 2 bands (or the \(\alpha\) and \(\beta\) counts implied by the cell multiplicity when F_terms_beta is given) and evaluates the formula above. With n_electrons_per_cell=None every band is occupied, the non-interacting Hcore limit, which is only useful as a shape check.

The runner computes the same matrix on every run that produces the DOS payload, on the GDF routes from the accepted SCF eigenstates, and writes the pairs with \(B_{AB} \ge 0.05\) into the QVF bond_orders section as a JSON member bond_orders with method = "mayer" and a pairs list; each pair records i, j, symbol_i, symbol_j, order and distance_ang (the raw home-cell distance, as for the COOP pairs). The CLI reads that section back:

vibeqc mayer h2-cube-coop-cohp.qvf                    # table, pairs with order >= 0.05
vibeqc mayer --threshold 0.1 h2-cube-coop-cohp.qvf    # raise the display cutoff
vibeqc mayer --json h2-cube-coop-cohp.qvf             # every stored pair, machine readable

For the H2 cube above the table holds one row, H0-H1  0.74 Å  1.0000: a single bond shared over the eight k points. A QVF without a bond_orders section (a route whose property payload was omitted, or a run without output_qvf=True) makes the command exit with contains no bond_orders section rather than printing zeros. For molecular Mayer and Wiberg indices, and the periodic Wiberg variant, see bond analysis.

Plotting an explicit result

The plotters accept a COOPCOHPResult directly:

from vibeqc.plot import cohp_figure, coop_figure


def save_pair_plots(analysis):
    coop = coop_figure(analysis, max_pairs=8, title="COOP")
    coop.savefig("coop.png", dpi=180)

    cohp = cohp_figure(analysis, max_pairs=8, title="-COHP")
    cohp.savefig("cohp.png", dpi=180)

bands_cohp_figure(bands, analysis) produces a combined view when bands is a valid BandStructure computed from the same operator and energy convention. The helper shifts each panel by its own stored reference, so first verify that those references agree. Otherwise align both objects to one stated reference and use shift_to_fermi=False. Do not combine Hcore bands with SCF-derived COHP and present them as one spectrum.

Convergence and interpretation

COOP and COHP inherit every approximation in the underlying SCF and add their own reconstruction, mesh, broadening, and pair-selection choices. Use this order:

  1. Converge the SCF basis, Coulomb route, k mesh, and numerical thresholds.

  2. Converge dos_kmesh for the pair curves and integrated values.

  3. Hold the mesh fixed while choosing Gaussian broadening for presentation.

  4. Inspect the available home-cell pairs and exclude claims that require an image-resolved or same-atom image bond.

  5. Compare atom-pair trends across at least one reasonable basis change.

  6. Retain the QVF and all provenance/citation sidecars.

An integrated value is not a universal bond energy. It depends on the basis, Hamiltonian, occupations, reference energy, pair partition, and sign convention. Use it for controlled comparisons with all of those choices held fixed.

For metals, the k mesh, occupations, electronic temperature, entropy term, and energy/free-energy convention must be converged together. See Smearing.

CPU time and memory

At each of \(N_\mathbf{k}\) k points, the analysis diagonalizes an \(N_\mathrm{bf}\times N_\mathrm{bf}\) operator and evaluates the requested pair weights. The leading diagonalization trend is

\[ T_\mathrm{diag}\sim\mathcal{O}(N_\mathbf{k}N_\mathrm{bf}^3). \]

Pair accumulation adds work proportional to the number of selected pairs and energy-grid points. In addition to \(\mathcal{O}(N_\mathrm{bf}^2)\) matrix workspace, the current implementation forms transient pair weights of shape \((N_\mathrm{pair},N_\mathbf{k},N_\mathrm{bf})\) per spin. Stored curves scale as \(\mathcal{O}(N_\mathrm{spin}N_\mathrm{pair}N_E)\).

Record CPU time, wall time, peak RSS, basis count, SCF and analysis meshes, pair-selection rule, number of pairs, energy-grid size, sigma, spin treatment, threads, host, and version.