Orbital Localization

Orbital localization methods.

Implements Intrinsic Bond Orbitals (IBO) (Reference: G. Knizia, JCTC 2013, 9, 4834-4843), Pipek-Mezey localization with Mulliken charges, and Foster-Boys localization for occupied orbital localization and orthogonal Projected Atomic Orbitals (OPAO) for virtual orbital localization.

compute_localization_rotations(EC::ECInfo; exponent=4)

Compute orbital localization rotation matrices for occupied (IBO) and virtual (OPAO) orbitals.

Returns (R_occ, R_virt) where:

• R_occ: nocc × nocc rotation matrix (canonical → IBO-localized occupied)
• R_virt: nvirt × nvirt rotation matrix (canonical → OPAO-localized virtual)

The rotation matrices are orthogonal/unitary and can be used to transform amplitude matrices before SVD decomposition.

localize_orbitals(EC::ECInfo)

Localize the current orbitals using IBO, PM, or Boys (occupied) and optionally OPAO (virtual).

Reads orbitals from the wavefunction dump, applies localization to the occupied orbitals and optionally OPAO rotation to the virtual orbitals, then stores the localized orbitals back to the dump file.

Options are read from EC.options.loc:

• method::String: "ibo" (default), "pm" (Pipek-Mezey with Mulliken charges), or "boys" (Foster-Boys, maximizes sum of squared orbital dipole moments).
• virtual::Bool: if true (default), also localize virtual orbitals via OPAO.
• exponent::Int: localization exponent (0=auto: 4 for IBO, 2 for PM; Boys ignores exponent).

Examples

@dfhf
@localize              # localize occupied (IBO) + virtual (OPAO)
@set loc method="pm"
@localize              # Pipek-Mezey localization (occupied) + OPAO (virtual)
@set loc method="boys"
@localize              # Foster-Boys localization (occupied) + OPAO (virtual)
@set loc virtual=false
@localize              # localize only occupied

_ao_compactness_sweep!(M, R, orbs; maxiter=50, tol=1e-10)

Maximize AO compactness $\sum_i \sum_\mu |M_{\mu i}|^4$ by 2×2 Jacobi rotations among the orbital indices in orbs.

Both M (nbas × norb representation matrix) and R (norb × norb rotation accumulator) are modified in-place. This can be used as a starting guess (pre-localization) or as a post-convergence refinement to break degeneracies.

_canonical_orbital_order(charges, M, nocc, natom)

Determine a canonical permutation for localized orbitals.

Sorts orbitals by: (1) dominant atom index (ascending), (2) charge on dominant atom (descending), (3) full charge vector (descending), and (4) an abs² coefficient fingerprint in the representation matrix M as a deterministic tiebreaker for orbitals localized on the same atom. The coefficient fingerprint is compared in the native basis order of M, which is stronger than using only sorted magnitudes and reliably distinguishes symmetry-related same-atom orbitals.

_compute_proto_iaos(EC::ECInfo, cMO_occ::AbstractMatrix)

Construct the non-orthogonal proto-IAOs using the revised Senjean/Knizia formulation based on the depolarized occupied orbitals from Eq. (S20) and the efficient matrix factorization from Eq. (S44) of the supplementary information to Senjean et al., JCTC 2021, 17, 1337-1354.

Returns (proto_iao, iao_atoms, natom, S) where S is the AO overlap matrix for the main basis.

_deterministic_group_rotation!(M, R, group)

Apply a deterministic rotation within a degenerate orbital group.

Replaces the iterative AO compactness sweep for degeneracy refinement with a non-iterative eigendecomposition, ensuring platform-independent results.

Computes a basis-index-weighted overlap matrix within the group and diagonalizes it. The eigenvectors (sorted by eigenvalue) define a unique rotation that breaks the degeneracy deterministically. This avoids the platform-dependent convergence behavior of Jacobi sweeps within degenerate subspaces.

_find_degenerate_groups(charges::AbstractMatrix, nocc::Int; charge_tol=1e-6)

Identify groups of orbitals with identical atom-resolved charges (degenerate under localization). Returns a vector of groups, where each group is a vector of orbital indices.

_fix_sign_convention!(R::AbstractMatrix)

Fix the sign ambiguity of column vectors in a rotation matrix.

For each column, ensures the element with the largest absolute value is positive. This makes the rotation matrix deterministic across platforms and Julia versions, where Jacobi sweeps or eigendecompositions may converge to equivalent solutions that differ by column signs.

_jacobi_localization!(R, charges, natom, nocc,
                      compute_qAij, rotate_workspace!, update_charges!;
                      exponent, maxiter, tol, method_name)

Shared Jacobi sweep engine for orbital localization.

Maximizes $\sum_i \sum_A (q_A^i)^p$ over 2×2 rotations in the occupied space. The method-specific operations (charge computation, workspace rotation, charge updates) are provided via typed callback arguments for both type stability and code reuse.

Arguments

• R: nocc × nocc rotation accumulator (modified in-place)
• charges: natom × nocc partial charge matrix (modified in-place)
• natom, nocc: dimensions
• compute_qAij(i, j, A): returns off-diagonal charge $q_A^{ij}$ for orbitals i,j on atom A
• rotate_workspace!(i, j, c, s): applies 2×2 rotation to method-specific workspace columns
• update_charges!(charges, i, j): recomputes charges for columns i and j after rotation

_metric_orthogonalize(C::AbstractMatrix, S::AbstractMatrix; tol=1e-12)

Löwdin-orthogonalize the columns of C with respect to the metric S.

compute_ao_atoms(EC::ECInfo)

Compute the atom assignment for each AO basis function.

Returns (ao_atoms, natom) where:

• ao_atoms: vector mapping each AO to its atom index (1-based, 0 for ghost atoms)
• natom: number of unique non-ghost atoms

compute_iaos(EC::ECInfo, cMO_occ::AbstractMatrix)

Compute Intrinsic Atomic Orbitals (IAOs) following Knizia, JCTC 2013, 9, 4834.

The proto-IAOs are built with the revised depolarized-orbital formulation recommended in the supplementary information to Senjean et al., JCTC 2021, 17, 1337-1354, and are then Löwdin-orthogonalized under the AO overlap metric.

The IAOs are constructed from a minimal basis and the occupied MOs. The minimal basis is determined by EC.options.loc.minao. If empty (default), the "minao" type from the basis Dict is used; otherwise "minao" as fallback. All nmin IAOs are returned (not just nocc), ensuring each atom contributes IAOs that can capture bonding character in the IBO charge analysis. Ghost atom minimal basis functions are excluded.

Returns (C_iao, iao_atoms, natom) where:

• C_iao: IAO coefficient matrix in AO basis (nAO × nmin), S-orthonormal
• iao_atoms: vector mapping each IAO to its atom index (1-based, contiguous)
• natom: number of unique (non-ghost) atoms represented in the IAOs

compute_opao_rotation(cMO_virt::AbstractMatrix, S::AbstractMatrix; relthr=3e-8)

Compute the rotation matrix for virtual orbitals based on orthogonalized PAOs.

PAOs (Projected Atomic Orbitals) are constructed by projecting AO basis functions onto the virtual space: $C_{\text{PAO}} = C_{\text{virt}} C_{\text{virt}}^T S$. The PAO overlap matrix is orthogonalized via select_lowdin_orth (relative- threshold rank detection + atom-centered pivot selection + symmetric Löwdin), which removes redundant PAOs while keeping the OPAOs local. relthr is the OPAO redundancy threshold, usually loc.opaofac * scf.redthr.

Returns R_virt (nvirt × nvirt) such that C_virt_loc = C_virt * R_virt.

localize_boys(cMO_occ::AbstractMatrix, S::AbstractMatrix,
              Dx::AbstractMatrix, Dy::AbstractMatrix, Dz::AbstractMatrix;
              maxiter=500, tol=1e-10)

Localize occupied orbitals using the Foster-Boys criterion.

Maximizes $\sum_i |\langle i | \mathbf{r} | i \rangle|^2 = \sum_i \sum_{\alpha=x,y,z} (d_\alpha^i)^2$ where $d_\alpha^i = \langle i | r_\alpha | i \rangle$ are the orbital dipole moments, using 2×2 Jacobi rotations.

Arguments

• cMO_occ: occupied MO coefficients (nAO × nocc)
• S: AO overlap matrix (nAO × nAO)
• Dx, Dy, Dz: AO dipole integral matrices $\langle \mu | r_\alpha | \nu \rangle$
• maxiter: maximum number of Jacobi sweeps
• tol: convergence threshold for rotation angles

Returns the rotation matrix R_occ (nocc × nocc) that transforms canonical to localized occupied MOs: C_occ_loc = C_occ * R_occ.

localize_ibo(cMO_occ::AbstractMatrix, S::AbstractMatrix, C_iao::AbstractMatrix,
             iao_atoms::Vector{Int}, natom::Int; exponent=4, maxiter=500, tol=1e-10)

Localize occupied orbitals using the IBO criterion.

Maximizes $\sum_i \sum_A (q_A^i)^p$ where $q_A^i = \sum_{\mu \in A} |\langle iao_\mu | \phi_i \rangle|^2$ using 2×2 Jacobi rotations.

Returns (R_occ, charges) where:

• R_occ is the rotation matrix (nocc × nocc) that transforms canonical to localized occupied MOs, C_occ_loc = C_occ * R_occ
• charges is the final atom-resolved IBO charge matrix (natom × nocc)

localize_pm(cMO_occ::AbstractMatrix, S::AbstractMatrix,
            ao_atoms::Vector{Int}, natom::Int; exponent=2, maxiter=500, tol=1e-10)

Localize occupied orbitals using the Pipek-Mezey criterion with Mulliken charges.

Maximizes $\sum_i \sum_A (q_A^i)^p$ where $q_A^i = \sum_{\mu \in A} \sum_\nu C_{\mu i} S_{\mu\nu} C_{\nu i}$ (Mulliken partial charges) using 2×2 Jacobi rotations.

Arguments

• cMO_occ: occupied MO coefficients (nAO × nocc)
• S: AO overlap matrix (nAO × nAO)
• ao_atoms: vector mapping each AO to its atom index (1-based, 0 for ghost)
• natom: number of unique (non-ghost) atoms
• exponent: power of charges in functional (default 2, the standard PM functional)
• maxiter: maximum number of Jacobi sweeps
• tol: convergence threshold for rotation angles

Returns the rotation matrix R_occ (nocc × nocc) that transforms canonical to localized occupied MOs: C_occ_loc = C_occ * R_occ.