Bond analysis: Wiberg, NBO, EDA, IBOs, and orbital entanglement

Chemists rarely stop at a total energy. The questions that follow a converged SCF – how strong is this bond, where did the charge go, what holds these fragments together – are answered by bond-analysis methods, and vibe-qc ships four complementary families of them. Each family looks at the same density matrix through a different lens: Wiberg indices count shared electron pairs, the provisional, exploratory NBO helper classifies orbitals, EDA splits an interaction energy into physical components, and entanglement measures read bonding straight off the quantum-information structure of the wavefunction. Bond orders, EDA and the entanglement measures are production on HF and DFT densities; NPA is gated and the NBO search provisional (see Methods: theory and status).

For energy-resolved periodic bonding analysis (COHP/COOP curves, band-by-band bonding character), use the canonical COOP/COHP module instead – the methods on this page give integrated per-bond scalars.

What ships where

Wiberg bond indices are computed automatically on every population dump. The .population.{txt,json} sidecars retain an NPA field, but it is empty and reports unavailable.npa until the full occupancy-weighted Natural Atomic Orbital construction is implemented. This prevents a Löwdin population from being reported under the NPA label. The NBO search is a provisional, exploratory Python helper; it, the EDA and the entanglement analyses are Python-API workflows you call on an SCF result.

Theory

All four families read the same converged one-particle density matrix \(\mathbf P\) (spin-summed, \(\mathbf P = \mathbf P^\alpha + \mathbf P^\beta\) for open shells) and the AO overlap \(\mathbf S\); what distinguishes them is which contraction they take and in which basis.

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)) contracts through the overlap in the raw AO basis,

\[ B_{AB} = \sum_{\mu\in A}\sum_{\nu\in B} (\mathbf P\mathbf S)_{\mu\nu}\,(\mathbf P\mathbf S)_{\nu\mu} \]

for a closed shell (open shells use \(2[(\mathbf P^\alpha\mathbf S)(\mathbf P^\alpha\mathbf S) + (\mathbf P^\beta\mathbf S)(\mathbf P^\beta\mathbf S)]\) summed the same way), while the Wiberg index (Wiberg, Kenneth B. (1968). *Application of the Pople-Santry-Segal CNDO method to the cyclopropylcarbinyl and cyclobutyl cation and to bicyclobutane.* Tetrahedron, **24**(3), 1083–1096. doi:[10.1016/0040-4020(68)88057-3](https://doi.org/10.1016/0040-4020(68)88057-3)) squares the elements of the density in the symmetrically orthogonalised (Loewdin) basis, \(\tilde{\mathbf P} = \mathbf S^{1/2}\mathbf P\mathbf S^{1/2}\):

\[ W_{AB} = \sum_{\mu\in A}\sum_{\nu\in B} \tilde P_{\mu\nu}^{\,2} . \]

Both count shared electron pairs, both give 1 for H2 in a minimal basis, and both depend on the basis because “atom A” means “the functions centred on A”; the Loewdin transformation makes the Wiberg index the less sensitive of the two. The delocalisation index of QTAIM (Matito, Eduard, Solà, Miquel, et al. (2007). *Electron sharing indexes at the correlated level. Application to aromaticity calculations.* Faraday Discussions, **135**, 325–345. doi:[10.1039/B605086G](https://doi.org/10.1039/B605086G), Outeiral, Carlos, Vincent, Mark A., et al. (2018). *Revitalizing the concept of bond order through delocalization measures in real space.* Chemical Science, **9**, 5517–5529. doi:[10.1039/C8SC01338A](https://doi.org/10.1039/C8SC01338A)) is defined without a basis, as the double integral of the exchange-correlation density over two atomic basins, which for a single determinant reduces to basin overlaps of the occupied orbitals:

\[ \delta(A,B) = -2\int_{\Omega_A}\!\!\int_{\Omega_B}\rho_{\mathrm{xc}}(\mathbf r_1,\mathbf r_2)\,d\mathbf r_1\,d\mathbf r_2 \;=\; 2\sum_{ij} S_{ij}(\Omega_A)\,S_{ij}(\Omega_B) . \]

vibe-qc has no basin integration yet (QTAIM guide), so delocalization_index replaces the basin overlaps by an AO partition, \(\delta^{\mathrm{AO}}_{AB} = 2\sum_{\mu\in A,\nu\in B}(\mathbf P\mathbf S)_{\mu\nu}(\mathbf P\mathbf S)_{\nu\mu} - \sum_{\mu\in A,\nu\in B}\tilde P_{\mu\nu}^{\,2}\). That is an approximation with no positivity guarantee: on H2O/STO-3G it returns \(\delta(\mathrm O,\mathrm H) = -0.17\) where the Wiberg index is 0.98 and the Mayer order 0.94 (vibe-qc 0.17.x development build), so read it as a diagnostic beside the two bond orders, not as a QTAIM number.

Energy decomposition. LMO-EDA (Su, Peifeng and Li, Hui (2009). *Energy decomposition analysis of covalent bonds and intermolecular interactions.* Journal of Chemical Physics, **131**(1), 014102. doi:[10.1063/1.3159673](https://doi.org/10.1063/1.3159673)) and the Kitaura-Morokuma scheme (Morokuma, Keiji (1971). *Molecular orbital studies of hydrogen bonds. III. C=O…H-O hydrogen bond in H2CO…H2O and H2CO…2H2O.* Journal of Chemical Physics, **55**(3), 1236–1244. doi:[10.1063/1.1676210](https://doi.org/10.1063/1.1676210)) write the interaction energy of two fragments as a sum of steps from the isolated fragments, through the frozen promolecule and the frozen orthogonalised state, to the relaxed dimer,

\[ \Delta E_{\mathrm{int}} = E_{AB} - E_A - E_B = \Delta E_{\mathrm{elstat}} + \Delta E_{\mathrm{exch}} + \Delta E_{\mathrm{rep}} + \Delta E_{\mathrm{pol+CT}} + \Delta E_{\mathrm{disp}} . \]

The terms are built from the fragment densities and the AO integrals and sum to \(\Delta E_{\mathrm{int}}\) exactly; the EDA section below defines each of them, gives the call signature and shows a counterpoise water-dimer example.

Orbital entanglement. Tracing a normalised wavefunction over all but one spatial orbital leaves a \(4\times4\) one-orbital reduced density matrix in the Fock-space basis \(\{|0\rangle, |\!\uparrow\rangle, |\!\downarrow\rangle, |\!\uparrow\downarrow\rangle\}\) with eigenvalues \(\omega_{i,\alpha}\); its von Neumann entropy and the mutual information of two orbitals (Legeza, Ors and Solyom, Jeno (2003). *Optimizing the density-matrix renormalization group method using quantum information entropy.* Physical Review B, **68**(19), 195116. doi:[10.1103/PhysRevB.68.195116](https://doi.org/10.1103/PhysRevB.68.195116), Szalay, Szilárd, Barcza, Gergely, et al. (2017). *The correlation theory of the chemical bond.* Scientific Reports, **7**, 2237. doi:[10.1038/s41598-017-02447-z](https://doi.org/10.1038/s41598-017-02447-z)) are

\[ s_i = -\sum_\alpha \omega_{i,\alpha}\ln\omega_{i,\alpha}, \qquad I_{ij} = s_i + s_j - s_{ij}, \]

with \(s_{ij}\) the entropy of the \(16\times16\) two-orbital density. \(s_i\) vanishes for an orbital that is empty or doubly occupied in every determinant and reaches \(\ln 4\) at maximal entanglement, so it measures how strongly an orbital takes part in correlation; \(I_{ij}\) is large for the pairs that must be correlated together, which makes it an active-space selector and, in the maximally entangled atomic-orbital picture (Ding, Lexin, Matito, Eduard, and Schilling, Christian (2026). *Chemical bonding concepts emerge naturally from maximally entangled atomic orbitals.* Nature Communications, **17**(1), 4732. doi:[10.1038/s41467-026-73527-w](https://doi.org/10.1038/s41467-026-73527-w)), a bonding index. For a single determinant the one-orbital density is diagonal in the natural occupations \(n_i\) and \(s_i = -n_i\ln n_i - (1-n_i)\ln(1-n_i)\), which is what entanglement_from_density evaluates from the eigenvalues of the orthogonalised density, estimating \(I_{ij}\) from products of fractional occupations. That shortcut is only informative for densities with fractional natural occupations (DFT, broken-symmetry UHF, correlated 1-RDMs); for a correlated wavefunction feed genuine one- and two-orbital RDMs to single_orbital_entropy and mutual_information. See the caveat in the entanglement section below about which orthogonalisation the shortcut applies.

Intrinsic bond orbitals (Knizia, Gerald (2013). *Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory and Chemical Concepts.* Journal of Chemical Theory and Computation, **9**(11), 4834–4843. doi:[10.1021/ct400687b](https://doi.org/10.1021/ct400687b)) and the natural-population gate (Reed, Alan E., Weinstock, Robert B., and Weinhold, Frank (1985). *Natural population analysis.* Journal of Chemical Physics, **83**(2), 735–746. doi:[10.1063/1.449486](https://doi.org/10.1063/1.449486)) are described in their own sections below.

Wiberg bond indices and the delocalization index

The Wiberg bond index counts the shared electron pairs between two atoms from the squared density-matrix elements in the Löwdin (symmetrically orthogonalised) basis. It is less basis-set sensitive than the Mayer bond order and complements it: Mayer contracts through the overlap matrix, Wiberg through the orthogonalised density.

import vibeqc as vq
from vibeqc.bond_analysis import wiberg_bond_orders, bond_order_summary

mol = vq.Molecule([
    vq.Atom(8, [0.0, 0.0, 0.0]),
    vq.Atom(1, [0.0, 1.43, -0.98]),
    vq.Atom(1, [0.0, -1.43, -0.98]),
])
basis = vq.BasisSet(mol, "sto-3g")
result = vq.run_rhf(mol, basis)

W = wiberg_bond_orders(result, basis, mol)   # (n_atoms, n_atoms)
print("O-H Wiberg index:", W[0, 1])

# All bond-order metrics side by side:
summary = bond_order_summary(
    result, basis, mol,
    compute_mayer=True,
    compute_wiberg=True,
    compute_delocalization=True,
)
print(summary.summary())

delocalization_index computes an AO-approximated delocalization index (DI) from the Löwdin-basis density; the exact DI is defined over QTAIM atomic basins (see the QTAIM guide). Periodic Gamma-point results go through the same code via vibeqc.bond_analysis.periodic_wiberg_bond_orders and periodic_delocalization_index.

References: Wiberg, Tetrahedron 24, 1083 (1968); Matito, Solà, Salvador & Duran, Faraday Discuss. 135, 325 (2007); Outeiral, Vincent, Martín Pendás & Popelier, Chem. Sci. 9, 5517 (2018).

Energy decomposition analysis (EDA)

Status: preview. The terms are real and add up to the interaction energy, for two closed-shell Hartree-Fock fragments only; other inputs are refused. The Kohn-Sham form, open-shell fragments and a run_job keyword are not implemented. The label vocabulary is in Methods: theory and status.

EDA answers “why do these two fragments bind” by partitioning the interaction energy of two closed-shell fragments into electrostatic, exchange, Pauli-repulsion, polarisation and dispersion terms. vibeqc.eda.eda_lmo implements the LMO-EDA partition of Su and Li (2009) and vibeqc.eda.eda_morokuma the Kitaura-Morokuma partition (1971, 1976) of the same states, both as post-processing of three converged Hartree-Fock runs: the dimer and the two fragments, each fragment computed in the dimer basis so that the interaction energy is counterpoise corrected. vibe-qc has no ghost-atom keyword, but a Molecule accepts atoms with Z = 0, and the explicit-shell BasisSet constructor puts the dimer’s shells on them. The default SAD guess has no atomic density for Z = 0, so the fragment runs start from the core Hamiltonian.

import vibeqc as vq
from vibeqc.eda import eda_lmo

ANG = 1.8897259886   # Angstrom to bohr
xyz = [  # water dimer, R(O..O) = 2.95 A
    (8, (0.000, 0.000, 0.000)), (1, (0.000, 0.000, 0.960)), (1, (0.000, 0.929, -0.240)),
    (8, (0.000, 0.000, 2.950)), (1, (0.000, 0.759, 3.538)), (1, (0.000, -0.759, 3.538)),
]
frag1, frag2 = [0, 1, 2], [3, 4, 5]

def molecule(real):
    """The dimer with every atom outside `real` turned into a ghost (Z = 0)."""
    return vq.Molecule([vq.Atom(z if i in real else 0, [c * ANG for c in p])
                        for i, (z, p) in enumerate(xyz)], 0, 1)

dimer = molecule(frag1 + frag2)
basis = vq.BasisSet(dimer, "6-31g*")
shells = basis.shells()

def scf(mol):
    opts = vq.RHFOptions()
    if any(a.Z == 0 for a in mol.atoms):
        opts.initial_guess = vq.InitialGuess.HCORE      # SAD has no density for Z = 0
        return vq.run_rhf(mol, vq.BasisSet(mol, shells, name="6-31g*"), opts)
    return vq.run_rhf(mol, basis, opts)

ab, a, b = scf(dimer), scf(molecule(frag1)), scf(molecule(frag2))
eda = eda_lmo(
    ab.energy, a.energy, b.energy,
    None, None, None,                        # Fock matrices: accepted, not used
    ab.density, a.density, b.density,
    vq.compute_overlap(basis),
    basis=basis, molecule=dimer, frag1_atoms=frag1, frag2_atoms=frag2,
)
print(eda.summary())

Energy Decomposition Analysis (LMO-EDA)
----------------------------------------------------------------
  E_total (supersystem)   =  -152.02576771 E_h
  E_fragment_1 (isolated) =   -76.00957528 E_h
  E_fragment_2 (isolated) =   -76.00960659 E_h
----------------------------------------------------------------
  dE_int (total)          =    -0.00658583 E_h      -4.13 kcal/mol

  Electrostatic           =    -0.01255786 E_h      -7.88 kcal/mol
  Exchange                =    -0.01136365 E_h      -7.13 kcal/mol
  Repulsion (Pauli)       =     0.01991724 E_h      12.50 kcal/mol
  Polarisation (+CT)      =    -0.00258156 E_h      -1.62 kcal/mol
----------------------------------------------------------------

The partition follows one path: the two fragments, the promolecule (P_1 + P_2, the fragment orbitals as they are), the frozen state (one determinant of the mutually orthogonalised occupied orbitals), and the relaxed dimer.

Term

Definition

e_elstat

Classical Coulomb energy of the two frozen charge distributions, nuclei included.

e_exch

Inter-fragment exchange energy of the promolecule state (LMO-EDA).

e_rep

Pauli repulsion beyond exchange: the frozen state minus the promolecule.

e_pol

Converged dimer minus frozen state: polarisation and charge transfer together.

e_disp

Difference of additive dispersion corrections, if you pass them (e_disp, e_disp_frag1, e_disp_frag2).

The five terms add up to e_int = e_total - e_frag1 - e_frag2 exactly. eda_morokuma reports the Kitaura-Morokuma terms of the same states: e_exch is then the full exchange repulsion EX (the sum of the LMO-EDA exchange and repulsion, +5.37 kcal/mol here), e_rep is zero, and e_pol is polarisation, charge transfer and coupling as one number. The restricted SCF variations that separate those three are not implemented.

What to expect: the electrostatic attraction of a hydrogen bond is larger than the binding energy, Pauli repulsion cancels most of it, and relaxation adds a smaller stabilisation. Only the sum is an observable; the split depends on the chosen intermediate states, and the polarisation and charge-transfer terms are not basis-set stable. Compare them within one basis.

Scope, enforced by the code. The terms are functionals of the fragment densities, so eda_lmo needs the integrals and takes basis, molecule and the atom indices of each fragment beside the SCF results. It evaluates the Hartree-Fock energy expression on the supplied densities and refuses, with a ValueError, energies that do not match it to consistency_tol (1e-6 Hartree): Kohn-Sham energies, ECP runs and fragments computed in another basis all fail here rather than return wrong numbers. It also refuses densities that are not closed-shell SCF densities. That excludes fragment_density_matrix, which returns a projection of the dimer orbitals and not a fragment density; use the density of a fragment run in the dimer basis. Open-shell fragments and the Kohn-Sham form of LMO-EDA are not implemented. The Fock-matrix arguments of the older signature are still accepted and are ignored.

References: Su & Li, J. Chem. Phys. 131, 014102 (2009); Morokuma, J. Chem. Phys. 55, 1236 (1971); Kitaura & Morokuma, Int. J. Quantum Chem. 10, 325 (1976).

Orbital entanglement measures

Quantum-information analysis reads bonding from the entanglement structure of the wavefunction: the single-orbital entropy measures how strongly an orbital participates in correlation, and the mutual information between two orbitals quantifies their bonding entanglement. vibe-qc offers an occupation-number approximation directly from a spin-restricted, spin-summed AO density matrix:

from vibeqc.entanglement import entanglement_from_density

ent = entanglement_from_density(P, S)
print("total quantum information:", ent.total_information)
print("correlation clusters:", ent.correlation_clusters)

The density is transformed as \(S^{1/2} P S^{1/2}\), using the same Loewdin convention as Wiberg bond orders. Its eigenvalues are the spin-summed natural occupations in \([0, 2]\). The binary entropy uses per-spin occupations \(n_i\), obtained by dividing those eigenvalues by two: \(s_i = -n_i\ln n_i - (1-n_i)\ln(1-n_i)\). Closed-shell RHF densities give zero entropy, mutual information and entanglement bond orders; fractional occupations, such as those of stretched H2 CASCI, give positive occupation entropy. Endpoint deviations within \(10^{-14}\) per spin are treated as eigensolver roundoff.

Warning

entanglement_from_density reports an occupation-based approximation. Its mutual information and bond orders are heuristic estimates; a one-particle density alone does not determine full orbital entanglement. Zero entropy for a closed-shell determinant refers to its natural orbitals, not to an arbitrary choice of localized orbitals. Supply the spin-summed density with Tr(P @ S) equal to the electron count.

For multi-determinantal wavefunctions (CAS, CC densities), feed the one- and two-orbital reduced density matrices to vibeqc.entanglement.single_orbital_entropy and vibeqc.entanglement.mutual_information directly.

References: Legeza & Sólyom, Phys. Rev. B 68, 195116 (2003); Szalay, Barcza, Szilvási, Veis & Legeza, Sci. Rep. 7, 2237 (2017); Ding, Matito & Schilling, Nat. Commun. 17, 4732 (2026), DOI 10.1038/s41467-026-73527-w.

Intrinsic bond orbitals (IBOs) and IAO charges

The methods above all read scalars off the density matrix. IBOs answer the complementary question – what do the bonds look like – by rotating the occupied orbitals into a set that reproduces the Lewis structure: one orbital per bond, core, or lone pair, with no Lewis pattern assumed anywhere in the construction.

import vibeqc as vq

# On by default for closed-shell molecular jobs: writes one localized
# wavefunction section per criterion for vibe-view to switch between.
res = vq.run_job(molecule=mol, basis="def2-svp", method="rhf")

res = vq.run_job(..., localize=False)     # off
res = vq.run_job(..., localize="ibo")     # just one criterion
res = vq.run_job(..., localize=["ibo", "boys"])

Three criteria, because they disagree

The default emits IBO, Foster-Boys and Pipek-Mezey. That is not redundancy: the criteria answer the same question differently, and on benzene they visibly part company.

benzene/def2-SVP

1-centre

2-centre

3-centre

IBO

6

12

3 (the pi system)

Foster-Boys

6

15

0

Pipek-Mezey

6

12

3

Boys maximises the orbital dipole spread and so mixes sigma with pi, turning the aromatic pi system into three more two-centre “banana” bonds. IBO and Pipek-Mezey maximise atomic populations and keep sigma and pi apart. Neither picture is wrong; which one you want depends on the question. Comparing them is the point.

Whichever criterion runs, the descriptors – atomic populations, centre counts, charges – are computed on the IAO yardstick, so the numbers are directly comparable across the three. The IAO charges are identical for all three by construction: localization is a unitary rotation inside the occupied space, and the charges are a property of that space.

Or directly, on a converged occupied set:

from vibeqc.iao import analyse_ibo

analysis = analyse_ibo(mol, basis, C_occ)
print(analysis.charges)           # IAO partial charges, one per atom
print(analysis.n_centres)         # 1 = core/lone pair, 2 = bond, 3+ = delocalised
print(analysis.atom_populations)  # n_A(i); each row sums to 1
print(analysis.centroids)         # <i|r|i> in bohr

For methane this gives one carbon core plus four equivalent C-H bonds; for benzene, six cores, twelve two-centre sigma bonds, and a pi system that genuinely cannot be squeezed onto two centres.

IAO charges instead of Mulliken. The partial charges that come with this are basis-set stable in a way Mulliken charges are not – the CH4 carbon moves only from -0.47 to -0.50 between def2-SVP and cc-pVTZ. That stability is the reason IBOs are built on IAO populations rather than the Mulliken populations that Pipek-Mezey uses, and it is why vibeqc.localise.pipek_mezey_localise is not the right tool for interpretation even though it is available (it exists for DLPNO domain construction, where the criterion is cost, not physical meaning).

The difference is not subtle. On the same CH4/def2-SVP wavefunction, Mulliken puts -0.144 on carbon where IAO puts -0.473. When localize="ibo" is set, the IAO charges are written into the QVF’s atom_properties section alongside the Mulliken, Löwdin and Hirshfeld rows, so vibe-view’s charge picker offers them directly.

Use iao_analysis=True to request IAO charges, spin populations and spin-resolved IAO-Wiberg bonds independently of localization or QVF. The requested section is then included in .out and the population text/JSON sidecars. See IAO population analysis for its definition, API, support limits and periodic integration plan.

Scope and cost

Closed-shell molecular SCF only. Measured against the SCF that produced them, all three criteria together cost about +16 % on a typical small job (H2O/def2-TZVPP) and +62 % on benzene/def2-SVP. Nearly all of the aromatic case is Foster-Boys, whose p=2 functional converges slowly on degenerate pi systems – ~89 Jacobi sweeps against IBO’s 14. If that matters for a batch, localize="ibo" keeps the cheapest and most interpretable one.

Jobs using ECPs are refused – an all-electron minimal reference basis cannot partition a valence-only target basis – as are elements past Z=86. In each case the job completes and warns; only the localized section is omitted.

Reference: Knizia, J. Chem. Theory Comput. 9, 4834 (2013), doi:10.1021/ct400687b. The minimal reference basis is Huzinaga MINI (Huzinaga et al., Gaussian Basis Sets for Molecular Calculations, Elsevier 1984).

Citations

Every method on this page carries its defining papers in the citation database. The DOIs are pinned mechanically by tests/test_citations.py::test_bond_analysis_routes_pin_dois. See Citations for the full provenance surface.