crystod.salc

crystod.salc#

Crystal-orbital SALC analysis: the Python face of the crystod command.

Everything the main command does with a structure file is available here with a phonopy.structure.atoms.PhonopyAtoms cell in place of the -c POSCAR argument: the irreducible representations of the crystal orbitals – the symmetry-adapted linear combinations (SALCs) of one atomic shell’s Bloch sums – at any k point, the SALC coefficient bases behind the interactive 3D viewer, the crystal-orbital diagrams built from extended-Hueckel or periodic PySCF overlaps, and the star of a k point. The classes take the same inputs as the command-line flags they mirror and return the numbers the command prints.

SALC irreps at a k point (crystod -c POSCAR --element EL --orbital ORB)

CrystalOrbital: Irrep decomposition of one element’s shell at a k point, or at every special point, with ISO-IR labels; --spinor is the spior flag of the constructor.

Crystal-orbital diagrams (crystod --diagram)

CrystalOrbitalDiagram: The symmetry + extended-Hueckel engine: fragment sublattices from --co-left/--co-right, full core + valence basis, one solve_at per special k point.
PySCFCrystalOrbitalDiagram: The quantitative engine of --pyscf: three periodic PySCF calculations sharing one AO space, deep-level column alignment (needs pip install "CrystOD[quantum]").
assign_bond_characters(): The COOP bonding/antibonding/nonbonding classification both engines apply to their levels.

Symmetry-adapted orbital bases (crystod --visualize)

SymmetryAdaptedOrbitalBasis: Explicit SALC coefficient vectors of one shell at a k point, from the projected representation matrices; the data of the SALC viewer.

Star of k (crystod --star-of-k)

compute_star(): The arms of the star of a k point and the operations reaching each.
format_star_lines(): The report lines of those arms.
resolve_kpoint_input(): A --kpoint argument (label or coordinates) as label plus coordinates.

Attributes resolve lazily (PEP 562): import crystod.salc is instant; the implementation modules, and with them phonopy and spgrep (and pyscf for the PySCF diagram), are imported the first time an attribute is used. Bad input that the command line reports as ERROR: ... and exits on is raised as ValueError from the functions of this namespace; the classes raise SystemExit from their constructors and methods, as the implementation modules do.

class crystod.salc.CrystalOrbital(cell, symprec=1e-05, spior=False)[source]#

Bases: object

Irreps of the crystal orbitals built from one atomic shell (crystod).

The Bloch sums of one element’s s, p, d, f, g, h or i shell at a k point span a (reducible) representation of the little group of k – the site-symmetry induced representation, or band representation, of that shell. Its character is the product of the permutation character of the element’s sites (with the Bloch phases of k) and the rotation character of the shell, and its decomposition into the irreps of the little group is what crystod -c POSCAR --element EL --orbital ORB [--kpoint K] prints: with --spinor the double-valued (spin-orbit) irreps, without --kpoint every special point of the space group. The spgrep irreps are labelled with the ISO-IR (Miller-Love) names, GM3+(2), X5-(2) and so on, the number in parentheses being the dimension.

The input cell is reduced to the spglib primitive cell, and the symmetry operations are reordered to match the ISO-IR tables (a ValueError is raised when that fails).

Parameters:

cell (PhonopyAtoms) – The crystal structure as phonopy.structure.atoms.PhonopyAtoms (any setting; it is converted to the primitive cell).
symprec (float) – Symmetry tolerance handed to spglib.
spior (bool) – True for the double-valued (spinor) irreps. The parameter is spelled this way in the signature; the attribute is spinor.

Variables:

primitive_cell – The standardized primitive PhonopyAtoms cell that every k point and atom index refers to.
spglib_dataset – The spglib symmetry dataset of the primitive cell ("international", "number", "wyckoffs", "site_symmetry_symbols", …).
transformation_matrix – The spglib transformation matrix of the primitive cell; it carries the ISO-IR k vectors and rotations, tabulated in the conventional setting, onto the primitive basis.
rotations – Integer rotation matrices in the primitive basis, in the ISO-IR table order.
translations – The matching fractional translations.
seitz_symbols – The Seitz symbol of every operation, same order.
spinor – Whether double-valued irreps are used.
symprec – The symmetry tolerance in use.
irt_character_table – The ISO-IR table of the space group (double-valued when spinor).
irt_kpoint_table – The single-valued ISO-IR table the special k points are enumerated from.
labels_from_isoir – Set by get_irrep_labels() at a k point the table does not list: whether the labels came from the full ISO-IR k-vector data.

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")
>>> co = salc.CrystalOrbital(cell)
>>> co.get_irt_special_points()[0]
['GM', 'R', 'X', 'M']
>>> _, _, counts, labels = co.irreducible_decomposition(
...     [0, 0, 0], "Sc", "d")
>>> {labels[key]: int(n) for key, n in counts.items() if n > 0}
{'GM3+(2)': 1, 'GM5+(3)': 1}

calc_reducible_characters(k, element, orbital, mapping_little_group)[source]#

Characters of the crystal-orbital representation at k.

Parameters:

k (list[float]) – Three primitive reciprocal coordinates.
element (str) – Chemical symbol.
orbital (str) – Shell letter ("s" to "i").
mapping_little_group (ndarray[tuple[Any, ...], dtype[int64]]) – Indices into rotations of the little-group operations.

Returns:

the product of the permutation character of the element’s atoms (get_permutation_characters()) and the rotation character of the shell (get_atomic_orbital_characters()).

Return type:

Complex array, one entry per little-group operation

get_atomic_orbital_characters(rotations, orbital)[source]#

Rotation characters of one atomic shell.

The character of a proper rotation by alpha on the 2l+1 orbitals of shell l is sin((l + 1/2) alpha) / sin(alpha / 2) (2l+1 for the identity); an improper operation multiplies it by (-1)^l. With spinor the double-group factor 2 cos(alpha / 2) (2 for the identity) is included.

Parameters:

rotations (ndarray[tuple[Any, ...], dtype[int64]]) – Integer rotation matrices, shape (n, 3, 3).
orbital (str) – Shell letter "s", "p", "d", "f", "g", "h" or "i".

Returns:

Real array of the n characters.

Raises:

ValueError – An unknown shell letter.

Return type:

ndarray[tuple[Any, …], dtype[float64]]

get_irrep_labels(k, irreps, mapping_little_group)[source]#

Map spgrep irreps at k to ISO-IR labels by comparing characters.

Parameters:

k (list[float]) – Three primitive reciprocal coordinates.
irreps – The spgrep irreps at k (arrays of shape (little_group_order, dim, dim)), as returned by irreducible_decomposition().
mapping_little_group (ndarray[tuple[Any, ...], dtype[int64]]) – Indices into rotations of the little-group operations, in the order of the irrep matrices.

Returns:

{generic: label} with the generic key "irrep_i(dim)" of every irrep and its ISO-IR label such as "GM3+(2)". A k point the table does not list is first mapped onto the tabulated arm of its star (characters transported by conjugation), otherwise the full ISO-IR k-vector data are consulted; an irrep that still finds no match keeps its generic key as the label.

Return type:

dict[str, str]

get_irt_irreps_at_k(k)[source]#

Tabulated ISO-IR irreps at a k point.

Parameters:: k (list[float]) – Three primitive reciprocal coordinates.
Returns:: The ISO-IR irrep records whose k vector equals k – each with kpname, k and its characters – or an empty list when the point is not tabulated as given (a non-tabulated arm of a star, or no special point at all).
Return type:: list[IsoTableIrrep]

get_irt_special_points()[source]#

Special k points of the space group from the ISO-IR tables.

Returns:: the tabulated names (GM, R, X, M for Pm-3m) and their primitive reciprocal coordinates, one entry per distinct point, in table order. These are the points crystod analyzes when --kpoint is omitted.
Return type:: (names, kpoints)

get_kpoint_name(k)[source]#

Name of a k point (GM, X, M, …), or None.

Any arm of a tabulated star is recognized, not only the tabulated arm; a k point that is no special point receives the ISO-IR k-vector type letters (GP, DT, …) when those data are available.

Parameters:: k (list[float]) – Three primitive reciprocal coordinates.
Returns:: The name, or None when nothing matches.
Return type:: str | None

get_little_group(k)[source]#

Little group of k: the operations that leave k invariant.

Parameters:: k (list[float]) – Three primitive reciprocal coordinates.
Returns:: (little_rotations, little_translations) in the ISO-IR operation order (spgrep’s get_little_group without the index mapping, which irreducible_decomposition() returns).

get_little_group_symbol(k)[source]#

Space-group-type symbol of the little group of k.

Parameters:: k (list[float]) – Three primitive reciprocal coordinates.
Returns:: The short Hermann-Mauguin symbol and number spglib assigns to the little-group operations, e.g. "Pm-3m (221)" at Gamma or "P4/mmm (123)" at X of a cubic perovskite.
Return type:: str

get_modified_permutation_rep(r, t, k)[source]#

Bloch-phased permutation matrix of one operation at k.

Parameters:

r (ndarray[tuple[Any, ...], dtype[int64]]) – Integer rotation matrix in the primitive basis.
t (ndarray[tuple[Any, ...], dtype[float64]]) – Its fractional translation.
k (list[float, float, float]) – Three primitive reciprocal coordinates.

Returns:

Complex array of shape (n_atoms, n_atoms) whose entry [j, i] is the Bloch phase exp(2 pi i k . (R^-1 (x_j - t) - x_j)) when the operation sends atom i onto atom j (modulo lattice translations), and 0 elsewhere.

Return type:

ndarray[tuple[Any, …], dtype[complex128]]

get_permutation_characters(element, permutation_matrices)[source]#

Permutation characters restricted to the atoms of one element.

Parameters:

element (str) – Chemical symbol.
permutation_matrices (ndarray[tuple[Any, ...], dtype[complex128]]) – Output of get_permutation_reps_at_k().

Returns:

Complex array with the trace of every matrix over the element’s block of atoms.

Return type:

ndarray[tuple[Any, …], dtype[complex128]]

get_permutation_reps_at_k(little_rotations, little_translations, k)[source]#

Permutation matrices of the little-group operations at k.

Parameters:

little_rotations – Rotations of the little group of k, shape (order, 3, 3).
little_translations – Their fractional translations, shape (order, 3).
k (list[float]) – Three primitive reciprocal coordinates.

Returns:

one get_modified_permutation_rep() matrix per operation.

Return type:

Complex array of shape (order, n_atoms, n_atoms)

get_target_element_positions(element)[source]#

Indices of the atoms of element in the primitive cell.

Parameters:: element (str) – Chemical symbol, e.g. "Sc".
Returns:: The atom indices in cell order (a contiguous block in the standardized primitive cell).
Raises:: ValueError – The element is not in the cell.
Return type:: list[int]

irreducible_decomposition(k, element, orbital)[source]#

Decompose the crystal orbitals of one shell at k into irreps.

The calculation behind every line crystod prints: the irreps of the little group of k from spgrep (double-valued with spinor), the reducible characters from calc_reducible_characters(), and the multiplicities from the character inner product.

Parameters:

k (list[float, float, float]) – Three primitive reciprocal coordinates.
element (str) – Chemical symbol, e.g. "Sc".
orbital (str) – Shell letter "s" to "i", e.g. "d".

Returns:

mapping_little_group are the indices into rotations of the little-group operations, irreps the spgrep irrep matrices, counts maps the generic key "irrep_i(dim)" of every irrep to its multiplicity (rounded to two decimals), and labels maps the same keys to the ISO-IR labels (get_irrep_labels()).

Return type:

(mapping_little_group, irreps, counts, labels)

Raises:

ValueError – An element that is not in the cell, or an unknown shell letter.

class crystod.salc.CrystalOrbitalDiagram(cell, left_tokens, right_tokens, symprec=1e-05, electrons=None, sketch_tokens=None, oxidation=None, conventional=False)[source]#

Bases: object

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, report_and_write(), is special_kpoints(), then solve_at() per k point, then the HTML writer, which draws the hover wave-function sketches with supercell_for() and sketch_partners().

Parameters:

cell – The crystal structure as phonopy.structure.atoms.PhonopyAtoms (converted to the spglib primitive cell).
left_tokens (list[str]) – --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 (list[str]) – --co-right formula tokens, e.g. ["O3"].
symprec (float) – Symmetry tolerance handed to spglib.
electrons (float | None) – Electrons per primitive cell for the aufbau filling of the crystal column (default: all electrons of the neutral atoms).
sketch_tokens (list[str] | None) – El-shell tokens ("Ti-3d", "O_2p") that restrict the drawn sketch components; None draws every component (what the CLI does).
oxidation (dict[str, float] | None) – {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 (bool) – Draw the hover sketches in the conventional cell instead of the k-commensurate primitive supercell (display only).

Variables:

builder – The 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 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')

bloch_overlap(kpoint)[source]#

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.

Parameters:: kpoint – Three primitive reciprocal coordinates.
Returns:: Hermitian complex array of shape (n_ao, n_ao).
Return type:: ndarray

hamiltonian(S)[source]#

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.

Parameters:: S (ndarray) – The Bloch overlap matrix from bloch_overlap().
Returns:: Hermitian complex array of shape (n_ao, n_ao), in eV.
Return type:: ndarray

little_group_data(kpoint)[source]#

Irreps, labels and the AO representation of the little group at k.

Parameters:: kpoint – Three primitive reciprocal coordinates.
Returns:: 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.
Return type:: (irreps, mapping, labels, representation)

sketch_partners(level, kpoint, sites)[source]#

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.

Parameters:

level (DiagramLevel) – A DiagramLevel from solve_at() (its vectors are used).
kpoint – The k point the level was solved at.
sites – The (atom index, translation) list of 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).

solve_at(kpoint)[source]#

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 (assign_bond_characters()); the outermost columns list the isolated on-site shell levels. Sets last_estimated and last_dependent.

Parameters:: kpoint – Three primitive reciprocal coordinates, for the diagram one of special_kpoints().
Returns:: 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.
Return type:: (levels, labels)
Raises:: SystemExit – The representation does not leave S or H invariant (a gauge or real-harmonics inconsistency; please report the case).

special_kpoints()[source]#

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).

supercell_for(kpoint)[source]#

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.

Parameters:: kpoint – Three primitive reciprocal coordinates.
Returns:: 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".
Return type:: (sites, symbols, cartesian, display_lattice, description)

class crystod.salc.PySCFCrystalOrbitalDiagram(cell, left_tokens, right_tokens, *, symprec=1e-05, basis='gth-dzvp-molopt-sr', pseudo='gth-pbe', xc='pbe', kmesh=None, ke_cutoff=200.0, oxidation=None, electrons=None, sigma=0.0, degeneracy_tol=None, no_ghost=False, symmetrize=True, max_l=None, projection='lowdin', chk=None, onsite=False, conventional=False, conv_tol=1e-08, max_cycle=100, max_memory=4000.0, verbose=0)[source]#

Bases: CrystalOrbitalDiagram

Crystal-orbital diagram from three periodic PySCF calculations.

The engine behind crystod --diagram --pyscf. It subclasses the extended-Hueckel engine only to reuse its k-point list, degenerate-group clustering, irrep projection, level filling and supercell helper; the Hamiltonian, the overlap and the orbitals all come from PySCF (the module docstring describes the construction). The CLI sequence, report_and_write(), is run() (the SCFs, or a chk restart), prepare_bands() for all special_kpoints(), solve_at() and site_symmetry_irreps() per k point, align_fragment_columns(), atomic_ion_levels() with attach_atomic_columns(), then the shared HTML writer. PySCF is an optional dependency (pip install "CrystOD[quantum]").

Parameters:

cell – The crystal structure as phonopy.structure.atoms.PhonopyAtoms (converted to the spglib primitive cell).
left_tokens – --co-left formula tokens, e.g. ["Sc"].
right_tokens – --co-right formula tokens, e.g. ["F3"].
symprec – Symmetry tolerance handed to spglib.
basis – GTH basis set name (--basis); GTH_BASIS_SETS lists what PySCF ships, and the coverage of every element is checked.
pseudo – GTH pseudopotential family (--pseudo).
xc – Exchange-correlation functional (--xc; "hf" for Hartree-Fock).
kmesh – SCF k mesh [n1, n2, n3] (--kmesh); default default_kmesh(), round(8 Angstrom / |a_i|) per axis.
ke_cutoff – FFT density-grid cutoff in Hartree (--ke-cutoff); below about 80 the GTH Gaussians are not resolved.
oxidation – {element: formal charge} (--oxidation); default pymatgen’s guess. Sets both the fragment ion charges and the point charges of the removed sublattice.
electrons – Electrons per cell filled into the crystal column (default: the crystal cell’s own count).
sigma – Fermi smearing width in eV (--sigma; 0 = integer occupations; an odd-electron cell always smears).
degeneracy_tol – Seed window in eV for clustering degenerate levels (--degeneracy-tol; default DEGENERACY_SEED_EV).
no_ghost – Exclude the removed sublattice’s basis functions from the fragment calculations (--no-ghost).
symmetrize – Re-diagonalize the group-averaged Fock so grid-broken degeneracies come out exact (--no-symmetrize turns it off).
max_l – Drop basis shells with l above this from every element (--max-l); None keeps all.
projection – "lowdin" or "mulliken", the population measure of the compositions and sketch lobe sizes (--projection).
chk – Restart file path (--chk): written after the SCFs when missing, read (skipping them) when present. report_and_write() defaults to CHK_{formula}.chk.
onsite – Single-Hamiltonian mode (--onsite): only the crystal SCF runs, and the fragment columns are the per-shell on-site multiplets of the crystal Fock operator.
conventional – Draw the hover sketches in the conventional cell (display only).
conv_tol – SCF convergence threshold in Hartree.
max_cycle – SCF iteration limit.
max_memory – PySCF memory limit in MB.
verbose – PySCF verbosity level (--verbose).

Variables:

builder – The SymmetryAdaptedOrbitalBasis of the cell.
symbols – Chemical symbols of the primitive-cell atoms; positions their fractional and cartesian their Cartesian coordinates, lattice the lattice vectors as rows (Angstrom).
cells – {"mo": crystal, "left": ..., "right": ...} PySCF pbc.gto.Cell objects that share one AO space of n_ao functions.
formula – {"left": ..., "right": ...}, the fragment formulas.
side_atoms – Atom indices of each fragment; side_charge their formal charges, side_electrons their electron counts, and crystal_electrons that of the crystal.
oxidation – The formal charges in use.
specs – One AOBlock per (element, shell) of the AO space; side_specs[column] those of one fragment, ao_blocks the unmerged per-atom blocks.
kmesh – The SCF mesh in use.
mean_field – After run(): {column: converged KRKS/KRHF}; density_matrix and scf_energy (Hartree) likewise.
smeared – The columns whose occupations ended up Fermi-smeared.
last_coupling – Set by solve_at(): (left level, right level, |H~|, gap, minority weight, |S|, |H|) tuples of the same-irrep fragment pairs, strongest mixing first; last_gauge_residual the D+ S D - S residual of the representation check.
atomic_ions – After atomic_ion_levels(): {element: {"charge", "nelec", "method", "shells", "caveats"}}.
chk_path – The restart file in use, if any.

Raises:

ImportError – PySCF is not installed.
SystemExit – A fragment formula that does not match the cell, a basis or functional PySCF does not ship for these elements, a projection other than lowdin/mulliken, a kmesh entry below 1, or oxidation states that are not charge-neutral.

Example

The diagram of ScF3, three SCFs on a 2x2x2 mesh (minutes, not seconds):

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.PySCFCrystalOrbitalDiagram(cell, ["Sc"], ["F3"])
diagram.run()
kpoints = diagram.special_kpoints()
diagram.prepare_bands([k for _, k in kpoints])
levels, labels = diagram.solve_at(kpoints[0][1])   # GM

align_fragment_columns(records)[source]#

Deep-level (XPS-style) alignment of the three energy columns.

Each calculation carries its own G = 0 average-potential reference, so the raw columns are offset by one rigid constant each. For every fragment column the deepest chemically inert level is located – the deepest fragment level some crystal level consists of to at least ALIGNMENT_PURITY – and its fragment -> crystal energy difference, averaged over the k points where the pair exists, is that column’s offset. The zero is then put at the deeper of the two anchors in its PRE-BONDING (fragment) value: that fragment column stays, the crystal column moves by -delta_ref, the other fragment column by delta_other - delta_ref.

Parameters:: records – The per-k-point list built by report_and_write(): dicts with "name", "kpoint" and "levels" (the solve_at() result), one per diagram k point. The level energies are shifted in place.
Returns:: shifts maps "left", "mo" and "right" to the applied energy shifts in eV and "reference" to the anchoring column; anchors[column] describes each fragment’s anchor level ("label", "fragment_energy", "purity", "spread", "n_k", "fallback") or is None. (None, anchors) when no chemically inert fragment level exists; nothing is shifted then.
Return type:: (shifts, anchors)

atomic_ion_levels()[source]#

Levels of one isolated ion per element, at its formal charge.

The atomic stage before the sublattice forms (the crystal analogue of MolOD’s ligand-ao column), computed with PySCF as the three-stage story charged atom -> charged sublattice -> crystal.

Same basis / pseudopotential / functional as the periodic calculations (GTH pseudopotentials work in PySCF’s molecular code). A cation whose formal charge removes every pseudo-valence electron (Al^3+ with GTH-q3) has nothing to converge: its levels are the bare-ion one-electron spectrum (hcore eigenvalues – for Al^3+ the 3s eigenvalue, -27.9 eV, reproduces the third ionization potential of Al, 28.4 eV). Anions are vacuum-unbound (positive eigenvalues); that raw offset is absorbed by the per-element deep-shell anchoring in attach_atomic_columns, which also bridges the molecular (vacuum) and periodic (G = 0) energy references.

Returns:: {"charge", "nelec", "method", "shells": [(shell_name, l, energy_eV), ...], "caveats"}, one shell per (element, shell) spec of the AO basis.
Return type:: None. Fills atomic_ions with one entry per element

attach_atomic_columns(records)[source]#

Outermost isolated-ion columns + splitting connector links.

Call AFTER align_fragment_columns: the molecular (vacuum) and periodic (G = 0) references share no common zero, so each element’s ion levels are shifted rigidly so that its DEEPEST shell matches the (degeneracy-weighted, all-k) mean energy of the fragment levels that shell dominates – the band’s center of gravity, which in an orthogonal-basis tight-binding picture IS the on-site energy. The anchor shell’s connector fan then shows pure intra-sublattice splitting; the other shells additionally carry the ion’s own level spacing against the environment’s.

Parameters:: records – The per-k-point list of align_fragment_columns(); the "left-ao" and "right-ao" columns are added to every record’s "levels" in place.
Returns:: {element: (anchor shell, shift)} for the report, the shift in eV; ("none", 0.0) for an element none of whose shells dominates a fragment level.

compound_formula()[source]#

Reduced formula in conventional chemical order (rutile: TiO2).

The same helper the symmetry-mode tables are named after: cations before anions, the cation on the most special Wyckoff site first. Falls back to a plain alphabetical composition if pymatgen’s oxidation-state guesser has nothing to say about the elements. report_and_write() names the default restart file CHK_{formula}.chk after it.

Returns:: The formula string.
Return type:: str

little_group_data(kpoint)[source]#

Irreps, labels and the AO representation of the little group at k.

Same construction as the extended-Hueckel engine, but written directly in PySCF’s AO ordering: D[(a', shell, m'), (a, shell, m)] = P[a', a] W^l[m', m].

Parameters:: kpoint – Three primitive reciprocal coordinates.
Returns:: (irreps, mapping, labels, representation) as in CrystalOrbitalDiagram.little_group_data().

overlap_at(kpoint)[source]#

PySCF AO overlap matrix S(k) of the crystal cell.

Parameters:: kpoint – Three primitive reciprocal coordinates.
Returns:: Hermitian complex array of shape (n_ao, n_ao).
Return type:: ndarray

prepare_bands(kpoints)[source]#

Diagonalize every calculation at every diagram k point in one call.

get_bands rebuilds the density on the FFT grid each time it is called, so asking for all the special points at once instead of one per k point removes that cost from all but the first. Call after run(); solve_at() then reads the cache.

Parameters:: kpoints – List of three-component primitive reciprocal coordinates.
Returns:: None.
Return type:: None

run(report=<built-in function print>)[source]#

Run the periodic SCF calculations, or restore them from chk.

The crystal is solved first; its converged density restricted to one sublattice’s AO block is the initial guess of that fragment. A calculation that does not converge is retried up a ladder of Fermi smearing widths, cycle counts and virtual-level shifts, each rung reported. The converged results land in mean_field, density_matrix and scf_energy, and are written to chk_path when it is set. With onsite only the crystal calculation runs.

Parameters:: report – Callable that receives the progress lines (default print).
Returns:: None.
Raises:: SystemExit – An SCF that does not converge even with smearing, or an explicit chk file whose parameters do not match.
Return type:: None

site_symmetry_irreps(kpoint, representation, irreps, labels)[source]#

Irrep content of every (element, shell) block at k.

The site-symmetry induced representation that crystod --element EL --orbital ORB reports, recomputed here from the very representation used for the labelling.

Parameters:

kpoint – Three primitive reciprocal coordinates (not used by the computation; kept for symmetry with little_group_data()).
representation – The AO representation matrices from little_group_data().
irreps – The spgrep irreps from the same call.
labels – Their ISO-IR labels.

Returns:

the irreps the shell’s Bloch orbitals span, multiplicities as prefixes.

Return type:

{(element, shell): ["GM1+", "2GM4-", ...]}

sketch_partners(level, kpoint, sites)[source]#

Real wave-function amplitudes of a level on the supercell atoms.

The hover sketch of one level from the PySCF AO coefficients, in the entry format of the extended-Hueckel sketch (CrystalOrbitalDiagram.sketch_partners()).

The PySCF eigenvectors are in the atomic Bloch gauge, so the cell-to-cell phase is exp(2 pi i k . T) without the site offset. Same-l shells of one atom accumulate with their contracted-GTO radial amplitude at r0 = 2 bohr (see _build_ao_blocks), which fixes the lobe signs and orientation; the lobe SIZE is then rescaled to the Loewdin population of that (atom, l) channel. Raw r0 amplitudes would misstate the sizes: a diffuse gth-dzvp Sc 4p is ~5x an F 2p at r0, so a 7%-population Sc admixture used to draw at 83% of the largest F lobe. f shells and higher are omitted from the drawing.

Parameters:

level – A DiagramLevel from solve_at().
kpoint – The k point the level was solved at.
sites – The (atom index, translation) list of supercell_for() at that k point.

Returns:

solve_at(kpoint)[source]#

Solve the fragment and crystal levels at one k point from the SCFs.

Reads the band energies and coefficients cached by prepare_bands() (or computes them), clusters degenerate levels until their irrep multiplicities are integral, labels them, drops ghost-dominated fragment states, fills them by aufbau, attaches the population rows and the COOP bond characters, and records the same-irrep fragment couplings in last_coupling.

Parameters:: kpoint – Three primitive reciprocal coordinates.
Returns:: (levels, labels) with the columns "left", "mo" and "right" as in CrystalOrbitalDiagram.solve_at() (the -ao columns are added later by attach_atomic_columns()); energies in eV on the raw per-calculation references until align_fragment_columns() shifts them.
Raises:: SystemExit – The AO representation does not leave the PySCF overlap invariant (please report the case).

class crystod.salc.SymmetryAdaptedOrbitalBasis(cell, symprec=1e-05, standardize=True)[source]#

Bases: SymmetryOnlyVibrations

SALC bases of one element’s shell at a k point (crystod --visualize).

Where CrystalOrbital counts irreps from characters, this class builds the representation matrices themselves – for every little-group operation the Kronecker product of the Bloch-phased site permutation with the real-orbital Wigner matrix of the shell – and projects them onto the irreps with spgrep, which gives the symmetry-adapted linear combinations (SALCs) as explicit coefficient vectors over the (atom, m) orbital components. The command crystod --visualize -c POSCAR --element EL --orbital ORB --kpoint K prints those coefficients and writes the interactive 3D HTML viewer from them (s, p, d and f shells; without --kpoint one page per special k point of the space group).

The symmetry machinery is inherited from SymmetryOnlyVibrations (the crystod-phonon engine): the cell is reduced to the spglib primitive cell, get_irrep_labels() supplies the ISO-IR labels of the spgrep irreps, and resolve_qpoint the seekpath k-point labels that resolve_kpoint_input() relies on.

Parameters:

cell (PhonopyAtoms) – The crystal structure as phonopy.structure.atoms.PhonopyAtoms.
symprec (float) – Symmetry tolerance handed to spglib.
standardize (bool) – Convert cell to the spglib primitive cell (the default); False keeps it as given, which must then already be a primitive cell.

Variables:

primitive_cell – The PhonopyAtoms cell all atom indices refer to.
spglib_dataset – The spglib symmetry dataset of that cell ("international", "number", "wyckoffs", …).
rotations – Integer rotation matrices in the primitive basis, shape (n_ops, 3, 3), in spglib order.
translations – The matching fractional translations, (n_ops, 3).
rotations_cartesian – The same rotations as Cartesian matrices, the input of the Wigner matrices.
symprec – The symmetry tolerance in use.

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")
>>> basis = salc.SymmetryAdaptedOrbitalBasis(cell=cell)
>>> k = [0, 0, 0]
>>> irreps, rep, mapping, atoms = basis.get_orbital_rep(k, "F", l=1)
>>> labels = basis.get_irrep_labels(k, irreps, mapping)
>>> counts = basis.decompose_orbital_rep(irreps, rep, labels)
>>> {label: n for label, n in counts.items() if n > 0}
{'GM4-(3)': 2.0, 'GM5-(3)': 1.0}
>>> spaces, space_labels = basis.get_orbital_basis(irreps, rep, labels)
>>> [space.shape for space in spaces], space_labels
([(3, 9), (3, 9), (3, 9)], ['GM4-(3)', 'GM4-(3)', 'GM5-(3)'])

decompose_orbital_rep(irreps, orbital_rep, irrep_labels)[source]#

Multiplicity of every irrep in the orbital representation.

Parameters:

irreps – The spgrep irreps from get_orbital_rep().
orbital_rep – The representation matrices from the same call.
irrep_labels (list[str]) – One label per irrep (get_irrep_labels()).

Returns:

{label: multiplicity} for every irrep, the multiplicities rounded to two decimals (zero entries included). The printed * Irreducible Decomposition * of the CLI lists the non-zero ones.

Return type:

dict[str, float]

get_element_indices(element)[source]#

Indices of the atoms whose orbitals enter the basis.

Parameters:: element (str) – Chemical symbol, or "all" for every atom of the cell.
Returns:: The atom indices in primitive-cell order.
Raises:: ValueError – The element is not in the cell.
Return type:: list[int]

get_orbital_basis(irreps, orbital_rep, irrep_labels)[source]#

Project the orbital representation onto its irreps: the SALCs.

Parameters:

irreps – The spgrep irreps from get_orbital_rep().
orbital_rep – The representation matrices from the same call.
irrep_labels (list[str]) – One label per irrep (get_irrep_labels()).

Returns:

one complex array of shape (dim, n_atoms (2l+1)) per occurrence of an irrep – its dim rows are the partner SALCs, the columns the (atom, m) coefficients in the order of get_orbital_rep() – and the irrep label of every space (repeated when an irrep occurs more than once). --mode-index N of the CLI selects the N-th space, 1-based.

Return type:

(basis_spaces, basis_labels)

get_orbital_rep(kpoint, element, l)[source]#

Representation of the little group on the shell’s Bloch orbitals.

Parameters:

kpoint (list[float]) – Three primitive reciprocal coordinates.
element (str) – Chemical symbol, or "all".
l (int) – Azimuthal quantum number of the shell (0 to 3).

Returns:

the spgrep irreps at kpoint; the representation matrices, a complex array of shape (order, n_atoms (2l+1), n_atoms (2l+1)) whose rows and columns run atom-major over the (atom, m) components, m in the real-orbital order of ORBITAL_COMPONENT_NAMES; the indices into rotations of the little-group operations; and the atom indices of the element.

Return type:

(irreps, orbital_rep, mapping_little_group, element_indices)

crystod.salc.assign_bond_characters(levels, overlap, rows_left, rows_right, spec_ranges, sqrt_overlap=None, hamiltonian=None)[source]#

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+).

Parameters:

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.

Return type:

None

crystod.salc.compute_star(rotations, translations, kpoint)[source]#

Compute the star of k: the orbit of a k point under the rotations.

Every rotation R of the space group sends k to k' = k R; the distinct images modulo reciprocal-lattice vectors are the arms of the star, and the rotations sending k to one arm form a coset of the little group of k. This is the computation behind crystod --star-of-k; crystod-phonon --modulation uses it to combine the arms of a multi-q modulation.

Parameters:

rotations (ndarray[tuple[Any, ...], dtype[int64]]) – Integer rotation matrices of the space group in the primitive basis, shape (n_ops, 3, 3) – for example SymmetryAdaptedOrbitalBasis.rotations.
translations (ndarray[tuple[Any, ...], dtype[float64]]) – The matching fractional translations, shape (n_ops, 3). Accepted for a uniform call signature; the star depends on the rotations only.
kpoint (list[float]) – Three primitive reciprocal coordinates of k.

Returns:

"kpoint" (the arm wrapped into [-0.5, 0.5)), "representative_index" (index of the first rotation reaching the arm, the coset representative) and "operation_indices" (indices of every rotation reaching the arm). With the identity listed first, as spglib does, the first arm is k itself and its "operation_indices" are the little group of k.

Return type:

One dict per arm, in order of first appearance along rotations

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")
>>> basis = salc.SymmetryAdaptedOrbitalBasis(cell=cell)
>>> arms = salc.compute_star(basis.rotations, basis.translations,
...                          [0.5, 0.5, 0])
>>> [arm["kpoint"].tolist() for arm in arms]
[[0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]]
>>> len(arms[0]["operation_indices"])
16

crystod.salc.format_star_lines(arms, seitz_symbols=None, indent=' ')[source]#

Format the arms of a star as the report lines of crystod --star-of-k.

Parameters:

arms (list[dict]) – The list returned by compute_star().
seitz_symbols (list[str] | None) – Seitz symbols of all rotations, indexed like the rotations given to compute_star(); when present, each line names the coset-representative operation.
indent (str) – Text put in front of every line.

Returns:

One string per arm, "arm 1: k = [+0.5, +0.5, +0]" and so on, with "(representative: 2_001)" appended when symbols are given.

Return type:

list[str]

Example

>>> for line in salc.format_star_lines(arms):   # arms of compute_star
...     print(line)
 arm 1: k = [+0.5, +0.5, +0]
 arm 2: k = [+0.5, +0, +0.5]
 arm 3: k = [+0, +0.5, +0.5]

crystod.salc.resolve_kpoint_input(structure, raw_kpoint)[source]#

Resolve a --kpoint argument into a label and primitive coordinates.

The two command-line spellings of a k point – one high-symmetry label such as GM/X/M/R, or three primitive reciprocal coordinates (fractions such as 1/2 allowed) – become (label, coordinates) with the seekpath labels of the structure: coordinates that hit a special point (or an arm of its star) receive that point’s name, any other point the ISO-IR k-vector type letters (GP for a general point). This is what crystod --star-of-k and crystod --visualize do with --kpoint. When the label lookup itself fails (seekpath unavailable), three coordinates are still accepted and labelled custom.

Parameters:

structure (SymmetryOnlyVibrations) – A SymmetryOnlyVibrations of the cell, or a subclass such as SymmetryAdaptedOrbitalBasis; its resolve_qpoint does the work.
raw_kpoint (list[str]) – The tokens as typed: ["M"] or ["1/2", "1/2", "0"].

Returns:

the k-point name and its three primitive reciprocal coordinates as floats.

Return type:

(label, kpoint)

Raises:

ValueError – An unknown label, or a token count that is neither one nor three (a ValueError from the implementation module as well as through crystod.salc).

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")
>>> basis = salc.SymmetryAdaptedOrbitalBasis(cell=cell)
>>> salc.resolve_kpoint_input(basis, ["M"])
('M', [0.5, 0.5, 0.0])
>>> salc.resolve_kpoint_input(basis, ["1/2", "1/2", "0"])
('M', [0.5, 0.5, 0.0])

crystod.salc

Contents

crystod.salc#