"""Crystal-orbital diagrams from symmetry + extended-Hueckel overlap
(crystod --diagram).
The crystalline counterpart of the molecular-orbital diagram of
``crystod-mol --diagram --ao-left ... --ao-right ...``: the two fragment
sublattices given by ``--co-left``/``--co-right`` (e.g. the SrTi cation
framework and the O3 anion framework of SrTiO3) are treated with their
full-electron basis -- every core and valence shell of every atom
(WIEN2k-style; core shells with Slater-rule exponents and the archived
neutral-atom PySCF Hartree-Fock levels of reference/atomic_level_*,
collected into crystod/atomic_levels.py; shells frozen into the def2
effective core potential beyond Kr are omitted) -- and each fragment
feels the removed
sublattice as a point-charge lattice with the formal oxidation states
(the Madelung ligand field; see crystod.point_charge_field), so their
Bloch states are the complete electronic states before chemical bond
formation. At every high-symmetry k point,
1. the Bloch orbitals of each fragment sublattice are symmetry-adapted
(the crystal-orbital irreps of ``crystod --atomic-orbital``, i.e. the
site-symmetry induced representations),
2. all inter- and intra-sublattice overlap integrals are evaluated exactly
as Bloch lattice sums of single/double-zeta STO overlaps,
3. the generalized eigenvalue problem with the Wolfsberg-Helmholz
Hamiltonian H_ij = K S_ij (H_ii + H_jj)/2 is solved -- fragment
orbitals sharing an irrep of the little group mix into bonding and
antibonding crystal orbitals, orbitals whose irrep finds no partner
remain nonbonding -- exactly the construction rule of the crystal
orbital diagram (COD),
4. the result is written as an interactive HTML diagram with one energy
diagram per k point (fragment | crystal | fragment columns; the energy
window opens on -20 .. 10 eV, "Show all energy levels" reveals the
deep shells).
Every level carries a hover wave-function sketch: the Re[psi] amplitudes
of all its atomic-orbital components, rendered on the k-commensurate
supercell (same-l shells accumulate with their radial weights, so the
drawn lobe signs are those of the real wave function in the bonding
region).
Reference: Y. Mochizuki, M. Nishibori and T. Fukushima, "Crystal Orbital
Diagram of Perovskites: A Revisit from Symmetry-Adapted Linear
Combination" (in preparation).
"""
from __future__ import annotations
import argparse
import re
from dataclasses import dataclass, field
import numpy as np
from .atomic_levels import ATOMIC_LEVELS
from .mo_diagram import (
ANGSTROM_TO_BOHR,
CORE_SHELLS,
EHT_PARAMETERS,
WOLFSBERG_HELMHOLZ_K,
AtomicOrbital,
make_aligned_cache,
pair_overlap,
render_diagram_page,
svg_sub_digits,
)
from .point_charge_field import (
COULOMB_EV_ANGSTROM,
ewald_site_potential,
point_charge_block,
radial_overlap,
slater_zeta,
)
from .runtime_compat import get_character, get_chemical_symbols, get_scaled_positions
from .spglib_compat import ensure_spglib_compat
ensure_spglib_compat()
from phonopy.structure.cells import get_primitive_matrix_by_centring
from spgrep.core import get_spacegroup_irreps_from_primitive_symmetry
from .irreptables_compat import load_irreptables
from .operations import wigner_D_real, snap_qpoint
from .visualize_basis import SymmetryAdaptedOrbitalBasis
IrrepTable, _Irrep = load_irreptables()
_DEGENERACY_TOL = 1e-5
# default view of the interactive energy window: +-8 eV around the HOMO/LUMO
# midpoint (the VBM/CBM region one usually inspects first); the "Show all
# energy levels" button reveals everything outside it. The fixed window below
# is only the fallback when no HOMO/LUMO pair exists.
_VIEW_HALF_WINDOW = 8.0
_VIEW_E_MIN = -20.0
_VIEW_E_MAX = 10.0
# Bloch combinations whose overlap eigenvalue falls below this floor are
# excluded from the variational solve. Diffuse cation valence shells
# (e.g. Sr 5s/5p, Ti 4s/4p) overlap so strongly in a dense sublattice that
# some Bloch combinations become nearly expressible by the rest of the
# basis; for those the extended-Hueckel H is no longer consistent and the
# energies diverge as (1-K) H_ii / eigenvalue (the well-known EHT overlap
# catastrophe), polluting even the occupied manifold (measured on
# ScF3/SrTiO3: catastrophic modes all have eigenvalue <= 0.19, physical
# ones >= 0.26). They are still genuine Bloch states, though: in a dense
# sublattice a SINGLE shell's own Bloch sum can fall below the floor
# (rocksalt AlN, cation fcc neighbours at 2.86 A: the one Al 3s X1+
# combination has eigenvalue 0.15), and deleting it removes a whole
# physical level and breaks the aufbau electron counts. Such modes are
# therefore kept as separate levels with a first-order Loewdin energy
# estimate (see _generalized_eigh); only truly linearly dependent
# combinations (eigenvalue below _DEPENDENT_TOL, no physical content)
# are removed outright.
_OVERLAP_FLOOR = 0.2
_DEPENDENT_TOL = 1e-6
# tooltip note attached to first-order-estimated levels (the terminal
# marks their energies with ~; the HTML keeps normal solid lines --
# user preference -- and carries this note in the tooltip)
_ESTIMATED_NOTE = (
"~ near-dependent Bloch combination (overlap eigenvalue {eps:.2f} < "
"floor {floor}): the energy is a first-order Loewdin estimate; the "
"variational extended-Hueckel value diverges (overlap catastrophe)")
@dataclass
class SublatticeSpec:
"""One (element, shell) block of a fragment sublattice.
Attributes:
element: Chemical symbol.
letter: Shell letter ``s``, ``p``, ``d`` or ``f``.
shell: Shell name such as ``"3d"``.
n: Principal quantum number.
l: Azimuthal quantum number.
zeta: STO exponent, a scalar or ``[(zeta, coefficient), ...]`` for a
double-zeta shell.
h_ii: On-site energy in eV.
sites: Indices of the atoms carrying the shell.
column: ``"left"`` or ``"right"``.
offset: First AO index of the block in the full basis.
"""
element: str
letter: str # s / p / d
shell: str # e.g. "3d"
n: int
l: int
zeta: object # scalar or [(zeta, coeff), ...]
h_ii: float
sites: list[int]
column: str # "left" / "right"
offset: int = 0 # first AO index of this spec
@property
def n_ao(self) -> int:
return len(self.sites) * (2 * self.l + 1)
@dataclass
class DiagramLevel:
"""One energy level of a diagram column, as returned by ``solve_at``.
Attributes:
level_id: Unique id within the k point, e.g. ``"mo3"``.
column: ``"left"``, ``"mo"`` (the crystal) or ``"right"``.
energy: Energy in eV.
degeneracy: Number of partners (the irrep dimension, or a multiple).
irrep: Bare irrep name, e.g. ``"GM4-"``.
label: Display label, e.g. ``"Sc 3d GM5+"`` or ``"GM5+ #1"``.
electrons: Electrons in the level after the aufbau filling.
vectors: AO coefficients, shape ``(n_ao, degeneracy)``,
S-orthonormal.
composition: ``[(level_id, weight)]`` of the fragment levels a
crystal level is made of (Loewdin-orthogonalized weights).
detail: Tooltip text: populations, bond character, notes.
estimated: The energy is a first-order Loewdin estimate of a
near-dependent Bloch combination, not the (divergent)
variational extended-Hueckel value.
"""
level_id: str
column: str # "left" / "mo" / "right"
energy: float
degeneracy: int
irrep: str # bare irrep name, e.g. GM4-
label: str # display label
electrons: int = 0
vectors: np.ndarray | None = None # (n_ao, degeneracy), S-orthonormal
composition: list = field(default_factory=list) # [(level_id, weight)]
detail: str = ""
estimated: bool = False # near-dependent Bloch combination: energy is
# a first-order Loewdin estimate, not the
# (divergent) variational EHT value
def parse_fragment_formula(tokens: list[str], flag: str) -> list[tuple[str, int | None]]:
"""Parse formula tokens (SrTi, O3, ...) into (element, count?) pairs."""
pairs = []
for token in tokens:
if not re.fullmatch(r"(?:[A-Z][a-z]?\d*)+", token):
raise SystemExit(
f"ERROR: invalid {flag} formula '{token}' (expected element "
"symbols with optional counts, e.g. SrTi or O3)."
)
for element, count in re.findall(r"([A-Z][a-z]?)(\d*)", token):
pairs.append((element, int(count) if count else None))
return pairs
def parse_oxidation_tokens(tokens: list[str]) -> dict[str, float]:
"""Parse El=Q tokens (Sr=+2, O=-2, Ti=4) into {element: charge}."""
oxidation = {}
for token in tokens:
match = re.fullmatch(r"([A-Z][a-z]?)=([+-]?\d+(?:\.\d+)?)", token)
if not match:
raise SystemExit(
f"ERROR: invalid --oxidation token '{token}' "
"(expected e.g. Sr=+2 Ti=+4 O=-2)."
)
oxidation[match.group(1)] = float(match.group(2))
return oxidation
def parse_sketch_tokens(tokens: list[str]) -> list[tuple[str, str]]:
"""Parse El-shell tokens (Ti-3d, Ti-d, O_2p) into (element, shell)."""
wanted = []
for token in tokens:
parts = re.split(r"[-_]", token)
if len(parts) != 2 or not re.fullmatch(r"\d?[spd]", parts[1]):
raise SystemExit(
f"ERROR: invalid --atomic-orbital token '{token}' "
"(expected e.g. Ti-3d, Ti-d or O_2p)."
)
wanted.append((parts[0], parts[1]))
return wanted
# |overlap population| below this is displayed as nonbonding: it is exactly
# 0 for symmetry-nonbonding states (irrep without a partner on the other
# sublattice/fragment), while genuine bonding/antibonding states come out
# at |P| = 0.02-6 (antibonding |P| is systematically the larger, the usual
# non-orthogonal COOP asymmetry)
BOND_CHARACTER_TOL = 0.01
# a fragment shell whose own occupied levels all lie this far below the
# fragment's HOMO is semicore: its orthogonality tails carry an
# antibonding-signed overlap population that masks the valence character
SEMICORE_DEPTH_EV = 10.0
[docs]
def assign_bond_characters(levels, overlap, rows_left, rows_right,
spec_ranges, sqrt_overlap=None,
hamiltonian=None) -> None:
"""COOP bonding character of every crystal (or molecular) orbital level.
``P = 2 Re[c_L+ S_LR c_R] / degeneracy`` is the electron weight
accumulated between the two fragments: ``P > 0`` in-phase (bonding,
drawn blue), ``P < 0`` out-of-phase with an internuclear node
(antibonding, red), ``P ~ 0`` nonbonding (black; exactly 0 when the
irrep has no partner on the other fragment). This is the
classification that colors the level connectors of ``crystod --diagram``
and ``crystod-mol --diagram``; the diagram engines call it from their
``solve_at`` after the aufbau filling. An ``(F - E S)`` energy
partition was rejected: Mulliken-like cross terms of the diffuse shells
give nonsense signs in a non-orthogonal basis.
Semicore handling. A filled semicore shell contributes to ``P`` in two
distinct ways: resonant filled-filled pairing (Sc 3p with F 2s of ScF3,
6 eV apart -- the He2-like closed-shell repulsion whose occupied upper
partner is genuinely antibonding) and far off-resonant orthogonality
tails (the same Sc 3p inside the F 2p band 23 eV above, or Sc 3s
against everything), which are not bonding physics and would flip the
sign of an otherwise donation-bonding state. With ``hamiltonian``
given (the crystal Fock/EHT operator), semicore shells are flagged
against the crystal valence-band maximum -- occupied fragment levels
whose expectation ``<phi|H|phi>`` (reference-consistent with the
crystal energies) tops out ``SEMICORE_DEPTH_EV`` below the VBM -- and a
flagged shell is excluded from a level's ``P`` only when the level is
more than ``SEMICORE_DEPTH_EV`` away from that shell's band top (the
semicore band and its resonant partners keep it). Without
``hamiltonian`` the legacy rule applies: flag against each fragment
column's own HOMO and keep the shell where it holds at least 40% of the
level's weight -- fine for molecules, but blind to a semicore that is
the fragment HOMO (Sc 3p of Sc3+).
Args:
levels: ``{"left": [...], "mo": [...], "right": [...]}`` level
records (``DiagramLevel`` or the molecular equivalent) with
``vectors``, ``energy``, ``degeneracy``, ``electrons`` and
``label`` (``"El shell irrep"``) filled in.
overlap: The overlap matrix ``S`` of the full basis at this k point.
rows_left: AO indices of the left fragment's basis functions.
rows_right: AO indices of the right fragment's basis functions.
spec_ranges: ``{(element, shell): AO indices}`` of every shell block.
sqrt_overlap: ``S^(1/2)`` for Loewdin weights in the legacy semicore
rule; ``None`` uses Mulliken gross populations there.
hamiltonian: The crystal one-electron operator (Fock or EHT ``H``)
at this k point; enables the VBM-referenced semicore rule.
Returns:
``None``. Every ``levels["mo"]`` record receives ``bond_character``
(``"bonding"``, ``"antibonding"`` or ``"nonbonding"``, threshold
``BOND_CHARACTER_TOL``) and ``overlap_population`` (``P``), and a
line stating them is appended to its ``detail``.
"""
semicore_tops: dict[tuple[str, str], float] = {}
if hamiltonian is not None:
vbm = max((lv.energy for lv in levels["mo"] if lv.electrons > 0),
default=0.0)
tops: dict[tuple[str, str], float] = {}
for column in ("left", "right"):
for lv in levels[column]:
if lv.electrons <= 0:
continue
parts = lv.label.split()
if len(parts) < 3 or (parts[0], parts[1]) not in spec_ranges:
continue
if getattr(lv, "estimated", False):
# the variational <H> of a near-dependent combination IS
# the divergent overlap-catastrophe value; its displayed
# first-order Loewdin estimate is the meaningful energy
# (EHT engine only, single Hamiltonian -- same reference)
expectation = lv.energy
else:
expectation = float(np.trace(
lv.vectors.conj().T @ hamiltonian @ lv.vectors
).real) / lv.degeneracy
key = (parts[0], parts[1])
tops[key] = max(tops.get(key, -1e30), expectation)
semicore_tops = {key: top for key, top in tops.items()
if top < vbm - SEMICORE_DEPTH_EV}
else:
for column in ("left", "right"):
occupied = [lv for lv in levels[column] if lv.electrons > 0]
if not occupied:
continue
homo_fragment = max(lv.energy for lv in occupied)
tops_column: dict[tuple[str, str], float] = {}
for lv in occupied:
parts = lv.label.split()
if len(parts) >= 3 and (parts[0], parts[1]) in spec_ranges:
key = (parts[0], parts[1])
tops_column[key] = max(tops_column.get(key, -1e30),
lv.energy)
semicore_tops.update({
key: top for key, top in tops_column.items()
if top < homo_fragment - SEMICORE_DEPTH_EV})
for level in levels["mo"]:
dropped: list = []
excluded = []
if semicore_tops and hamiltonian is not None:
for key, top in sorted(semicore_tops.items()):
if abs(level.energy - top) <= SEMICORE_DEPTH_EV:
continue # resonant filled-filled pair: genuine physics
indices = np.asarray(spec_ranges[key])
dropped.extend(indices.tolist())
excluded.append(" ".join(key))
elif semicore_tops:
if sqrt_overlap is not None:
gross = (np.abs(sqrt_overlap @ level.vectors) ** 2).sum(axis=1)
else:
gross = (np.conj(level.vectors)
* (overlap @ level.vectors)).real.sum(axis=1)
total = float(gross.sum()) or 1.0
for key in sorted(semicore_tops):
indices = np.asarray(spec_ranges[key])
if float(gross[indices].sum()) < 0.4 * total:
dropped.extend(indices.tolist())
excluded.append(" ".join(key))
keep_left = np.setdiff1d(np.asarray(rows_left), dropped)
keep_right = np.setdiff1d(np.asarray(rows_right), dropped)
cross = (level.vectors[keep_left, :].conj().T
@ overlap[np.ix_(keep_left, keep_right)]
@ level.vectors[keep_right, :])
population = 2.0 * float(np.trace(cross).real) / level.degeneracy
level.overlap_population = population
if population > BOND_CHARACTER_TOL:
level.bond_character = "bonding"
elif population < -BOND_CHARACTER_TOL:
level.bond_character = "antibonding"
else:
level.bond_character = "nonbonding"
level.detail += (f"\n{level.bond_character}: left-right "
"overlap population 2 Re c_L+ S c_R = "
f"{population:+.3f}"
+ (f" (semicore {', '.join(excluded)} tails "
"excluded)" if excluded else ""))
def assign_fragment_compositions(levels, overlap) -> None:
"""Crystal-level compositions in the LOEWDIN-ORTHOGONALIZED fragment
basis (sets level.composition, normalized, and
level.absolute_composition, the raw orthogonal weights).
The fragment eigenstates of the two columns are mutually
non-orthogonal (a diffuse Sc 4p against a compact F 2s combination
overlaps by 0.69 in ScF3), so plain projections |<phi|S|psi>|^2
double-count shared density: the occupied F 2s band state summed to
144% before renormalization and displayed as "Sc 4p 34%", and diffuse
virtuals (F 3s at +48 eV) showed up at 9% inside the semicore bands.
Symmetric (Loewdin) orthogonalization of the fragment-level basis --
the orthonormal set closest to the original fragment states, so the
labels keep their meaning -- removes the double counting: the weights
are |(G^{-1/2} Phi+ S psi)|^2 with G = Phi+ S Phi, they sum to the
span completeness (<= 100%), and the ScF3 examples become
F 2s 77 / Sc 4p 16 / Sc 3p 7 and F 3d 0.9 -> 0.1 for the t2g LUMO.
"""
fragments = levels["left"] + levels["right"]
for crystal in levels["mo"]:
crystal.composition = []
crystal.absolute_composition = []
if not fragments:
return
basis = np.hstack([fragment.vectors for fragment in fragments])
gram = basis.conj().T @ overlap @ basis
values, vectors = np.linalg.eigh(gram)
keep = values > 1e-8 * float(values.max())
inverse_half = ((vectors[:, keep] / np.sqrt(values[keep]))
@ vectors[:, keep].conj().T)
for crystal in levels["mo"]:
projections = basis.conj().T @ (overlap @ crystal.vectors)
tilde = inverse_half @ projections
weights = []
raw_weights = []
start = 0
for fragment in fragments:
count = fragment.vectors.shape[1]
rows = slice(start, start + count)
weight = float(np.sum(np.abs(tilde[rows]) ** 2))
weight /= crystal.degeneracy
raw = float(np.sum(np.abs(projections[rows]) ** 2))
raw /= crystal.degeneracy
start += count
if weight > 1e-6:
weights.append((fragment.level_id, weight))
if raw > 1e-6:
raw_weights.append((fragment.level_id, raw))
total = sum(weight for _, weight in weights) or 1.0
crystal.composition = [(i, w / total) for i, w in weights]
crystal.composition_completeness = float(
sum(weight for _, weight in weights))
# the RAW (double-counting) projections keep serving the alignment
# anchors, whose purity criteria were validated on them
crystal.absolute_composition = raw_weights
def _composition_string(symbols: list[str]) -> str:
counts: dict[str, int] = {}
for symbol in symbols:
counts[symbol] = counts.get(symbol, 0) + 1
return "".join(
f"{element}{count if count > 1 else ''}" for element, count in counts.items()
)
[docs]
class CrystalOrbitalDiagram:
"""Crystal-orbital diagram engine: symmetry + extended-Hueckel overlaps.
The engine behind ``crystod --diagram -c POSCAR --co-left A --co-right B``.
The two fragment sublattices are given as element formulas (every atom
of the primitive cell must belong to exactly one of them), each atom
carries its full core + valence shell basis, and at every high-symmetry
k point the fragment Bloch orbitals are symmetry-adapted, the Bloch
overlap lattice sums are evaluated, and the Wolfsberg-Helmholz
eigenproblem is solved for the two fragment columns and the crystal
column of the diagram (the module docstring describes the physics and
the caveats). The CLI sequence, :func:`report_and_write`, is
:meth:`special_kpoints`, then :meth:`solve_at` per k point, then the
HTML writer, which draws the hover wave-function sketches with
:meth:`supercell_for` and :meth:`sketch_partners`.
Args:
cell: The crystal structure as ``phonopy.structure.atoms.PhonopyAtoms``
(converted to the spglib primitive cell).
left_tokens: ``--co-left`` formula tokens, e.g. ``["SrTi"]`` or
``["Sr", "Ti"]``; a count such as ``O3`` is optional and, when
given, checked against the primitive cell.
right_tokens: ``--co-right`` formula tokens, e.g. ``["O3"]``.
symprec: Symmetry tolerance handed to spglib.
electrons: Electrons per primitive cell for the aufbau filling of
the crystal column (default: all electrons of the neutral atoms).
sketch_tokens: ``El-shell`` tokens (``"Ti-3d"``, ``"O_2p"``) that
restrict the drawn sketch components; ``None`` draws every
component (what the CLI does).
oxidation: ``{element: formal charge}`` for the point-charge lattice
of the removed sublattice; must be charge-neutral over the cell.
Default: pymatgen's oxidation-state guess.
conventional: Draw the hover sketches in the conventional cell
instead of the k-commensurate primitive supercell (display only).
Attributes:
builder: The :class:`SymmetryAdaptedOrbitalBasis` of the cell
(symmetry operations, irrep labels, ``spglib_dataset``).
symbols: Chemical symbols of the primitive-cell atoms.
positions: Their fractional coordinates, shape ``(n_atoms, 3)``.
lattice: Primitive lattice vectors as rows, in Angstrom.
formula: ``{"left": ..., "right": ...}``, the fragment formulas.
oxidation: The formal charges in use.
specs: One ``SublatticeSpec`` per (element, shell) block of the AO
basis, fragment-major; ``side_specs[column]`` lists those of
one fragment.
n_ao: Size of the AO basis; ``side_slice[column]`` is the slice of
a fragment's contiguous block.
orbitals: One ``AtomicOrbital`` per basis function, in the
representation order (spec-major, site-major, then ``m``).
electrons: Electrons per cell in the crystal column;
``side_electrons[column]`` those of the fragment columns.
h_raw: Bare atomic on-site energies (eV) of every basis function;
``h_bar`` the same after the shell-averaged point-charge shift,
``v_onsite`` the full on-site ligand-field matrix, and
``site_potential`` the Ewald monopole potential at every atom.
sketch_specs: Indices into ``specs`` of the drawn shells, or
``None`` for all.
last_estimated: Set by :meth:`solve_at`: the number of
near-dependent Bloch combinations whose energies are first-order
Loewdin estimates (marked ``~`` in the report), and
``last_dependent`` the number of linearly dependent
combinations removed.
Raises:
SystemExit: A fragment formula that does not match the cell, an
element without extended-Hueckel parameters or archived atomic
levels, or oxidation states that are not charge-neutral.
Example:
>>> from phonopy.interface.calculator import read_crystal_structure
>>> from crystod import salc
>>> from crystod.examples import example_path
>>> cell, _ = read_crystal_structure(
... str(example_path("221_PPOSCAR_ScF3")), interface_mode="vasp")
>>> diagram = salc.CrystalOrbitalDiagram(cell, ["Sc"], ["F3"])
>>> diagram.formula, diagram.oxidation
({'left': 'Sc', 'right': 'F3'}, {'Sc': 3.0, 'F': -1.0})
>>> levels, labels = diagram.solve_at([0, 0, 0])
>>> occupied = [lv for lv in levels["mo"] if lv.electrons]
>>> homo = max(occupied, key=lambda lv: lv.energy)
>>> homo.label, round(homo.energy, 2), homo.bond_character
('GM5- #1', -17.4, 'nonbonding')
"""
def __init__(self, cell, left_tokens: list[str], right_tokens: list[str],
symprec: float = 1e-5, electrons: float | None = None,
sketch_tokens: list[str] | None = None,
oxidation: dict[str, float] | None = None,
conventional: bool = False):
# --conventional: draw the hover sketches in the conventional cell
# (display-only; the diagram itself is unchanged)
self.conventional = bool(conventional)
self.builder = SymmetryAdaptedOrbitalBasis(cell=cell, symprec=symprec)
primitive = self.builder.primitive_cell
self.symbols = get_chemical_symbols(primitive)
self.positions = np.array(get_scaled_positions(primitive))
self.lattice = np.array(primitive.cell) # rows, Angstrom
formulas = {
"left": parse_fragment_formula(left_tokens, "--co-left"),
"right": parse_fragment_formula(right_tokens, "--co-right"),
}
self.formula = {
"left": "".join(left_tokens),
"right": "".join(right_tokens),
}
composition: dict[str, int] = {}
for symbol in self.symbols:
composition[symbol] = composition.get(symbol, 0) + 1
comp_str = _composition_string(self.symbols)
flags = {"left": "--co-left", "right": "--co-right"}
assigned: dict[str, str] = {}
for side in ("left", "right"):
for element, count in formulas[side]:
if element in assigned:
raise SystemExit(
f"ERROR: element {element} is listed more than once "
"across --co-left/--co-right."
)
if element not in composition:
raise SystemExit(
f"ERROR: element {element} is not in the structure "
f"(primitive-cell composition: {comp_str})."
)
if count is not None and count != composition[element]:
raise SystemExit(
f"ERROR: {flags[side]} lists {element}{count} but the "
f"primitive cell has {composition[element]} "
f"{element} atom(s) (composition: {comp_str})."
)
if element not in EHT_PARAMETERS:
raise SystemExit(
f"ERROR: no extended-Hueckel parameters for element "
f"{element} (supported: {', '.join(EHT_PARAMETERS)})."
)
assigned[element] = side
missing = [el for el in composition if el not in assigned]
if missing:
raise SystemExit(
f"ERROR: element(s) {', '.join(missing)} not assigned to "
f"--co-left/--co-right (primitive-cell composition: "
f"{comp_str}; every atom must belong to one fragment)."
)
# formal oxidation states of the ions: the removed sublattice enters
# each fragment as a point-charge lattice with these charges
if oxidation is None:
from pymatgen.core import Composition
guesses = Composition(comp_str).oxi_state_guesses()
if not guesses:
raise SystemExit(
"ERROR: could not guess the oxidation states of "
f"{comp_str}; pass them explicitly, e.g. "
"--oxidation Sr=+2 Ti=+4 O=-2."
)
oxidation = {el: float(q) for el, q in guesses[0].items()}
missing_ox = [el for el in composition if el not in oxidation]
if missing_ox:
raise SystemExit(
f"ERROR: --oxidation misses element(s) "
f"{', '.join(missing_ox)}."
)
net = sum(oxidation[el] * composition[el] for el in composition)
if abs(net) > 1e-6:
raise SystemExit(
f"ERROR: oxidation states are not charge-neutral "
f"(net {net:+g} per cell)."
)
self.oxidation = oxidation
# full-electron basis (core + valence shells of every atom),
# fragment-major so each fragment is one contiguous AO block; core
# shells get Slater-rule exponents and the archived PySCF
# neutral-atom Hartree-Fock levels (reference/atomic_level_*),
# valence shells the extended-Hueckel parameters. Shells frozen
# into the def2 effective core potential (beyond Kr) do not exist
# in the atomic data and are omitted -- the pseudopotential
# picture.
self.specs: list[SublatticeSpec] = []
self.side_specs = {"left": [], "right": []}
offset = 0
for side in ("left", "right"):
for element, _count in formulas[side]:
if element not in ATOMIC_LEVELS:
raise SystemExit(
f"ERROR: no archived atomic levels for {element}; "
"run script/generate_atomic_levels.py and "
"script/collect_atomic_levels.py."
)
sites = [
index for index, symbol in enumerate(self.symbols)
if symbol == element
]
levels = ATOMIC_LEVELS[element]["levels"]
shells = [
(shell, int(shell[0]), "spdf".index(shell[-1]),
slater_zeta(element, shell), levels[shell])
for shell in CORE_SHELLS[element] if shell in levels
] + list(EHT_PARAMETERS[element])
for shell, n, l, zeta, h_ii in shells:
spec = SublatticeSpec(element, shell[-1], shell, n, l,
zeta, h_ii, sites, side, offset)
self.specs.append(spec)
self.side_specs[side].append(spec)
offset += spec.n_ao
self.n_ao = offset
self.side_slice = {}
start = 0
for side in ("left", "right"):
width = sum(spec.n_ao for spec in self.side_specs[side])
self.side_slice[side] = slice(start, start + width)
start += width
# AO objects in the representation ordering (spec-major, site-major)
self.orbitals: list[AtomicOrbital] = []
for spec in self.specs:
for site in spec.sites:
for m in range(2 * spec.l + 1):
self.orbitals.append(AtomicOrbital(
site, spec.element, spec.shell, spec.n, spec.l, m,
spec.zeta, spec.h_ii,
))
self.side_electrons = {
side: sum(
ATOMIC_LEVELS[element]["electrons"] * composition[element]
for element, _count in formulas[side]
)
for side in ("left", "right")
}
if electrons is None:
electrons = sum(self.side_electrons.values())
self.electrons = float(electrons)
# --atomic-orbital: the sketch filter (indices into self.specs)
self.sketch_specs: frozenset[int] | None = None
self.sketch_tokens = sketch_tokens
if sketch_tokens is not None:
available = []
for spec in self.specs:
token = f"{spec.element}-{spec.shell}"
if token not in available:
available.append(token)
selected = set()
for element, shell in parse_sketch_tokens(sketch_tokens):
hits = [
index for index, spec in enumerate(self.specs)
if spec.element == element
and (spec.shell == shell
or (len(shell) == 1 and spec.letter == shell))
]
if not hits:
raise SystemExit(
f"ERROR: --atomic-orbital {element}-{shell} matches no "
f"basis orbital (available: {', '.join(available)})."
)
selected.update(hits)
self.sketch_specs = frozenset(selected)
self._aligned = make_aligned_cache()
self._cutoffs = self._pair_cutoffs()
self._images = self._lattice_images(max(self._cutoffs.values()))
self._build_ligand_field()
# -------------------------------------------------- point-charge field
def _build_ligand_field(self, r_cut_angstrom: float = 7.0):
"""Same-site matrices of the removed-sublattice point-charge field.
Every atom feels the complementary sublattice as a lattice of point
charges with the formal oxidation states: charges within
r_cut_angstrom enter as exact <phi_i|q/|r-R||phi_j> STO integrals
(monopole shift + multipole ligand-field splitting), the long-range
rest as the Ewald site potential (neutralizing-background
convention, as in charged periodic DFT cells). The blocks are
added identically to the fragment and the crystal Hamiltonians, so
the three diagram columns share one energy reference."""
complementary = {"left": "right", "right": "left"}
side_of_atom = {}
for side in ("left", "right"):
for spec in self.side_specs[side]:
for site in spec.sites:
side_of_atom[site] = side
charge_lattice = {
side: [
(self.oxidation[self.symbols[site]], self.positions[site])
for site in range(len(self.symbols))
if side_of_atom[site] == side
]
for side in ("left", "right")
}
bounds = []
volume = abs(np.linalg.det(self.lattice))
for i in range(3):
j, k = (i + 1) % 3, (i + 2) % 3
perpendicular = volume / np.linalg.norm(
np.cross(self.lattice[j], self.lattice[k])
)
bounds.append(int(np.ceil(r_cut_angstrom / perpendicular)) + 1)
images = [
np.array([n1, n2, n3])
for n1 in range(-bounds[0], bounds[0] + 1)
for n2 in range(-bounds[1], bounds[1] + 1)
for n3 in range(-bounds[2], bounds[2] + 1)
]
self.h_raw = np.array([o.h_ii for o in self.orbitals], dtype=float)
self.h_bar = self.h_raw.copy()
self.v_onsite = np.zeros((self.n_ao, self.n_ao))
self.site_potential = np.zeros(len(self.symbols))
for atom in range(len(self.symbols)):
charges = charge_lattice[complementary[side_of_atom[atom]]]
near = []
near_monopole = 0.0
for q, frac in charges:
for image in images:
vector = (frac + image - self.positions[atom]) @ self.lattice
distance = float(np.linalg.norm(vector))
if distance <= r_cut_angstrom:
near.append((q, vector * ANGSTROM_TO_BOHR))
near_monopole += q / distance
v_ewald = ewald_site_potential(
self.lattice, charges, self.positions[atom]
)
self.site_potential[atom] = -COULOMB_EV_ANGSTROM * v_ewald
e_far = -COULOMB_EV_ANGSTROM * (v_ewald - near_monopole)
indices = [
index for index, orbital in enumerate(self.orbitals)
if orbital.atom == atom
]
orbitals = [self.orbitals[index] for index in indices]
block = point_charge_block(orbitals, near)
# The EHT basis treats same-site shells of one l (2p/3p/4p,
# ...) as orthonormal, but the raw STO radials are not: the
# point-charge matrix must be expressed in the Loewdin-
# orthogonalized on-site basis, so that a constant potential
# maps exactly onto the identity (and the far-field term is
# exactly diagonal).
s_site = np.eye(len(indices))
for a_local, oa in enumerate(orbitals):
for b_local in range(a_local + 1, len(indices)):
ob = orbitals[b_local]
if (oa.l, oa.m) == (ob.l, ob.m) and oa.shell != ob.shell:
s_site[a_local, b_local] = s_site[b_local, a_local] \
= radial_overlap((oa.n, oa.zeta), (ob.n, ob.zeta))
s_values, s_vectors = np.linalg.eigh(s_site)
o_half = (s_vectors / np.sqrt(s_values)) @ s_vectors.T
block = o_half @ block @ o_half
block += np.eye(len(indices)) * e_far
# The monopole (Madelung) part of the site shift is omitted: it
# depends on the neutralizing-background convention of the
# charged sublattice array (it flips the SrTiO3 band ordering),
# and it largely cancels against the intra-atomic charging
# energy not present in the extended-Hueckel VSIPs -- the
# standard argument why neutral-atom VSIPs work in ionic
# crystals. What remains is background-independent and
# absolutely convergent: the anisotropic multipole ligand field
# (t2g/eg splittings, ...) and the near-shell penetration
# corrections. The omitted jellium-referenced monopole is kept
# in self.site_potential for the report.
block -= np.eye(len(indices)) * self.site_potential[atom]
self.v_onsite[np.ix_(indices, indices)] = block
# shell-averaged (rotation-invariant) shift for the W-H h_bar
shells = {}
for local, orbital in enumerate(orbitals):
shells.setdefault(orbital.shell, []).append(local)
for members in shells.values():
average = float(np.mean([block[m, m] for m in members]))
for m in members:
self.h_bar[indices[m]] += average
# ------------------------------------------------------------ Bloch sums
def _pair_cutoffs(self, tol: float = 2e-5, r_max: float = 44.0):
"""Per shell-pair lattice-sum cutoff from the actual STO tails.
Diffuse cation shells (e.g. Sr 5s, Ti 4p) still overlap at 30+ bohr;
a fixed short cutoff truncates their Bloch sums so badly that S_k
loses positive semidefiniteness. The sigma overlap along z is
probed outward until it falls below tol."""
sigma_m = {0: 0, 1: 2, 2: 2} # s / pz / dz2: the slowest-decaying
kinds = []
seen = set()
for spec in self.specs:
key = (spec.element, spec.shell)
if key not in seen:
seen.add(key)
kinds.append(spec)
cutoffs = {}
for index, a_spec in enumerate(kinds):
for b_spec in kinds[index:]:
key_a = (a_spec.element, a_spec.shell)
key_b = (b_spec.element, b_spec.shell)
if a_spec.l > 2 or b_spec.l > 2:
# f shells appear only as ultra-compact cores (Pb/Bi
# 4f); their inter-site overlap is neglected
cutoffs[key_a, key_b] = cutoffs[key_b, key_a] = 6.0
continue
a = AtomicOrbital(0, a_spec.element, a_spec.shell, a_spec.n,
a_spec.l, sigma_m[a_spec.l], a_spec.zeta,
a_spec.h_ii)
b = AtomicOrbital(0, b_spec.element, b_spec.shell, b_spec.n,
b_spec.l, sigma_m[b_spec.l], b_spec.zeta,
b_spec.h_ii)
r = 6.0
while r < r_max and abs(pair_overlap(
a, b, np.array([0.0, 0.0, r]), self._aligned
)) > tol:
r += 2.0
cutoffs[key_a, key_b] = r
cutoffs[key_b, key_a] = r
return cutoffs
def _lattice_images(self, cutoff_bohr: float):
"""Integer lattice translations n with any-atom pair within cutoff."""
lattice_bohr = self.lattice * ANGSTROM_TO_BOHR
volume = abs(np.linalg.det(lattice_bohr))
bounds = []
for i in range(3):
j, k = (i + 1) % 3, (i + 2) % 3
perpendicular = volume / np.linalg.norm(
np.cross(lattice_bohr[j], lattice_bohr[k])
)
# +1: margin for the in-cell offset x_j - x_i
bounds.append(int(np.ceil(cutoff_bohr / perpendicular)) + 1)
return np.array([
[n1, n2, n3]
for n1 in range(-bounds[0], bounds[0] + 1)
for n2 in range(-bounds[1], bounds[1] + 1)
for n3 in range(-bounds[2], bounds[2] + 1)
])
[docs]
def bloch_overlap(self, kpoint) -> np.ndarray:
"""Bloch overlap matrix ``S(k)`` of the full AO basis, atom gauge.
``S_k(i, j) = sum_n exp(2 pi i k . (n + x_j - x_i)) s(i at 0, j at
n)`` over the lattice translations ``n`` within the pair cutoffs;
the on-site term of an orbital with itself is 1, and the compact
``f`` cores carry no inter-site overlap.
Args:
kpoint: Three primitive reciprocal coordinates.
Returns:
Hermitian complex array of shape ``(n_ao, n_ao)``.
"""
k = np.asarray(kpoint, dtype=float)
S = np.zeros((self.n_ao, self.n_ao), dtype=complex)
lattice_bohr = self.lattice * ANGSTROM_TO_BOHR
for i in range(self.n_ao):
for j in range(i, self.n_ao):
a, b = self.orbitals[i], self.orbitals[j]
cutoff = self._cutoffs[
(a.element, a.shell), (b.element, b.shell)
]
offset = self.positions[b.atom] - self.positions[a.atom]
fractionals = offset + self._images
distances = np.linalg.norm(fractionals @ lattice_bohr, axis=1)
total = 0.0 + 0.0j
for image_index in np.nonzero(distances <= cutoff)[0]:
fractional = fractionals[image_index]
if distances[image_index] < 1e-9:
overlap = 1.0 if (a.shell, a.m) == (b.shell, b.m) else 0.0
elif a.l > 2 or b.l > 2:
overlap = 0.0 # compact f cores: no inter-site overlap
else:
overlap = pair_overlap(
a, b, fractional @ lattice_bohr, self._aligned
)
if overlap == 0.0:
continue
phase = np.exp(2j * np.pi * float(k @ fractional))
total += phase * overlap
S[i, j] = total
S[j, i] = np.conj(total)
return S
[docs]
def hamiltonian(self, S: np.ndarray) -> np.ndarray:
"""Wolfsberg-Helmholz Hamiltonian over the Bloch overlaps.
``H_ij = K S_ij (h_i + h_j) / 2`` with the shell-averaged shifted
energies ``h_bar`` (rotation-invariant, so the symmetry of ``H``
stays exact); the on-site blocks carry the bare atomic energies plus
the full anisotropic point-charge ligand-field matrices
(``v_onsite``). The diagonal of ``S_k`` is 1 plus the same-orbital
neighbour-cell Bloch sum, so only the on-site ``R = 0`` term is the
bare atomic energy: the diagonal correction ``(1 - K)`` restores
``h + K h_bar (S_kk - 1) + V_ii``.
Args:
S: The Bloch overlap matrix from :meth:`bloch_overlap`.
Returns:
Hermitian complex array of shape ``(n_ao, n_ao)``, in eV.
"""
H = 0.5 * WOLFSBERG_HELMHOLZ_K * (
self.h_bar[:, None] + self.h_bar[None, :]
) * S
H += self.v_onsite
H[np.diag_indices(self.n_ao)] += (
self.h_raw - WOLFSBERG_HELMHOLZ_K * self.h_bar
)
return H
# -------------------------------------------------------- representation
[docs]
def little_group_data(self, kpoint):
"""Irreps, labels and the AO representation of the little group at ``k``.
Args:
kpoint: Three primitive reciprocal coordinates.
Returns:
``(irreps, mapping, labels, representation)``: the spgrep irreps,
the indices of the little-group operations into
``builder.rotations``, their ISO-IR labels, and one complex
``(n_ao, n_ao)`` matrix per operation -- block-diagonal over the
(element, shell) specs, each block the Kronecker product of the
Bloch-phased site permutation with the real-orbital Wigner
matrix of the shell.
"""
irreps, mapping = get_spacegroup_irreps_from_primitive_symmetry(
rotations=self.builder.rotations,
translations=self.builder.translations,
kpoint=kpoint,
)
labels = self.builder.get_irrep_labels(kpoint, irreps, mapping)
little_rotations = self.builder.rotations[mapping]
little_translations = self.builder.translations[mapping]
permutations = self.builder.get_permutation_reps_at_k(
little_rotations=little_rotations,
little_translations=little_translations,
kpoint=kpoint,
)
representation = []
for op_index, op in enumerate(mapping):
blocks = []
wigners = {}
for spec in self.specs:
if spec.l not in wigners:
wigners[spec.l] = wigner_D_real(
spec.l,
np.real(self.builder.rotations_cartesian[op]),
)
grid = np.ix_(spec.sites, spec.sites)
blocks.append(np.kron(
permutations[op_index][grid], wigners[spec.l]
))
matrix = np.zeros((self.n_ao, self.n_ao), dtype=complex)
row = 0
for block in blocks:
size = block.shape[0]
matrix[row:row + size, row:row + size] = block
row += size
representation.append(matrix)
return irreps, mapping, labels, representation
# -------------------------------------------------------------- solving
@staticmethod
def _generalized_eigh(H: np.ndarray, S: np.ndarray, h_bar: np.ndarray):
"""Canonically orthogonalized generalized eigenproblem (complex
Hermitian). Near-dependent Bloch combinations (overlap eigenvalue
below _OVERLAP_FLOOR, see there) are excluded from the variational
solve -- their EHT energies diverge as (1-K) h / eigenvalue (the
overlap catastrophe) -- but they are genuine Bloch states, so they
are solved separately with the first-order Loewdin-orthogonalized
Hamiltonian H~ = H - (S - I) (h_bar_i + h_bar_j)/2 (the same
correction as the coupling tables' |H~|; h_bar keeps the symmetry
of H~ exact). H~ is bounded for any overlap and reduces to the
variational result at small overlap; for a single shell's Bloch
sum it gives E = h (1 + (K-1) sigma) instead of the divergent
h (K + (1-K)/eps). Returns (energies, vectors,
(estimated energies, estimated plain vectors), n_dependent) where
n_dependent counts truly linearly dependent combinations
(eigenvalue below _DEPENDENT_TOL) that are removed outright."""
s_values, s_vectors = np.linalg.eigh(S)
keep = s_values > _OVERLAP_FLOOR
X = s_vectors[:, keep] / np.sqrt(s_values[keep])
H_orth = X.conj().T @ H @ X
energies, coefficients = np.linalg.eigh(H_orth)
estimate = (~keep) & (s_values > _DEPENDENT_TOL)
est_energies = np.empty(0)
est_vectors = np.empty((S.shape[0], 0), dtype=complex)
if np.any(estimate):
V = s_vectors[:, estimate]
mean = 0.5 * (h_bar[:, None] + h_bar[None, :])
H_first = H - (S - np.eye(S.shape[0])) * mean
block = V.conj().T @ H_first @ V
est_energies, w = np.linalg.eigh(block)
est_vectors = V @ w
return (energies, X @ coefficients, (est_energies, est_vectors),
int(np.sum(s_values <= _DEPENDENT_TOL)))
def _group_levels(self, energies, vectors):
"""Cluster eigenvalues into degenerate groups."""
groups = []
start = 0
for i in range(1, len(energies) + 1):
if i == len(energies) or energies[i] - energies[start] > max(
_DEGENERACY_TOL, 1e-6 * max(1.0, abs(energies[start]))
):
groups.append((float(np.mean(energies[start:i])),
vectors[:, start:i]))
start = i
return groups
def _irrep_split(self, vectors, S, representation, irreps, labels,
subspace=None):
"""Split an S-orthonormal degenerate group into its irrep components.
Robust against accidental degeneracies (several irreps at one
energy, e.g. nearly uncoupled sublattice shells): the group space is
decomposed with the character projectors and one (label, vectors)
entry is returned per contributing irrep."""
if subspace is not None:
embedded = np.zeros((self.n_ao, vectors.shape[1]), dtype=complex)
embedded[subspace] = vectors
vectors = embedded
gram = vectors.conj().T @ S @ vectors
values, basis = np.linalg.eigh(gram)
V = vectors @ (basis / np.sqrt(values))
characters = np.array([
np.trace(V.conj().T @ S @ D @ V) for D in representation
])
order = len(representation)
components = []
for irrep, label in zip(irreps, labels):
chi = np.array(get_character(irrep), dtype=complex)
multiplicity = float(np.real(
np.sum(characters * np.conj(chi)) / order
))
count = int(round(multiplicity))
if count <= 0:
continue
dimension = irrep.shape[1]
projector = np.zeros((self.n_ao, self.n_ao), dtype=complex)
for g, D in enumerate(representation):
projector += np.conj(chi[g]) * D
projector *= dimension / order
projected = projector @ V
# S-orthonormal basis of the projected span
gram_p = projected.conj().T @ S @ projected
p_values, p_basis = np.linalg.eigh(gram_p)
keep = p_values > 1e-6
space = projected @ (p_basis[:, keep] / np.sqrt(p_values[keep]))
if space.shape[1] != count * dimension:
# numerical safety: fall back to the expected count
space = space[:, : count * dimension]
components.append((label, space))
if not components:
components.append(("?", V))
return components
def _dominant_spec(self, space, S, column) -> SublatticeSpec:
"""Fragment (element, shell) block with the largest Mulliken gross
population of an S-orthonormal level space (for the level label)."""
Sv = S @ space
weights = []
for spec in self.side_specs[column]:
rows = slice(spec.offset, spec.offset + spec.n_ao)
weights.append(float(np.real(
np.sum(np.conj(space[rows]) * Sv[rows])
)))
return self.side_specs[column][int(np.argmax(weights))]
[docs]
def solve_at(self, kpoint):
"""Solve the fragment and crystal eigenproblems at one k point.
The fragment columns are the generalized eigenproblems of the two
sublattice blocks, the crystal column that of the full basis. Every
level is labelled by projecting its degenerate space onto the irreps,
filled by aufbau, and given its Loewdin fragment composition and
its COOP bond character (:func:`assign_bond_characters`); the
outermost columns list the isolated on-site shell levels. Sets
:attr:`last_estimated` and :attr:`last_dependent`.
Args:
kpoint: Three primitive reciprocal coordinates, for the diagram
one of :meth:`special_kpoints`.
Returns:
``(levels, labels)``: ``levels`` maps the columns ``"left"``,
``"mo"`` (the crystal), ``"right"``, ``"left-ao"`` and
``"right-ao"`` (the outermost isolated-shell columns) to lists
of ``DiagramLevel`` records -- ``label``, ``irrep``, ``energy``
(eV), ``degeneracy``, ``electrons``, ``vectors`` (shape
``(n_ao, degeneracy)``, S-orthonormal), ``composition``,
``bond_character``, ``overlap_population``, ``estimated`` and
``detail`` -- and ``labels`` are the ISO-IR irrep labels at
``kpoint``.
Raises:
SystemExit: The representation does not leave ``S`` or ``H``
invariant (a gauge or real-harmonics inconsistency; please
report the case).
"""
S = self.bloch_overlap(kpoint)
H = self.hamiltonian(S)
irreps, mapping, labels, representation = self.little_group_data(kpoint)
# gauge self-check: the representation must leave S invariant
worst = max(
float(np.max(np.abs(D.conj().T @ S @ D - S)))
for D in representation
)
if worst > 1e-5:
raise SystemExit(
"ERROR: Bloch-overlap gauge inconsistency "
f"(residual {worst:.2e}); please report this case."
)
# ... and H (checks the point-charge ligand-field blocks against
# the site-symmetry representation, i.e. the real-harmonics
# conventions)
h_scale = float(np.max(np.abs(H)))
worst_h = max(
float(np.max(np.abs(D.conj().T @ H @ D - H)))
for D in representation
) / max(h_scale, 1.0)
if worst_h > 1e-6:
raise SystemExit(
"ERROR: Hamiltonian symmetry inconsistency "
f"(relative residual {worst_h:.2e}); please report this case."
)
def strip(label):
return label.split("(")[0]
def merged_groups(energies, vectors, est_energies, est_vectors,
S_block):
"""Exact and first-order-estimated level groups, energy-sorted.
The third entry is None for variational levels and the mean
overlap eigenvalue of the group for estimated ones (for the
annotation)."""
groups = [(energy, group, None)
for energy, group in self._group_levels(energies,
vectors)]
for energy, group in self._group_levels(est_energies,
est_vectors):
overlap_eigenvalue = float(np.mean(np.real(np.sum(
np.conj(group) * (S_block @ group), axis=0))))
groups.append((energy, group, overlap_eigenvalue))
groups.sort(key=lambda item: item[0])
return groups
levels = {"left": [], "mo": [], "right": []}
self.last_estimated = 0
self.last_dependent = 0
# fragment (sublattice) levels: the full valence problem of one side
for column in ("left", "right"):
block = self.side_slice[column]
indices = np.arange(block.start, block.stop)
energies, vectors, (est_energies, est_vectors), dependent = \
self._generalized_eigh(
H[block, block], S[block, block], self.h_bar[block]
)
self.last_estimated += est_energies.size
self.last_dependent += dependent
for energy, group, overlap_eigenvalue in merged_groups(
energies, vectors, est_energies, est_vectors, S[block, block]
):
for irrep_label, space in self._irrep_split(
group, S, representation, irreps, labels, subspace=indices
):
name = strip(irrep_label)
spec = self._dominant_spec(space, S, column)
level = DiagramLevel(
level_id=f"{column}{len(levels[column])}",
column=column,
energy=float(energy),
degeneracy=space.shape[1],
irrep=name,
label=f"{spec.element} {spec.shell} {name}",
vectors=space,
)
if overlap_eigenvalue is not None:
level.estimated = True
level.detail = _ESTIMATED_NOTE.format(
eps=overlap_eigenvalue, floor=_OVERLAP_FLOOR)
levels[column].append(level)
# two fragment levels can share (element, shell, irrep) -- e.g.
# the two F 2p GM4- combinations; number them so the crystal
# compositions stay readable
seen: dict[str, int] = {}
for level in levels[column]:
seen[level.label] = seen.get(level.label, 0) + 1
repeated = {label for label, n in seen.items() if n > 1}
occurrence: dict[str, int] = {}
for level in levels[column]:
if level.label in repeated:
occurrence[level.label] = occurrence.get(level.label, 0) + 1
level.label = f"{level.label}#{occurrence[level.label]}"
# crystal levels
energies, vectors, (est_energies, est_vectors), dependent = \
self._generalized_eigh(H, S, self.h_bar)
self.last_estimated += est_energies.size
self.last_dependent += dependent
counts: dict[str, int] = {}
for energy, group, overlap_eigenvalue in merged_groups(
energies, vectors, est_energies, est_vectors, S
):
for irrep_label, space in self._irrep_split(
group, S, representation, irreps, labels
):
name = strip(irrep_label)
counts[name] = counts.get(name, 0) + 1
occurrence = counts[name]
level = DiagramLevel(
level_id=f"mo{len(levels['mo'])}",
column="mo",
energy=float(energy),
degeneracy=space.shape[1],
irrep=name,
# "GM4- #2" = second GM4- multiplet from the bottom, the
# same #N numbering as the fragment columns; "GM4-(2)"
# read like a degeneracy count
label=f"{name} #{occurrence}",
vectors=space,
)
if overlap_eigenvalue is not None:
level.estimated = True
level.detail = _ESTIMATED_NOTE.format(
eps=overlap_eigenvalue, floor=_OVERLAP_FLOOR)
levels["mo"].append(level)
# compositions: crystal levels in the Loewdin-orthogonalized
# fragment-level basis (plain |<phi|S|psi>|^2 double-counts the
# strongly overlapping fragment states; see the helper)
assign_fragment_compositions(levels, S)
# electron filling (aufbau, 2 electrons per orbital)
self._fill(levels["mo"], self.electrons)
for column in ("left", "right"):
self._fill(levels[column], self.side_electrons[column])
# COOP bonding character of every crystal level (needs occupations)
spec_ranges: dict[tuple[str, str], np.ndarray] = {}
for column in ("left", "right"):
for spec in self.side_specs[column]:
key = (spec.element, spec.shell)
indices = np.arange(spec.offset, spec.offset + spec.n_ao)
spec_ranges[key] = (np.concatenate([spec_ranges[key], indices])
if key in spec_ranges else indices)
# the ONE displayed composition (panel bars, hover tooltip and the
# terminal): Loewdin AO-shell populations of each crystal state --
# the same partial-charge measure as the PySCF engine. The fragment
# projection above keeps positioning the connector lines but its
# weights are not displayed (fragment eigenstates mix AO shells
# among themselves, so the two measures disagree).
eigenvalues, eigenvectors = np.linalg.eigh(S)
sqrt_overlap = (eigenvectors * np.sqrt(np.clip(eigenvalues.real,
0.0, None))
) @ eigenvectors.conj().T
# fragment columns included: a repeated label like "F 2p GM4-#2"
# only names the dominant shell of a same-irrep mixture, and the
# bars/tooltip must reveal the mixture itself. Fragment levels
# list their own sublattice's shells only (they are sublattice
# states; the Loewdin attribution reaching the other side is
# aggregated into one closing note).
side_keys = {
column: {(spec.element, spec.shell)
for spec in self.side_specs[column]}
for column in ("left", "right")
}
# outermost columns (the MolOD "ligand-ao"/"center-ao" analogue):
# one level per (element, shell) at the rotation-invariant on-site
# energy h_bar (VSIP + spherical part of the point-charge ligand
# field), so the sublattice column reads as "how the equivalent
# atoms' Bloch combinations split each atomic shell at this k
# point" -- with several atoms per cell (wurtzite AlN: 2 Al) the
# intra-sublattice overlap splits one shell into several levels
# and can even push 3p combinations below 3s; the splitting
# connector lines carry the same Loewdin shell populations as the
# tooltip. k-independent by construction.
ao_ids: dict[str, dict[tuple[str, str], str]] = {}
for column in ("left", "right"):
ao_column = f"{column}-ao"
levels[ao_column] = []
ao_ids[column] = {}
for spec in self.side_specs[column]:
key = (spec.element, spec.shell)
if key in ao_ids[column]:
continue
count = len(spec.sites)
# h_bar carries a per-SITE ligand-field shift: for an element
# on several Wyckoff positions (K2SeO4: two K classes, 1 eV
# apart) the single AO level shows the mean, and the spread
# is disclosed instead of silently drawing the first site
block = self.h_bar[spec.offset:spec.offset + spec.n_ao]
onsite = float(np.mean(block))
raw_vsip = float(np.mean(
self.h_raw[spec.offset:spec.offset + spec.n_ao]))
field_shift = onsite - raw_vsip
spread = float(np.max(block) - np.min(block))
equivalent = self.builder.spglib_dataset["equivalent_atoms"]
orbits = len({int(equivalent[site]) for site in spec.sites})
prefix = f"{count}" if count > 1 else ""
level = DiagramLevel(
level_id=f"{ao_column}{len(levels[ao_column])}",
column=ao_column,
energy=onsite,
degeneracy=2 * spec.l + 1,
irrep="",
label=f"{prefix}{spec.element} {spec.shell}",
vectors=None,
)
level.electrons = None
level.display_composition = [
(f"{spec.element} {spec.shell}", 1.0)]
site_note = ""
if spread > 0.02:
site_note = (f" (mean over the {count} sites on {orbits} "
"crystallographically inequivalent "
f"positions; on-site spread {spread:.2f} eV)")
if count > 1:
inequivalent = (
"" if orbits == 1 else
f" ({orbits} inequivalent positions)")
tail = (f"the {self.formula[column]} column shows how "
f"the {count} {spec.element} atoms'{inequivalent}"
" Bloch combinations split this shell at each "
"k point")
else:
tail = (f"the {self.formula[column]} column shows this "
"shell's Bloch combination at each k point")
cage_note = ""
if abs(field_shift) > 8.0:
cage_note = (
"\nWARNING: a point-charge shift this large is no "
"small correction: the shell's STO is diffuse "
"enough to engulf the surrounding charge cage, so "
"this on-site energy -- and every level the shell "
"dominates -- is an artifact of the point-charge "
"model, not chemistry")
level.detail = (
f"isolated {spec.element} {spec.shell} on-site level: "
f"VSIP {raw_vsip:.2f} + point-charge ligand field "
f"(spherical part) {field_shift:+.2f} "
f"= {onsite:.2f} eV{site_note}\n{tail}{cage_note}")
ao_ids[column][key] = level.level_id
levels[ao_column].append(level)
channel_of_ao = np.empty(self.n_ao, dtype=object)
for spec in self.specs:
w = 2 * spec.l + 1
for site_pos, site in enumerate(spec.sites):
start = spec.offset + site_pos * w
for index in range(start, start + w):
channel_of_ao[index] = (site, spec.l)
for column in ("mo", "left", "right"):
for level in levels[column]:
gross = (np.abs(sqrt_overlap @ level.vectors) ** 2
).sum(axis=1) / level.degeneracy
# per-(atom, l) Loewdin channel populations: the hover
# sketch calibrates its lobe sizes to these -- raw
# coefficient x STO-amplitude lobes misstate the mix (a
# 90%-Ti-4s level used to be DRAWN as d lobes: the compact
# 3d weighs ~4x the diffuse 4s at the probe radius, and
# the semicore 3s tail cancels most of the s channel)
pops: dict = {}
for index, value in enumerate(gross):
key = channel_of_ao[index]
if key is not None:
pops[key] = pops.get(key, 0.0) + float(value)
level.channel_pop = pops
shares = []
linked = []
for (element, shell), indices in spec_ranges.items():
if column != "mo" and (element, shell) \
not in side_keys[column]:
continue
value = float(gross[indices].sum())
if value >= 0.001:
shares.append((value, f"{element} {shell}"))
if column != "mo":
linked.append(
(ao_ids[column][(element, shell)], value))
if column != "mo":
# splitting connector lines into the atomic-shell column
level.composition = linked
shares.sort(key=lambda item: -item[0])
level.display_composition = [
(f"{name} {level.irrep}", value)
for value, name in shares
]
row = "Loewdin: " + ", ".join(
f"{name} {100 * value:.1f}%" for value, name in shares)
if column != "mo":
remainder = 1.0 - sum(value for value, _ in shares)
if remainder >= 0.005:
row += (f" (+{100 * remainder:.1f}% Loewdin-"
"attributed to the other sublattice's "
"basis: overlap density; the state has no "
"coefficients there)")
level.detail = (row if not level.detail
else f"{row}\n{level.detail}")
assign_bond_characters(
levels, S,
np.arange(self.side_slice["left"].start,
self.side_slice["left"].stop),
np.arange(self.side_slice["right"].start,
self.side_slice["right"].stop),
spec_ranges,
hamiltonian=H,
)
return levels, labels
@staticmethod
def _fill(column_levels, electrons):
"""Aufbau filling (2 electrons per orbital).
Returns the level left PARTIALLY filled, if any: an insulator's
count exhausts exactly at a level boundary, so a partial level is
either a genuinely metallic k point or a level-ordering artifact
worth flagging (rutile TiO2 in the extended-Hueckel engine: the
near-dependent Ti 4p Bloch sums inherit a point-charge-shifted
on-site energy, land inside the O 2p band, swallow 4 electrons
and leave the true O 2p top half-filled)."""
remaining = float(electrons)
integer_total = abs(remaining - round(remaining)) < 1e-9
odd_total = integer_total and int(round(remaining)) % 2 == 1
partial = None
for level in sorted(column_levels, key=lambda lv: lv.energy):
capacity = 2 * level.degeneracy
want = min(remaining, capacity)
take = int(round(want))
level.electrons = max(take, 0)
fractional = want > 0 and abs(want - take) > 1e-6
if 0 < level.electrons < capacity or fractional:
partial = level
if odd_total:
# parity-forced: an odd electron count half-fills one
# level at EVERY k point -- not a warning, a property
level.partial_parity = True
note = ("NOTE: this column's electron count "
f"({int(round(electrons))}) is odd, so one "
"level is necessarily half-filled at every "
"k point (spin-restricted display)")
elif fractional:
level.partial_parity = True
note = ("NOTE: aufbau ends with a fractional "
f"occupation (requested {electrons:g} "
"electrons; this level is drawn with "
f"{level.electrons})")
else:
level.partial = True
note = ("NOTE: aufbau leaves this level PARTIALLY "
f"filled ({level.electrons} of {capacity} "
"electrons) at this k point -- a genuinely "
"metallic k point, or a level-ordering "
"artifact of the model")
level.detail = (note if not level.detail
else f"{level.detail}\n{note}")
remaining -= take
if remaining <= 0:
break
return partial
# ---------------------------------------------------- wave-function sketch
[docs]
def supercell_for(self, kpoint):
"""Display supercell of the hover sketch at ``k``.
Default: the k-commensurate diagonal supercell of the primitive
cell. With ``conventional`` the display cell is the conventional
cell of the detected centring (times the diagonal multiples that
make ``exp(2 pi i k . T) = 1`` for its edge vectors, as in the SALC
viewer), and each atom is wrapped into it with its own
primitive-lattice translation -- the Bloch phases stay exact because
every conventional-cell position is a primitive-lattice translate
of a basis atom.
Args:
kpoint: Three primitive reciprocal coordinates.
Returns:
``(sites, symbols, cartesian, display_lattice, description)``:
one ``(primitive atom index, integer lattice translation)`` pair
per drawn atom, the chemical symbols, the Cartesian positions in
Angstrom (shape ``(n, 3)``), the display lattice vectors as
rows, and a text such as ``"primitive cell, 2 x 2 x 2"``.
"""
from fractions import Fraction
from itertools import product
if self.conventional:
from .phonon_vector import (get_commensurate_supercell_matrix,
get_conventional_matrix)
centring = self.builder.spglib_dataset["international"][0]
base_matrix = get_conventional_matrix(centring)
cell_matrix = np.array(
get_commensurate_supercell_matrix(kpoint, base_matrix),
dtype=int)
multiples = np.rint(np.diag(
cell_matrix @ np.linalg.inv(np.array(base_matrix, dtype=float))
)).astype(int)
description = (f"conventional cell ({centring} centring), "
f"{multiples[0]} x {multiples[1]} x {multiples[2]}")
else:
repetitions = [
Fraction(float(value)).limit_denominator(12).denominator
for value in kpoint
]
cell_matrix = np.diag(repetitions).astype(int)
description = (f"primitive cell, {repetitions[0]} x "
f"{repetitions[1]} x {repetitions[2]}")
display_lattice = np.array(cell_matrix, dtype=float) @ self.lattice
inverse_cell = np.linalg.inv(np.array(cell_matrix, dtype=float))
corners = np.array([np.asarray(shift) @ cell_matrix
for shift in product((0, 1), repeat=3)])
t_low = corners.min(axis=0) - 1
t_high = corners.max(axis=0) + 1
eps = 1e-6
sites, symbols, cartesian = [], [], []
for t1 in range(int(t_low[0]), int(t_high[0]) + 1):
for t2 in range(int(t_low[1]), int(t_high[1]) + 1):
for t3 in range(int(t_low[2]), int(t_high[2]) + 1):
translation = np.array([t1, t2, t3], dtype=float)
for index, symbol in enumerate(self.symbols):
frac = ((self.positions[index] + translation)
[docs]
@ inverse_cell)
if np.any(frac < -eps) or np.any(frac >= 1 - eps):
continue
sites.append((index, translation))
symbols.append(symbol)
cartesian.append((self.positions[index] + translation)
@ self.lattice)
return (sites, symbols, np.array(cartesian), display_lattice,
description)
def sketch_partners(self, level: DiagramLevel, kpoint, sites):
"""Real wave-function amplitudes of a level on the supercell atoms.
The hover sketch of one level, one entry per degenerate partner.
Only the ``sketch_specs`` components are drawn (all shells unless
``sketch_tokens`` restricted them).
The amplitudes are Re[psi] (or Im[psi] when Re vanishes) of the
Bloch crystal orbital, so the sign alternation between the cells of
the k-commensurate supercell is displayed faithfully; degenerate
partners are realified and RREF-canonicalized like the molecular
sketch.
Same-l shells of one atom (e.g. Sc 2p/3p/4p) share their slots and
accumulate, each weighted by its STO radial amplitude at a probe
radius, so the drawn lobe signs are the signs of the real wave
function there. (A bare coefficient of one shell is wrong: a
semicore level like Sc 3p would be drawn from the tiny
orthogonalization tail of the 4p shell, whose sign is inverted --
the crystal analogue of the contracted-GTO compression in the
molecular PySCF sketch.)
In the full-basis --diagram mode the lobe SIZE of each (atom, l)
channel is additionally calibrated to its Loewdin population
(level.channel_pop, attached by solve_at), and the probe radius is
chosen per channel as the one (1.5/2.0/2.5/3.0 bohr) where the
accumulated amplitude is largest -- the visualize_eht/PySCF-viewer
recipe. Raw amplitude x coefficient lobes misstate the mix badly:
the 89.9%-Ti-4s GM1+ level of rutile TiO2 was DRAWN as d lobes
(drawn d:s = 2.7:1), because the compact 3d weighs ~4x the diffuse
4s at a fixed 2-bohr radius while the semicore 3s orthogonality
tail cancels most of the s channel. A level without channel_pop
(or a sketch_specs-filtered sketch, an API-only mode -- the CLI
rejects --atomic-orbital with --diagram) keeps the legacy
fixed-radius raw amplitudes.
Args:
level: A ``DiagramLevel`` from :meth:`solve_at` (its ``vectors``
are used).
kpoint: The k point the level was solved at.
sites: The ``(atom index, translation)`` list of
:meth:`supercell_for` at that k point.
Returns:
One list per partner; each holds one sketch entry
``[atom, s, px, py, pz, dxy, dyz, dz2, dxz, dx2-y2]`` per
supercell atom (index into ``sites``, then the real amplitudes
of the nine ``s``/``p``/``d`` components).
"""
from .point_charge_field import _primitives
from .visualize_basis import realify_basis_space
rows, _ = realify_basis_space(level.vectors.T)
rows = np.asarray(rows)
width = 9
slot_of = {0: 0, 1: 1, 2: 4}
radii = (1.5, 2.0, 2.5, 3.0)
angular = {0: 0.28209479, 1: 0.48860251, 2: 0.63078313}
profiles = {}
for spec_index, spec in enumerate(self.specs):
if spec.l in slot_of:
profiles[spec_index] = tuple(
angular[spec.l] * sum(
c * radius ** (n - 1) * np.exp(-z * radius)
for c, n, z in _primitives(spec.n, spec.zeta)
)
for radius in radii
)
pops = (getattr(level, "channel_pop", None)
if self.sketch_specs is None else None)
# no populations available (a level built outside solve_at, or the
# filtered sketch): keep the legacy fixed-radius raw amplitudes
# entirely -- a best-radius pick without the Loewdin sizes would be
# an inconsistent hybrid
calibrate = pops is not None
fixed_radius = radii.index(2.0)
# pass 1: per-partner per-(site, l) channel amplitudes at each
# probe radius, in the primitive cell (no supercell phases yet)
prim = []
amp2: dict = {}
for vector in rows:
channels: dict = {}
for spec_index, spec in enumerate(self.specs):
if (self.sketch_specs is not None
and spec_index not in self.sketch_specs):
continue
if spec.l not in slot_of:
continue
w = 2 * spec.l + 1
for site_pos, site in enumerate(spec.sites):
start = spec.offset + site_pos * w
block = np.asarray(vector[start:start + w])
key = (site, spec.l)
entry = channels.setdefault(
key, [np.zeros(w, dtype=complex)
for _ in radii])
for radius_index in range(len(radii)):
entry[radius_index] = (
entry[radius_index]
+ block * profiles[spec_index][radius_index])
for key, per_radius in channels.items():
norms = np.array([float(np.linalg.norm(c)) ** 2
for c in per_radius])
amp2[key] = amp2.get(key, 0.0) + norms
prim.append(channels)
# per (atom, l): the probe radius with the largest multiplet
# amplitude (away from any semicore orthogonalization node), and
# the Loewdin size calibration
best = {key: (int(np.argmax(values)) if calibrate
else fixed_radius)
for key, values in amp2.items()}
scale = {}
for key, values in amp2.items():
reference = float(values[best[key]])
if reference < 1e-24:
continue
if pops is not None:
scale[key] = np.sqrt(max(pops.get(key, 0.0), 0.0)
/ reference)
else:
scale[key] = 1.0
# pass 2: place the calibrated channels on the supercell sites
# with their Bloch phases
partner_rows = []
for channels in prim:
amp_re = np.zeros((len(sites), width))
amp_im = np.zeros_like(amp_re)
for (site, l_channel), per_radius in channels.items():
key = (site, l_channel)
if key not in scale:
continue
base = per_radius[best[key]] * scale[key]
slot = slot_of[l_channel]
w = 2 * l_channel + 1
for row_index, (atom, translation) in enumerate(sites):
if atom != site:
continue
phase = np.exp(2j * np.pi * float(np.dot(
kpoint, translation + self.positions[site]
)))
values = base * phase
amp_re[row_index, slot:slot + w] += values.real
amp_im[row_index, slot:slot + w] += values.imag
choice = (amp_re if np.linalg.norm(amp_re) >= np.linalg.norm(amp_im)
else amp_im)
partner_rows.append(choice.reshape(-1))
rows_arr = np.array(partner_rows)
if len(partner_rows) > 1:
from .molecular_salc import _rref_orthogonal
work = rows_arr.copy()
peak = np.max(np.abs(work)) or 1.0
work[:, np.max(np.abs(work), axis=0) < 0.05 * peak] = 0.0
canonical = _rref_orthogonal(list(work))
if len(canonical) == len(partner_rows):
rows_arr = np.array(canonical, dtype=float)
partners = []
for row in rows_arr:
grid = row.reshape(len(sites), width)
peak = np.max(np.abs(grid)) or 1.0
entries = []
for atom_index in range(grid.shape[0]):
values = grid[atom_index] / peak
if np.max(np.abs(values)) >= 0.04:
entries.append(
[atom_index] + [round(float(x), 3) for x in values]
)
partners.append(entries)
return [p for p in partners if p] or None
# ------------------------------------------------------------- k points
[docs]
def special_kpoints(self):
"""Tabulated special k points of the space group.
Returns:
``[(name, kpoint), ...]`` with the ISO-IR names (``GM``, ``R``,
``X``, ``M`` for Pm-3m) and primitive reciprocal coordinates, in
table order: the k points ``crystod --diagram`` draws
(``--kpoint NAME`` restricts the run to one of them).
"""
table = IrrepTable(self.builder.spglib_dataset["number"], spinor=False)
primitive_matrix = get_primitive_matrix_by_centring(
self.builder.spglib_dataset["international"][0]
)
names, kpoints = [], []
for irrep in table.irreps:
kpoint = snap_qpoint(np.array(irrep.k) @ primitive_matrix)
if kpoint not in kpoints:
kpoints.append(kpoint)
names.append(irrep.kpname)
return list(zip(names, kpoints))
# --------------------------------------------------------------------- output
def _format_kpoint(kpoint) -> str:
from fractions import Fraction
parts = []
for value in kpoint:
fraction = Fraction(float(value)).limit_denominator(12)
parts.append(str(fraction))
return "(" + ",".join(parts) + ")"
def _detail_html(level: DiagramLevel, names: dict) -> str:
# consumed as the SVG <title> textContent (the native hover tooltip),
# which renders newlines but shows HTML tags literally
# atomic-shell (ao) levels carry no irrep and electrons=None
rows = [f"E = {level.energy:.2f} eV"]
if level.irrep:
rows.append(f"irrep: {level.irrep} (degeneracy {level.degeneracy})")
else:
rows.append(f"degeneracy {level.degeneracy}")
if level.electrons is not None:
rows.append(f"electrons: {level.electrons}")
if level.detail:
rows.append(level.detail)
return "\n".join(rows)
def _periodic_sketch_geometry(display_lattice, symbols, positions,
description=None):
"""VESTA-style periodic geometry for the orbital sketch.
Atoms sitting on the display-cell boundary are drawn at every
translationally equivalent boundary position (an atom at fractional 0
also appears at 1, a corner atom at all eight corners), and the
display-cell outline is passed along as ``geometry["cell"]`` (origin +
three lattice vectors, in the centered coordinates the sketch uses)
for the dashed frame; ``description`` becomes ``geometry["desc"]``
(shown under the sketch). Returns (geometry, replica_map) with
replica_map[source_row] = [rows of its boundary images].
"""
from itertools import product
from .mo_diagram import diagram_geometry
super_lattice = np.asarray(display_lattice, dtype=float)
fractional = np.asarray(positions, dtype=float) @ np.linalg.inv(super_lattice)
tolerance = 1e-6
all_symbols = list(symbols)
all_positions = [np.asarray(p, dtype=float) for p in positions]
replica_map: dict[int, list[int]] = {}
for row, frac in enumerate(fractional):
choices = []
for axis in range(3):
choice = [0.0]
if frac[axis] < tolerance:
choice.append(1.0)
elif frac[axis] > 1.0 - tolerance:
choice.append(-1.0)
choices.append(choice)
for shift in product(*choices):
if not any(shift):
continue
replica_map.setdefault(row, []).append(len(all_symbols))
all_symbols.append(symbols[row])
all_positions.append(
all_positions[row] + np.asarray(shift) @ super_lattice)
geometry = diagram_geometry(all_symbols, np.array(all_positions))
# frame in the same centered frame as the atoms (diagram_geometry
# subtracts the centroid of the positions it is given)
center = np.array(all_positions).mean(axis=0)
geometry["cell"] = (
[[round(float(x), 4) for x in -center]]
+ [[round(float(x), 4) for x in row] for row in super_lattice]
)
corners = np.array([
np.asarray(shift) @ super_lattice - center
for shift in product((0.0, 1.0), repeat=3)
])
geometry["radius"] = max(
float(geometry["radius"]),
float(np.max(np.linalg.norm(corners, axis=1))),
)
if description:
geometry["desc"] = description
return geometry, replica_map
def _with_periodic_images(partners, replica_map):
"""Copy each sketch entry onto the boundary images of its atom.
The images carry IDENTICAL amplitudes: within the k-commensurate
supercell the Bloch wave function is exactly periodic
(exp(i k . T_super) = 1 by construction), which is the whole point of
drawing that supercell.
"""
if not partners or not replica_map:
return partners
extended = []
for entries in partners:
extra = [[image] + entry[1:]
for entry in entries
for image in replica_map.get(entry[0], [])]
extended.append(entries + extra)
return extended
def write_crystal_diagram_html(diagram: CrystalOrbitalDiagram,
k_entries: list, output_path: str,
structure_label: str) -> None:
"""One interactive page with one energy diagram per k point."""
# both engines add the outermost atomic-shell columns (MolOD's
# "ligand-ao" analogue): the EHT engine as VSIP + ligand-field on-site
# levels, the PySCF engine as isolated formal-charge-ion calculations;
# only --onsite pages keep the three-column layout
has_ao = bool(k_entries) and "left-ao" in k_entries[0][2]
if has_ao:
# left-ao at 155 keeps its end-anchored labels ("2Al 3p" ends at
# x = 155 - 26 - 9 = 120) clear of the energy axis (line at 66,
# tick numbers up to 57), mirroring MolOD's ligand-ao spacing
columns = {"left-ao": 155, "left": 325, "mo": 495,
"right": 675, "right-ao": 865}
half = {"left-ao": 26, "left": 30, "mo": 34,
"right": 30, "right-ao": 26}
order = ["left-ao", "left", "mo", "right", "right-ao"]
side = {"left-ao": -1, "left": -1, "mo": 1,
"right": 1, "right-ao": 1}
else:
columns = {"left": 200, "mo": 480, "right": 760}
half = {"left": 34, "mo": 34, "right": 34}
order = ["left", "mo", "right"]
side = {"left": -1, "mo": 1, "right": 1}
def ao_header(column):
counts: dict[str, set] = {}
for spec in diagram.side_specs[column]:
counts.setdefault(spec.element, set()).update(spec.sites)
return " + ".join(
(f"{len(sites)}{el}" if len(sites) > 1 else el)
for el, sites in counts.items()
) + " AOs"
headers = {
"left": svg_sub_digits(diagram.formula["left"]),
"mo": "crystal orbitals",
"right": svg_sub_digits(diagram.formula["right"]),
}
if has_ao:
headers["left-ao"] = ao_header("left")
headers["right-ao"] = ao_header("right")
variants = []
for name, kpoint, levels in k_entries:
(sites, super_symbols, super_positions, display_lattice,
cell_description) = diagram.supercell_for(kpoint)
geometry, replica_map = _periodic_sketch_geometry(
display_lattice, super_symbols, super_positions,
cell_description)
names = {
level.level_id: level.label
for column in ("left", "right") for level in levels[column]
}
levels_json = []
bond_letter = {"bonding": "b", "nonbonding": "n", "antibonding": "a"}
from .mo_diagram import element_color
for column in order:
for level in levels.get(column, []):
character = getattr(level, "bond_character", None)
# sublattice levels carry the VESTA color of their dominant
# element ("Sc 3d GM5+" -> the Sc color); atomic-shell
# levels ("2Al 3s") strip the leading site count
elc = None
if column.endswith("-ao"):
elc = element_color(
level.label.split()[0].lstrip("0123456789"))
elif column != "mo":
parts = level.label.split()
if len(parts) >= 3:
elc = element_color(parts[0])
levels_json.append({
**({"elc": elc} if elc else {}),
"id": level.level_id,
"col": level.column,
"e": round(level.energy, 4),
"deg": level.degeneracy,
# near-dependent Bloch combination: first-order Loewdin
# energy estimate (machine-readable marker; the human-
# facing note lives in the tooltip detail)
**({"est": 1} if getattr(level, "estimated", False)
else {}),
**({"bond": bond_letter[character]} if character else {}),
"label": level.label,
# atomic-shell levels carry electrons=None (no arrows;
# occupation is a sublattice/crystal-column concept)
"el": level.electrons,
"occ": bool(level.electrons),
"links": [
[i, round(w, 4)]
for i, w in level.composition if w >= 0.02
],
# crystal levels: the panel bars quote the SAME AO-shell
# populations as the hover tooltip (display_composition);
# the fragment projection only draws the links above
"comp": (
[[label, round(100 * w, 1)]
for label, w in level.display_composition]
if getattr(level, "display_composition", None)
is not None else
[[names[i], round(100 * w, 1)]
for i, w in sorted(level.composition,
key=lambda kv: -kv[1])
if w > 0.005]
),
"detail": _detail_html(level, names),
# None (not []) for the atomic-shell levels: an empty
# array is truthy in JS and would open a lobe-less
# sketch pane
"orb": (None if level.vectors is None else
_with_periodic_images(
diagram.sketch_partners(level, kpoint, sites),
replica_map)),
})
occupied = [lv for lv in levels["mo"] if lv.electrons > 0]
empty = [lv for lv in levels["mo"] if lv.electrons == 0]
homo_level = max(occupied, key=lambda lv: lv.energy) if occupied else None
lumo_level = min(empty, key=lambda lv: lv.energy) if empty else None
homo = homo_level.level_id if homo_level else None
lumo = lumo_level.level_id if lumo_level else None
energies = [lv.energy for column in order
for lv in levels.get(column, [])]
e_min, e_max = min(energies), max(energies)
padding = 0.08 * (e_max - e_min) or 1.0
# the interactive view opens on the frontier states: +-8 eV around the
# HOMO/LUMO midpoint ("Show all energy levels" reveals the deep shells
# outside it); without a HOMO/LUMO pair, fall back to a fixed window
if homo_level is not None and lumo_level is not None:
center = 0.5 * (homo_level.energy + lumo_level.energy)
view_lo, view_hi = center - _VIEW_HALF_WINDOW, center + _VIEW_HALF_WINDOW
else:
view_lo = max(e_min - padding, _VIEW_E_MIN)
view_hi = min(e_max + padding, _VIEW_E_MAX)
if view_hi - view_lo < 1.0:
view_lo, view_hi = e_min - padding, e_max + padding
variants.append({
"key": f"{name} {_format_kpoint(kpoint)}",
"levels": levels_json,
"homo": homo,
"lumo": lumo,
"eMin": round(view_lo, 2),
"eMax": round(view_hi, 2),
"geom": geometry,
})
first = variants[0]
fragment_names = (f"{svg_sub_digits(diagram.formula['left'])} + "
f"{svg_sub_digits(diagram.formula['right'])}")
chips = [
structure_label,
f"{diagram.builder.spglib_dataset['international']} "
f"(No. {diagram.builder.spglib_dataset['number']})",
f"{fragment_names} "
+ getattr(diagram, "basis_chip", "(full-electron basis)"),
f"{int(diagram.electrons)} electrons / cell",
# --onsite pages replace the point-charge chip: there is no
# point-charge embedding in the single-Hamiltonian mode
getattr(diagram, "embedding_chip", "") or (
"point charges: " + " ".join(
f"{element}{diagram.oxidation[element]:+g}"
for element in dict.fromkeys(diagram.symbols)
)
),
getattr(diagram, "method_chip", "extended Hückel + SALC"),
]
sketch_cell = (
"the conventional cell (× the multiples that make the Bloch "
"phase commensurate)" if getattr(diagram, "conventional", False)
else "the k-commensurate supercell"
)
sketch_foot = (
" Click or hover a level for its composition and the real-space "
"wave function Re[ψ] of all its atomic-orbital components drawn "
f"on {sketch_cell} (drag to rotate; degenerate "
"partners switchable)."
)
ao_foot = ""
if has_ao:
# engines override ao_foot (the PySCF columns are isolated-ion
# calculations, not VSIP levels)
ao_foot = getattr(diagram, "ao_foot", "") or (
" The outermost columns are each element's isolated atomic "
"shells at their on-site energies (VSIP + spherical part of "
"the point-charge ligand field, k-independent); their "
"connector lines into the sublattice columns show how the "
"atoms' Bloch combinations split each shell at the "
"chosen k point (with several atoms per cell the "
"intra-sublattice overlap can reorder shells, e.g. 3p "
"combinations below 3s)."
)
estimated_foot = ""
if any(getattr(level, "estimated", False)
for _, _, levels in k_entries
for column in ("left", "mo", "right") for level in levels[column]):
estimated_foot = (
" Some levels are near-dependent diffuse Bloch combinations "
f"(overlap eigenvalue below {_OVERLAP_FLOOR}) whose variational "
"extended-Hückel energy diverges (overlap catastrophe): "
"their energies are first-order Löwdin estimates, marked "
"~ in the terminal report and noted in the level's tooltip."
)
bond_foot = ""
if any(getattr(level, "bond_character", None)
for _, _, levels in k_entries for level in levels["mo"]):
bond_foot = (
" Crystal-orbital line colors: "
"<span style=\"color:#1565c0\">bonding</span> / "
"<span style=\"color:#333\">nonbonding</span> / "
"<span style=\"color:#d32f2f\">antibonding</span>, from the "
"left–right overlap population 2 Re c<sub>L</sub>†"
"S c<sub>R</sub> of each state (COOP-style; semicore "
"orthogonality tails excluded; value in the level's tooltip)."
)
render_diagram_page(
output_path,
title=f"Crystal orbital diagram: {structure_label}",
heading_html=(
f"Crystal-orbital diagram: {structure_label} "
f"<span style=\"color:#90a4ae;font-size:14px\">{fragment_names}, "
"per k point</span>"
),
chips=chips,
columns=columns, half=half, order=order, side=side, headers=headers,
levels_json=first["levels"],
homo_id=first["homo"],
lumo_id=first["lumo"],
e_min=first["eMin"], e_max=first["eMax"],
foot_html=(
# engines override foot_intro (the PySCF pages used to reuse the
# extended-Hueckel energetics sentence below, which was wrong)
(getattr(diagram, "foot_intro", "") or (
"Crystal-orbital diagram (COD): the fragment-sublattice "
"Bloch orbitals (columns, full-electron basis: every core "
"and valence shell) are the electronic states before "
"chemical bond formation, in the point-charge ligand field "
"of the removed sublattice (formal oxidation states; exact "
"multipole + penetration terms, background-dependent "
"monopole omitted); states sharing an irrep of the little "
"group at k mix into bonding/antibonding crystal orbitals "
"(center), states without a partner remain nonbonding. "
"Energies: symmetry-adapted extended Hückel "
"(VSIP/core-level diagonal + Wolfsberg-Helmholz "
"off-diagonal over exact Bloch STO overlap sums). "
"<b>CAUTION</b>: this is a non-self-consistent one-electron "
"model with fixed atomic parameters. What is rigorous is the "
"SYMMETRY — the irrep labels, which states may mix and "
"which stay nonbonding. The level ORDER can be qualitatively "
"wrong (no charge-transfer response; a diffuse cation shell "
"in the point-charge field can drop a whole band into the "
"wrong place, as the Ti 4p band does inside the O 2p band of "
"rutile TiO2, breaking its insulating filling), and a wrong "
"order inverts the compositions it feeds (in that same TiO2 "
"the bonding GM3+ comes out Ti-dominated and its antibonding "
"partner O-dominated — the reverse of the heteropolar "
"rule and of PySCF). Cross-check with the same command plus "
"--pyscf wherever it runs; it cannot run for the lanthanides "
"(no GTH basis covers them), where this engine is the only "
"one and the diagram should be read as a symmetry analysis."))
+ " The energy window opens on the frontier states; use \"Show "
"all energy levels\" for the deep shells. Switch the k point "
"with the buttons above." + ao_foot + estimated_foot + bond_foot
+ sketch_foot
),
geometry=variants[0]["geom"],
variants=variants,
)
def report_and_write(cell, *, left, right, symprec, electrons,
kpoint_filter, output_path, structure_label,
oxidation=None, conventional=False):
"""Terminal report + HTML for the crystal-orbital diagram."""
diagram = CrystalOrbitalDiagram(
cell, left, right, symprec=symprec, electrons=electrons,
oxidation=oxidation, conventional=conventional,
)
dataset = diagram.builder.spglib_dataset
print("\n * Space group *")
print(f" {dataset['international']} ({dataset['number']})\n")
print(" * Fragments (full-electron basis: core + valence shells;"
" atomic levels from reference/atomic_level_*) *")
for column in ("left", "right"):
by_element: dict[str, list] = {}
for spec in diagram.side_specs[column]:
by_element.setdefault(spec.element, []).append(spec)
parts = " | ".join(
f"{element} " + " ".join(spec.shell for spec in specs)
+ f" x{len(specs[0].sites)} site(s)"
+ (f" [ECP-{ATOMIC_LEVELS[element]['ecp_core']} core frozen]"
if ATOMIC_LEVELS[element]["ecp_core"] else "")
for element, specs in by_element.items()
)
print(f" {column:<5} {diagram.formula[column]:<6}: {parts}, "
f"{diagram.side_electrons[column]} electrons")
# the outermost diagram columns: isolated atomic shells at the
# (k-independent) on-site energies; a +-x marker discloses the
# per-Wyckoff-site spread of multi-position elements
onsites = []
seen: set = set()
cage_shells = []
for spec in diagram.side_specs[column]:
key = (spec.element, spec.shell)
if key not in seen:
seen.add(key)
block = diagram.h_bar[spec.offset:spec.offset + spec.n_ao]
spread = float(np.max(block) - np.min(block))
onsite = float(np.mean(block))
raw_vsip = float(np.mean(
diagram.h_raw[spec.offset:spec.offset + spec.n_ao]))
field_shift = onsite - raw_vsip
onsites.append(
f"{spec.element} {spec.shell} {onsite:.2f}"
+ (f"(+-{spread / 2:.2f})" if spread > 0.02 else "")
+ (f" [VSIP {raw_vsip:.2f}, field {field_shift:+.2f}]"
if abs(field_shift) >= 2.0 else ""))
if abs(field_shift) > 8.0:
cage_shells.append(
(spec.element, spec.shell, field_shift))
print(" on-site atomic levels (VSIP + ligand field, eV): "
+ ", ".join(onsites))
for element, shell, field_shift in cage_shells:
print(f" WARNING: the {element} {shell} point-charge "
f"shift ({field_shift:+.1f} eV) is no small correction: "
"this STO is diffuse enough to engulf the surrounding "
"charge cage, so its on-site energy -- and every level "
"it dominates -- is an artifact of the point-charge "
"model, not chemistry")
print(f" electrons per cell in the diagram: {int(diagram.electrons)}"
+ (" (all electrons of the neutral atoms; override with"
" --electrons)" if electrons is None else ""))
print(" * Ligand-field point charges (removed sublattice) *")
shifts: dict[str, list] = {}
for atom, symbol in enumerate(diagram.symbols):
shifts.setdefault(symbol, []).append(diagram.site_potential[atom])
other = {"left": "right", "right": "left"}
for column in ("left", "right"):
felt = " + ".join(
f"{element}^{diagram.oxidation[element]:+g}"
for element in dict.fromkeys(
spec.element for spec in diagram.side_specs[other[column]]
)
)
own = ", ".join(
f"{element} {np.mean(shifts[element]):+.2f} eV"
for element in dict.fromkeys(
spec.element for spec in diagram.side_specs[column]
)
)
print(f" {column:<5} {diagram.formula[column]:<6} feels the {felt} "
f"lattice (multipole ligand field + penetration; "
f"jellium-referenced monopole {own} omitted)")
print(" hover wave-function sketches: all atomic-orbital components "
"of every level")
print(" * CAUTION: extended Hueckel gets the LEVEL ORDER wrong "
"sometimes *")
print(" This is a non-self-consistent one-electron model with fixed "
"atomic parameters.\n"
" What is rigorous here is the SYMMETRY: the irrep labels, "
"which states may mix,\n"
" and which stay nonbonding. The energies carry no "
"charge-transfer response, so\n"
" a whole band can land in the wrong place -- rutile TiO2 puts "
"an artificial Ti 4p\n"
" band inside the O 2p band and loses its insulating filling. "
"A wrong ORDER also\n"
" inverts the COMPOSITIONS it feeds: in the same TiO2 the "
"bonding GM3+ comes out\n"
" Ti-dominated and the antibonding one O-dominated, the reverse "
"of the\n"
" heteropolar rule and of PySCF. Cross-check with --pyscf "
"(same command +\n"
" --pyscf) wherever it runs -- it cannot for the lanthanides, "
"which no GTH basis\n"
" covers; there this engine is the only one, read it as a "
"symmetry analysis.\n")
entries = []
kpoints = diagram.special_kpoints()
if kpoint_filter is not None:
available = [name for name, _ in kpoints]
kpoints = [
(name, kpoint) for name, kpoint in kpoints
if name == kpoint_filter
]
if not kpoints:
raise SystemExit(
f"ERROR: k point '{kpoint_filter}' is not a special point of "
f"this space group (available: {', '.join(available)})."
)
for name, kpoint in kpoints:
levels, _ = diagram.solve_at(kpoint)
entries.append((name, kpoint, levels))
print(f" * k point {name} {_format_kpoint(kpoint)} *")
if diagram.last_estimated:
print(f" ({diagram.last_estimated} near-dependent diffuse Bloch "
f"combination(s) below overlap floor {_OVERLAP_FLOOR}: "
"energies marked ~ are first-order Loewdin estimates; the "
"variational extended-Hueckel values diverge)")
if diagram.last_dependent:
print(f" ({diagram.last_dependent} linearly dependent Bloch "
"combination(s) removed by canonical orthogonalization)")
for col_name, col_levels in (("crystal", levels["mo"]),
(diagram.formula["left"],
levels["left"]),
(diagram.formula["right"],
levels["right"])):
for lv in col_levels:
if getattr(lv, "partial", False):
marker = "~" if getattr(lv, "estimated", False) else ""
print(f" WARNING: aufbau leaves the {col_name} "
f"column's {lv.irrep} level at "
f"{marker}{lv.energy:.2f} eV "
f"partially filled ({lv.electrons} of "
f"{2 * lv.degeneracy} electrons) -- a genuinely "
"metallic k point, or a level-ordering artifact "
"(an occupied level dominated by a strongly "
"shifted diffuse shell is the usual culprit; "
"see the on-site warnings above)")
for column in ("left", "right"):
parts = ", ".join(
f"{lv.label} ({'~' if lv.estimated else ''}{lv.energy:.2f})"
for lv in sorted(levels[column], key=lambda lv: lv.energy)
)
print(f" {diagram.formula[column]:<10}: {parts}")
print(" crystal :")
for lv in sorted(levels["mo"], key=lambda lv: lv.energy):
# the same AO-population list as the HTML panel and tooltip
composition = " ".join(
f"{label} {100 * w:.1f}%"
for label, w in getattr(lv, "display_composition", [])
)
occupancy = f"{lv.electrons}e" if lv.electrons else " "
energy_str = f"{'~' if lv.estimated else ''}{lv.energy:.2f}"
print(f" {lv.label:<10} {energy_str:>9} eV x{lv.degeneracy}"
f" {occupancy:<4} {composition}")
print("")
write_crystal_diagram_html(diagram, entries, output_path, structure_label)
print(f"Crystal-orbital diagram written to {output_path}")
def main(argv: list[str] | None = None) -> None:
parser = argparse.ArgumentParser(
description="Crystal-orbital diagram from symmetry + extended-Hueckel "
"overlap."
)
parser.add_argument("--poscar", default="POSCAR")
parser.add_argument("--co-left", nargs="+", required=True,
metavar="FORMULA",
help="left fragment sublattice, e.g. SrTi")
parser.add_argument("--co-right", nargs="+", required=True,
metavar="FORMULA",
help="right fragment sublattice, e.g. O3")
parser.add_argument("--oxidation", nargs="+", default=None,
metavar="EL=Q",
help="formal oxidation states for the removed-"
"sublattice point charges, e.g. Sr=+2 Ti=+4 O=-2 "
"(default: guessed with pymatgen)")
parser.add_argument("--kpoint", default=None,
help="restrict to one special k point label (e.g. GM)")
parser.add_argument("--electrons", type=float, default=None,
help="electrons per primitive cell "
"(default: neutral-atom valence counts)")
parser.add_argument("--conventional", action="store_true",
help="draw the hover wave-function sketches in the "
"conventional cell instead of the primitive "
"k-commensurate supercell")
parser.add_argument("--output", default=None)
parser.add_argument("--tolerance", type=float, default=1e-5)
args = parser.parse_args(argv)
from .star_of_k import read_poscar_or_exit
cell = read_poscar_or_exit(args.poscar)
from pathlib import Path
stem = Path(args.poscar).name
for extension in (".vasp", ".poscar"):
if stem.lower().endswith(extension):
stem = stem[: -len(extension)]
output_path = args.output or f"CrystOD_{stem}.html"
report_and_write(
cell,
left=args.co_left,
right=args.co_right,
symprec=args.tolerance,
electrons=args.electrons,
kpoint_filter=args.kpoint,
output_path=output_path,
structure_label=stem,
oxidation=(parse_oxidation_tokens(args.oxidation)
if args.oxidation else None),
conventional=args.conventional,
)
if __name__ == "__main__":
main()
