Orbital tools

OrbTools

A collection of tools for working with orbitals

canonical_orthogonalization(sao::AbstractMatrix, redthr; verbose=false)

Construct the canonical orthogonalization transformation from the AO overlap matrix sao. Eigenvectors of sao whose eigenvalues are below redthr are considered linearly dependent (redundant) and are removed from the orbital basis. This makes the SCF robust for redundant basis sets (e.g. Cartesian basis sets with 6 instead of 5 d functions), where sao is (near-)singular.

Return (X, Xredundant) where X (nAO × nMO) fulfills X' sao X = I and is used to solve the SCF equations in an orthonormal basis, and Xredundant (nAO × nredundant) spans the removed near-null space. The latter is kept only to pad the MO coefficient matrices to square shape (nAO × nAO).

eigen_orth(fock::AbstractMatrix, X::AbstractMatrix, Xredundant::AbstractMatrix; large=1.0e6)

Solve the generalized eigenvalue problem fock C = sao C ϵ in the orthonormal basis defined by the canonical orthogonalization X (see canonical_orthogonalization), i.e. diagonalize X' fock X and back-transform. The redundant directions Xredundant are appended as virtual orbitals with energy large, so that they remain unoccupied and the returned coefficient matrix stays square (nAO × nAO).

Return (ϵ, cMO) (orbital energies and coefficients).

extend_classes_to_completed(classes_old, kept, n_new, nredundant) -> (classa, classb)

Build orbital classes for a completed new-basis orbital set (see project_onto_basis_complete) from classes_old = (classa, classb). kept[ispin] are the original orbital indices that survived as the leading columns of the completed set (they keep their original class); the intermediate complement orbitals are labelled "Virtual"; and the last nredundant (linearly-dependent) orbitals are labelled "Deleted", so they are excluded from the correlation treatment.

freeze_orbitals!(EC::ECInfo; MO="mo", redundant=true, verbose=true) -> (; occ_a, occ_b, virt_a, virt_b)

Apply the full frozen-core, redundant- and deleted-virtual selection for a correlated calculation, honoring the orbital classes stored in the dump (e.g. from @region) while letting the user override them:

• Occupied core: if wf.freeze_nocc ≥ 0, freeze that many lowest orbitals; otherwise if wf.core == :auto, freeze the orbitals tagged "Core" in the dump (falling back to the :large chemical core when the dump carries no class information); otherwise freeze the chemical core selected by wf.core.
• Virtuals: the linearly-dependent (redundant) orbitals and the (e.g. @region) deleted virtuals are both stored as class "Deleted". If wf.freeze_nvirt < 0 (auto), every "Deleted" orbital is frozen by its actual index (never by top index), so a region's active virtuals are never frozen by mistake; with no class info the recomputed AO redundancy (n_redundant_orbitals) is frozen as the highest virtuals. If wf.freeze_nvirt ≥ 0, exactly that many highest virtuals are frozen (the user decides), with a warning when this leaves redundant orbitals in the correlation space.

Pass redundant=false when the redundant orbitals were already excluded (e.g. a 3-index FCIDUMP). Returns the number of occupied/virtual orbitals removed, resolved by spin.

guess_orb(EC::ECInfo, guess::Symbol)

Calculate starting guess for MO coefficients. Type of initial guess for MO coefficients is given by guess.

See ScfOptions.guess for possible values.

guessposorb(EC::ECInfo, guess::Symbol)

Calculate starting guess for MO positron coefficients. Type of initial guess for MO coefficients is given by guess.

See ScfOptions.guess for possible values.

left_from_right_rotations(cMOr::SpinMatrix)

Calculate left biorthogonal rotation coefficients from right BO coefficients.

load_left_right_rotations(EC::ECInfo) -> (left::SpinMatrix, right::SpinMatrix)

Load (last) left and right orbital rotations from file WfOptions.dump.

If the type of the rotations does not contain the word biorthogonal, the same rotation is returned for left and right (can be checked with ===).

load_orbitals(EC::ECInfo; start::Bool=false)

Load (last) orbitals from file WfOptions.dump.

If start=true, load from wf.start instead. If the basis has changed, the orbitals will be projected onto the new basis. Returns ::SpinMatrix.

load_orbitals_for_correlation(EC::ECInfo; start::Bool=false) -> (cMO::SpinMatrix, classes)

Load orbitals for building a correlation FCIDUMP, honoring a geometry and/or basis change.

When the AO basis is genuinely unchanged (same geometry and basis), the stored orbitals are reused verbatim and classes = nothing is returned (freezing then uses the dump's own, already-matching classes). Otherwise — a geometry displacement and/or a different/larger/smaller basis, as in a dump=""+start restart — the stored orbitals are projected onto the current basis and Löwdin-orthonormalized (and completed for a size change) via project_onto_basis_complete, and the matching orbital classes for the resulting set are returned so freezing (frozen core / linearly-dependent orbitals) uses classes that describe the actual orbital set.

load_positron_orbitals(EC::ECInfo)

Load (last) positron orbitals from file WfOptions.dump.

Returns ::SpinMatrix.

load_rotations(EC::ECInfo)

Load (last) orbital rotations from file WfOptions.dump.

Returns ::SpinMatrix.

n_deleted_orbitals(EC::ECInfo; MO="mo")

Number of orbitals marked "Deleted" in the wavefunction dump, i.e. the linearly-dependent orbitals that the preceding (DF-)HF projected out. Post-HF methods freeze these out of the correlation treatment exactly like freeze_nvirt.

The stored count is the authoritative number actually used by the HF. It is cross-checked against the present redundancy of the AO basis (n_redundant_orbitals) and a warning is issued if they disagree (e.g. when scf.redthr was changed between the HF and the correlation step). If the dump contains no class/energy information, the recomputed redundancy is used as a fallback.

Redundant orbitals are identified as those both marked "Deleted" and parked at the sentinel energy REDUNDANT_ORBITAL_ENERGY; this distinguishes them from ordinary frozen virtuals, which are also stored as "Deleted" but keep their physical orbital energy.

n_deleted_orbitals_all(EC::ECInfo; MO="mo")

Number of orbitals stored as class "Deleted", i.e. ALL orbitals kept out of the correlation treatment: the linearly-dependent (redundant) ones AND virtuals removed on purpose (e.g. by @region), which keep their physical orbital energy. Use this where only the count of non-correlated orbitals matters (cost estimates); use n_deleted_orbitals when the linearly-dependent ones specifically are meant. Falls back to n_deleted_orbitals when no class information is stored (e.g. imported orbitals).

n_redundant_orbitals(EC::ECInfo)

Return the number of linearly-dependent (redundant) orbitals of the AO basis set, i.e. the number of eigenvalues of the AO overlap matrix below scf.redthr.

These orbitals are projected out by the (DF-)HF (see canonical_orthogonalization) and parked as the highest (unoccupied) orbitals. Post-HF methods freeze them out exactly like freeze_nvirt, so they do not enter the correlation treatment. Returns 0 if no molecular system is set up (e.g. for a plain FCIDump).

normalize_phase!(cMO)

Normalize the phase of the MO coefficients in cMO.

The phase is chosen such that the first largest coefficient is positive.

opao_relthr(EC::ECInfo)

Effective relative redundancy threshold for the OPAO (orthogonalized-PAO) construction in select_lowdin_orth: loc.opaofac × scf.redthr. Tying it to the AO basis redundancy threshold keeps the two consistent — only directions that are (near-)redundant by the same standard the basis uses are dropped — and basis-adaptive: lowering scf.redthr (a more complete basis) automatically loosens the PAO pruning.

orbital_classes_with_deleted(occ, norb, ndeleted)

Build a vector of TREXIO orbital classes of length norb for storing in the wavefunction dump: "Inactive" for the occupied orbitals occ, "Virtual" for the remaining orbitals, and "Deleted" for the last ndeleted orbitals (the linearly-dependent orbitals projected out by the (DF-)HF, see canonical_orthogonalization).

orbital_energies(EC::ECInfo, spincase::Symbol=:α)

Return orbital energies for a given spincase∈{:α,:β}.

project_onto_basis(cMO::SpinMatrix, old_basis::BasisSet, new_basis::BasisSet; check=false, redthr=1.0e-8)

Project the MO coefficients onto a new basis.

The projector $S_{new}^{-1} S_{new,old}$ uses the Moore-Penrose pseudo-inverse of the new-basis overlap, obtained from the canonical orthogonalization (see canonical_orthogonalization): with X' S_{new} X = I after dropping overlap eigenvalues below redthr, $S_{new}^{-1} = X X'$. This handles redundant / linearly-dependent bases (singular $S_{new}$) and reduces to the ordinary inverse when the basis is not redundant.

If check is true, the function will check whether the projection is needed and return the same array cMO if it is not (i.e., it can be checked with ===).

project_onto_basis_complete(cMO::SpinMatrix, old_basis::BasisSet, new_basis::BasisSet; redthr=1.0e-8)

Project cMO (given in old_basis) onto new_basis and complete it to a full orthonormal orbital set spanning the (non-redundant) new_basis AO space.

Unlike project_onto_basis, which keeps the original number of MOs, this returns a full set for the new basis, so a restarted correlation calculation actually uses the new basis:

• the projected orbitals are symmetric-Löwdin orthonormalized (kept as close as possible to the originals) and placed as the leading columns, preserving their occupied/virtual ordering;
• if new_basis is larger, the remaining space is filled with the orthogonal complement of the projected orbitals (built from the canonical orthogonalization of $S_{new}$) and appended as extra, higher virtual orbitals;
• if new_basis is smaller, the highest projected orbitals that no longer fit are dropped.

Consequently the returned set has n_ao(new_basis) columns (like a fresh (DF-)HF), and its leading min(n_old, n_new) columns correspond one-to-one to the (lowest) original orbitals — which lets a caller rebuild orbital classes as [classes_old truncated to the kept count ; "Virtual"…].

Redundant (linearly-dependent) new_basis sets are handled exactly like a fresh HF: the redundant directions (Xredundant from the canonical orthogonalization) are appended as the last columns so downstream freezing can drop them from the correlation treatment.

Returns (cMO_new::SpinMatrix, kept::Vector{Vector{Int}}, nredundant::Int), where kept[ispin] lists (in order) the original orbital indices that survived as the leading columns and nredundant is the number of trailing linearly-dependent columns — so a caller can rebuild orbital classes as [classes_old[kept] ; "Virtual"… ; "Deleted" × nredundant]. The survivor bookkeeping is exact even when linearly-dependent projected orbitals had to be dropped.

replace_core_from_dump!(EC::ECInfo, cMO::SpinMatrix, ncore::Int) -> SpinMatrix

Replace the ncore frozen-core orbitals (the leading columns) of cMO with the corresponding core orbitals read from the WfOptions.dump file (e.g., a fresh HF at the current geometry), and orthonormalize the remaining (correlating) orbitals against the new core (and each other) by projection — modified Gram–Schmidt in the current AO overlap metric, seeded with the fixed core. Used for the dump4core_only geometry-change restart: the correlating orbitals in cMO come from start, but their frozen core is stale (stuck at the previous geometry); this swaps in the correct new-geometry core. Returns the modified cMO.

rotate_orbs(EC::ECInfo, orb1, orb2, angle=90; spin::Symbol=:α)

Rotate orbitals orb1 and orb2 from WfOptions.dump by angle degrees. For unrestricted orbitals, spin can be :α or :β.

rotate_orbs!(cMO::Matrix, orb1, orb2, angle=90)

Rotate orbitals orb1 and orb2 from cMO by angle degrees.

cMO is a matrix of MO coefficients.

select_lowdin_orth(S::AbstractMatrix; relthr=3e-8, verbose=false)

Locality-preserving orthogonalizer of a Hermitian metric S (typically a projected atomic-orbital (PAO) overlap S_PAO = C_PAO' S C_PAO), used to build orthogonal PAOs (OPAOs) for the virtual space without redundancies and without delocalizing them.

A single Hermitian eigendecomposition drives both the rank and the column selection:

1. Rank r from the eigenvalues of S with a relative threshold: directions with eigenvalue > relthr * λmax are kept. With relthr tied to the basis redundancy threshold (relthr = loc.opaofac * scf.redthr, a conservative ~1e-8) only directions that are (near-)redundant by the same standard the basis uses are dropped; small-but-real directions (e.g. the virtual residual of a frozen core AO, which sits at ~1e-7 λmax, well above redthr) are kept as genuine degrees of freedom.
2. Selection of exactly r actual (atom-centered) PAOs spanning the retained subspace, via a rank-revealing column-pivoted QR of the r kept eigenvectors Vₖ (the r most linearly-independent rows of Vₖ). Working on the clean, orthonormal kept subspace — instead of a pivoted Cholesky of the rank-deficient metric S — keeps the selection well-conditioned; pinning the count to r from step 1 avoids relying on a fragile pivot threshold.
3. Symmetric Löwdin orthogonalization S_sub^{-1/2} of the selected (clean, well-conditioned) r×r block, which keeps each OPAO as close as possible to its parent PAO (maximally local; cf. canonical_orthogonalization, which instead rotates into overlap eigenvectors and delocalizes), followed by one symmetric Löwdin refinement on the (≈ identity) residual metric M' S M to restore machine-precision orthonormality when a low-presence direction is amplified.

Returns M (size(S,1) × r) with M' S M = I, a drop-in for sqrtinvchol(S; max_rank=...) in C_OPAO = C_PAO * M. relthr is the OPAO redundancy threshold, opao_relthr(EC) = loc.opaofac * scf.redthr.

show_orbitals(EC::ECInfo, range=nothing)

Print the MO coefficients in WfOptions.dump with respect to the atomic orbitals.

range is a range of molecular orbitals to be printed.

show_orbitals(EC::ECInfo, cMO::Matrix, basis::BasisSet, range=1:size(cMO,2)

Print the MO coefficients in cMO with respect to the atomic orbitals in basis.

range is a range of molecular orbitals to be printed.

try_load_starting_orbitals(EC::ECInfo) -> (SpinMatrix, Bool)

Try to load starting orbitals from wf.start file.

If wf.start is set, load and project orbitals from that file.

Returns (cMO, loaded) where loaded indicates if orbitals were successfully loaded. If no orbitals are available, returns (SpinMatrix{Float64}(), false).

_independent_columns_in_metric(C::AbstractMatrix, S::AbstractMatrix, thr) -> Vector{Int}

Return the indices of the columns of C that are linearly independent in the metric S, processing columns left-to-right (so earlier columns are preferred/protected). A column is kept if its S-norm after orthogonalization against the already-kept columns exceeds thr.

_remove_orbitals_spin!(EC::ECInfo, occ_a, occ_b, virt_a, virt_b)

Remove the given (spin-resolved) occupied and virtual orbital indices from the active subspaces and rebuild the derived spin spaces (d/s/S), mirroring setup_space!.

_try_orbital_classes_energies(EC::ECInfo; MO="mo")

Best-effort fetch of (classa, classb, ea, eb) from the dump, returning empty vectors when no class/energy information is available (e.g. an FCIDUMP-only run).

guess_hcore(EC::ECInfo)

Guess MO coefficients from core Hamiltonian.

guess_pos_hcore(EC::ECInfo)

Guess MO coefficients for positron from core Hamiltonian.

guess_positron(EC::ECInfo)

Initialize positron MO coefficients as zeroes.

guess_sad(EC::ECInfo)

Guess MO coefficients from atomic densities.