Post-SCF properties: charges, bonds, dipole

Once an SCF is converged, vibe-qc gives you the standard chemical interpretation tools: atomic charges (Mulliken and Löwdin), Mayer bond orders, and the molecular dipole moment. This tutorial walks through them on water, then a short comparison on H-bonded systems.

These properties are emitted automatically by run_job when it has enough information; you can also call the helpers directly to get the raw arrays for plotting or downstream analysis.

A single call

from vibeqc import Atom, Molecule, run_job

mol = Molecule([
    Atom(8, [ 0.0,  0.00,  0.00]),
    Atom(1, [ 0.0,  1.43, -0.98]),
    Atom(1, [ 0.0, -1.43, -0.98]),
])
run_job(mol, basis="6-31g*", method="rhf", output="water")

(Positions are in bohr. Molecule.from_xyz("examples/h2o.xyz") works too if you have an XYZ file.)

Open water.out — after the SCF trace and orbital-energy table you’ll find (HF/6-31G*, experimental geometry):

  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

What the numbers mean

  • Mulliken charges partition the electron density according to which basis function it was expanded in. Simple, interpretable, but basis-set dependent — a polarisation function added to H shifts charge onto it.

  • Löwdin charges partition in the orthonormalised AO basis (S^{-1/2} transform). Less basis-set sensitive, generally preferred for trend analysis across systems.

  • Mayer bond orders generalise the Wiberg index to non-orthogonal bases. B_ij ≈ 1 for single bonds, ≈ 2 for double, ≈ 3 for triple. Values under ~0.1 are hidden by default.

  • Dipole moment is computed with the origin at the centre of mass, which makes it origin-invariant for neutral molecules. For ions pass an explicit origin= to dipole_moment() if you want a different reference.

Pulling the arrays yourself

When you need raw data for plotting or fitting, call the helpers directly — they take the same result object you get back from run_rhf / run_job.

import numpy as np
from vibeqc import (
    Atom, BasisSet, Molecule, run_rhf,
    mulliken_charges, loewdin_charges,
    mayer_bond_orders, dipole_moment,
)
from vibeqc.properties import prominent_bonds

mol = Molecule([
    Atom(8, [ 0.0,  0.00,  0.00]),
    Atom(1, [ 0.0,  1.43, -0.98]),
    Atom(1, [ 0.0, -1.43, -0.98]),
])
basis = BasisSet(mol, "6-31g*")
result = run_rhf(mol, basis)

q_mull = mulliken_charges(result, basis, mol)      # (n_atoms,) np.ndarray
q_lowd = loewdin_charges(result, basis, mol)       # same shape

bonds = mayer_bond_orders(result, basis, mol)      # (n_atoms, n_atoms)
print(prominent_bonds(bonds, mol, threshold=0.10)) # [(i, j, B_ij), ...]

mu = dipole_moment(result, basis, mol)
print(f"|μ| = {mu.total_debye:.3f} D")

Basis-set dependence: a quick sweep

Mulliken charges drift significantly between basis sets — run the same water molecule at three levels and print the oxygen charge:

mol = Molecule([
    Atom(8, [ 0.0,  0.00,  0.00]),
    Atom(1, [ 0.0,  1.43, -0.98]),
    Atom(1, [ 0.0, -1.43, -0.98]),
])

for basis_name in ("sto-3g", "6-31g*", "def2-tzvp"):
    basis = BasisSet(mol, basis_name)
    r = run_rhf(mol, basis)
    q = mulliken_charges(r, basis, mol)
    print(f"{basis_name:12s}  q_O = {q[0]:+.3f}  μ (D) = "
          f"{dipole_moment(r, basis, mol).total_debye:.3f}")

Output (HF, experimental geometry):

sto-3g        q_O = -0.418  μ (D) = 1.725
6-31g*        q_O = -0.910  μ (D) = 2.061
def2-tzvp     q_O = -0.685  μ (D) = 2.008

Two lessons:

  1. Mulliken is sensitive to basis. q_O swings half an electron across these sets — never compare absolute Mulliken charges between different bases, different codes, or different molecules unless you hold the basis fixed. Löwdin charges are more stable but not bulletproof either.

  2. Dipole moment is well-behaved. 6-31g* already overshoots slightly; TZVP pulls it back toward 2.0 D. The experimental gas-phase water dipole is 1.85 D. The residual over-estimate comes from missing electron correlation — rerun the sweep with method="rks", functional="PBE" and the dipole drops further.

Following the H-bond in a water dimer

Run the dimer example (examples/input-h2o-dimer-opt.py) and load output-h2o-dimer-opt.molden in Avogadro. The .out file already lists the Mulliken charges at the optimised geometry — note the donor-hydrogen charge is noticeably larger than the free hydrogens, and the oxygen of the donor water is slightly more negative than the acceptor. That’s the H-bond polarisation, visible in the numbers.

Same in Python:

from vibeqc import Molecule, BasisSet, run_rhf, mulliken_charges
mol = Molecule.from_xyz("h2o_dimer_optimised.xyz")   # after optimisation
basis = BasisSet(mol, "6-31g*")
r = run_rhf(mol, basis)
q = mulliken_charges(r, basis, mol)
print("Charges on:", q)

Next

  • Dipole moment tracks the molecular geometry; see geometry optimisation for relaxations.

  • For vibrational contributions (IR intensities = dipole derivatives) see vibrational frequencies once the Hessian infrastructure lands.