Post-SCF properties

Every run_job calculation now reports a properties block at the end of the .out file with atomic charges, bond orders, and the dipole moment — the standard sanity-check output every QC program emits.

What you’ll see in the .out

Atomic charges
----------------------------------------------------
   #  elem    Mulliken      Löwdin
----------------------------------------------------
   1     O   -0.910152   -0.646928
   2     H    0.455076    0.323464
   3     H    0.455076    0.323464
----------------------------------------------------
        sum  -0.000000

Bond orders (Mayer)
----------------------------------------
        atoms    bond order
----------------------------------------
O1   — H2          0.7628
O1   — H3          0.7628

Dipole moment (origin: centre of mass)
----------------------------------------------------
   component     a.u. (e·bohr)       Debye
----------------------------------------------------
           x       -0.00000000     -0.0000
           y       -0.00000000     -0.0000
           z       -0.81069504     -2.0606
         |μ|        0.81069504      2.0606

Bond orders below 0.10 are suppressed from the table to keep the output compact; the full matrix is available programmatically.

Programmatic access

import vibeqc as vq

mol = vq.Molecule.from_xyz("h2o.xyz")
basis = vq.BasisSet(mol, "6-31g*")
result = vq.run_rhf(mol, basis)

# Per-atom charges (numpy array, length = n_atoms)
mul = vq.mulliken_charges(result, basis, mol)
low = vq.loewdin_charges(result, basis, mol)

# Bond-order matrix (n_atoms × n_atoms numpy)
B = vq.mayer_bond_orders(result, basis, mol)
print("O–H bond order:", B[0, 1])

# Dipole moment (with conversions)
dip = vq.dipole_moment(result, basis, mol)
print(f"|μ| = {dip.total_debye:.3f} Debye  (origin: {dip.origin})")
print(f"μ = ({dip.x:.4f}, {dip.y:.4f}, {dip.z:.4f}) a.u.")

All four methods accept the result from any of the four molecular SCF drivers (run_rhf, run_uhf, run_rks, run_uks). UHF / UKS results are handled correctly — the total density D_α + D_β is used for populations and dipole, and Mayer bond orders use the spin-resolved formula 2·[(P_α S)(P_α S) + (P_β S)(P_β S)] which reduces to the closed-shell expression when D_α = D_β.

Validation

vibe-qc’s Mulliken charges and dipole-moment components match PySCF’s to better than 1e-6 on H₂O / 6-31G* (verified in tests/test_properties.py). The Mayer bond order on H₂ at the equilibrium bond length comes out near 1.0 as expected for a single covalent bond; in H₂O at HF/6-31G* the O–H bond orders are about 0.76 and the H…H “bond order” is below 0.005.

Caveats

  • Mulliken charges depend strongly on the basis set, especially on diffuse functions. Use them for trends along a series, not for quantitative charge assignment. Löwdin charges are usually less basis-sensitive.

  • The default origin for the dipole is the centre of mass, computed from a built-in atomic-mass table covering H–Kr. For heavier elements pass an explicit origin= (a 3-vector in bohr).

  • Dipole values for charged species are origin-dependent; callers should choose an origin that matches the convention they want to report (centre of charge, nuclear centroid, etc.).

  • Periodic systems are not yet covered by the properties block. run_rhf_periodic / run_rks_periodic results don’t currently flow through the .out output infrastructure; periodic dipole moments themselves are subtle anyway (require a Berry-phase polarisation treatment for unambiguous values), and arrive in a later phase.