"""Hamiltonian construction from AO integrals.
Provides a clean interface for building :class:`Hamiltonian` objects
from vibe-qc's native C++ integral bindings, with optional orbital
transformation.
"""
from __future__ import annotations
from typing import Optional
import numpy as np
from .._vibeqc_core import (
BasisSet,
Molecule,
compute_eri,
compute_kinetic,
compute_nuclear,
compute_overlap,
)
from ._common import Hamiltonian, _chemist_to_physicist
def build_hamiltonian_ao(
molecule: Molecule,
basis: BasisSet,
*,
description: str = "",
) -> Hamiltonian:
"""Build a :class:`Hamiltonian` in the **AO basis**.
The returned ``h1e`` and ``h2e`` are in the original (non-orthogonal)
atomic-orbital basis. Most solvers require an orthonormal basis;
transform with :func:`transform_hamiltonian_to_mo` after orthogonalising.
Parameters
----------
molecule : Molecule
Molecular geometry, charge, multiplicity.
basis : BasisSet
AO basis set.
description : str
Label for logging.
Returns
-------
Hamiltonian
With ``h1e``, ``h2e``, ``nuclear_repulsion``, ``nelec``, ``ms2``,
and ``norb`` populated.
"""
S = np.asarray(compute_overlap(basis), order="C")
T = np.asarray(compute_kinetic(basis), order="C")
V = np.asarray(compute_nuclear(basis, molecule), order="C")
h1e_ao = T + V # (norb, norb)
# Chemist's notation (muν|ls) from libint -> physicist's g_{muνls}
eri_chem_ao = np.asarray(compute_eri(basis), order="C")
h2e_ao = _chemist_to_physicist(eri_chem_ao) # g_{muνls} = (mul|νs)
norb = h1e_ao.shape[0]
nelec = molecule.n_electrons()
multiplicity = molecule.multiplicity
ms2 = multiplicity - 1 # 2*S_z
return Hamiltonian(
h1e=h1e_ao,
h2e=h2e_ao,
nuclear_repulsion=molecule.nuclear_repulsion(),
norb=norb,
nelec=nelec,
ms2=ms2,
description=description or f"{molecule.n_electrons()}e, {norb} orbitals",
)
[docs]
def build_hamiltonian_mo(
molecule: Molecule,
basis: BasisSet,
mo_coeffs: np.ndarray,
*,
description: str = "",
) -> Hamiltonian:
"""Build a :class:`Hamiltonian` in the **MO basis**.
Transforms the one- and two-electron integrals from the AO basis
into the molecular-orbital basis given by ``mo_coeffs``.
Parameters
----------
molecule : Molecule
basis : BasisSet
mo_coeffs : (n_ao, n_mo) ndarray
AO -> MO coefficient matrix (columns are MOs).
description : str
Returns
-------
Hamiltonian
With integrals in the MO basis. ``norb`` = n_mo.
"""
ham_ao = build_hamiltonian_ao(molecule, basis, description=description)
return transform_hamiltonian(ham_ao, mo_coeffs)
[docs]
def get_hf_orbital_provider(
molecule: Molecule,
basis: BasisSet,
*,
method: str = "rhf",
) -> np.ndarray:
"""Run vibe-qc's RHF and return the MO coefficient matrix.
This is a **convenience** function -- the solver layer does not
require HF orbitals. Any source of orthonormal orbitals works.
Parameters
----------
molecule : Molecule
basis : BasisSet
method : str
``"rhf"`` (default) or ``"uhf"``.
Returns
-------
mo_coeffs : (n_ao, n_ao) ndarray
Columns are canonical HF molecular orbitals.
"""
from .._vibeqc_core import RHFOptions, UHFOptions, run_rhf, run_uhf
if method == "rhf":
opts = RHFOptions()
opts.max_iter = 200
opts.conv_tol_energy = 1e-12
opts.conv_tol_grad = 1e-10
result = run_rhf(molecule, basis, opts)
return np.asarray(result.mo_coeffs, order="C")
elif method == "uhf":
opts = UHFOptions()
opts.max_iter = 200
opts.conv_tol_energy = 1e-12
opts.conv_tol_grad = 1e-10
result = run_uhf(molecule, basis, opts)
return np.asarray(result.mo_coeffs_alpha, order="C")
elif method == "uno":
# UHF natural orbitals (UNO-CAS) -- P. Pulay & T. P. Hamilton,
# J. Chem. Phys. 88, 4926 (1988), doi:10.1063/1.454704: diagonalize
# the total UHF density in the symmetrically-orthogonalized AO
# basis,
# n, U = eigh(S^{1/2} (Pa + Pb) S^{1/2}), C = S^{-1/2} U,
# ordered by descending occupation n. Occupations cluster near 2
# (doubly occupied), 1 (open shell) and 0 (virtual), so a single
# spin-restricted orbital set carries the open-shell character and
# the descending order matches the CAS family's lowest-core /
# active-window convention. This is the recommended starting
# reference for open-shell CASCI / CASSCF / MR-PT2 (the plain
# "uhf" choice uses the a orbitals only, leaving the b space -- and
# hence the doubly-occupied core -- only approximately represented).
from .._vibeqc_core import compute_overlap
opts = UHFOptions()
opts.max_iter = 200
opts.conv_tol_energy = 1e-12
opts.conv_tol_grad = 1e-10
result = run_uhf(molecule, basis, opts)
dm_total = np.asarray(result.density_alpha) + np.asarray(
result.density_beta
)
s_ao = np.asarray(compute_overlap(basis))
s_val, s_vec = np.linalg.eigh(s_ao)
if float(np.min(s_val)) < 1e-10:
raise ValueError(
"AO overlap is near-singular (min eigenvalue "
f"{np.min(s_val):.2e}); UNO construction needs a "
"non-degenerate basis."
)
s_half = (s_vec * np.sqrt(s_val)) @ s_vec.T
s_inv_half = (s_vec / np.sqrt(s_val)) @ s_vec.T
occ, u_no = np.linalg.eigh(s_half @ dm_total @ s_half)
order = np.argsort(occ)[::-1] # descending occupation
return np.ascontiguousarray(s_inv_half @ u_no[:, order])
elif method == "rohf":
# Restricted open-shell HF reference (Roothaan, Rev. Mod. Phys. 32,
# 179 (1960)). ROHF gives ONE spin-restricted orbital set already
# ordered closed (doubly occ) -> open (singly occ) -> virtual, so
# the doubly-occupied core and the active window are well-defined
# AND the determinant is spin-pure (<S^2> = S(S+1) exactly) -- a
# cleaner CAS starting reference than the spin-contaminated UHF
# alpha set, and an alternative to UNO that needs no natural-orbital
# diagonalisation.
from ..rohf import ROHFOptions, run_rohf
opts = ROHFOptions()
opts.max_iter = 200
opts.conv_tol_energy = 1e-12
opts.conv_tol_grad = 1e-9
result = run_rohf(molecule, basis, opts)
return np.asarray(result.mo_coeffs, order="C")
else:
raise ValueError(f"Unknown HF method: {method!r}")
def canonical_orthogonalize(
S: np.ndarray,
threshold: float = 1e-8,
) -> np.ndarray:
"""Return the canonical (Löwdin-symmetric) orthogonalisation matrix X.
X = U s^{-1/2} U^T where S = U s U^T is the overlap eigen-decomposition.
This produces orthonormal orbitals closest to the original AO basis.
Parameters
----------
S : (n, n) ndarray
Overlap matrix.
threshold : float
Eigenvalue threshold for linear-dependence removal.
Returns
-------
X : (n, n_active) ndarray
Orthogonalisation matrix: C_orth = X^T . C_ao.
"""
evals, evecs = np.linalg.eigh(S)
mask = evals > threshold
s_inv_sqrt = np.diag(1.0 / np.sqrt(evals[mask]))
return evecs[:, mask] @ s_inv_sqrt