r"""Quantitative crystal-orbital diagrams via PySCF PBC
(``crystod --diagram --pyscf``).
The crystalline counterpart of ``crystod-mol --diagram --pyscf``. The
symmetry + extended-Hueckel diagram of :mod:`crystod.crystal_orbital_diagram`
is replaced by three periodic self-consistent-field calculations that share
one atomic-orbital space:
======================= ============================== ==================
calculation real atoms point charges
======================= ============================== ==================
left fragment the ``--co-left`` sublattice the right sublattice
right fragment the ``--co-right`` sublattice the left sublattice
crystal everything none
======================= ============================== ==================
The removed sublattice stays in the basis as **ghost atoms** (basis functions
without a nucleus), so all three calculations span the same AO space and the
crystal orbitals can be projected exactly onto the fragment Bloch orbitals --
the counterpoise-consistent construction of the molecular version, carried
over to a periodic system. In addition the removed sublattice acts on the
fragment through its **formal-charge point lattice**, evaluated as a
jellium-referenced Ewald/FFT potential, so each fragment is the electronic
state of one sublattice in the Madelung field of the other: the state before
chemical bond formation.
Charges
-------
Both the ions that remain and the point charges that replace the removed
sublattice carry the formal oxidation states (``--oxidation`` to override).
For ScF3 that makes the left fragment Sc(3+) with three F(-1) point charges
and the right fragment 3 F(-1) with one Sc(3+) point charge, so
* every cell is **electrically neutral**, which removes the monopole
divergence a charged fragment array would introduce (the residual
constant offset between the calculations is handled by the deep-level
alignment below);
* every electron count is even, so a restricted (KRKS/KRHF) treatment is
consistent;
* the fragment electron counts add up to the electron count of the crystal.
k points
--------
A crystal-orbital diagram needs the crystal orbitals at the high-symmetry
points where the irreducible representations are tabulated -- not a full band
path. The diagram is therefore built at the special points of the space
group (the same list as ``crystod-bz --show-kpoint``), and the SCF runs on a
small regular mesh whose default is taken from the lattice constants,
n_i = round(8 Angstrom / |a_i|) (``--kmesh`` to override).
Energy reference (deep-level alignment)
---------------------------------------
The eigenvalues of a periodic calculation have no absolute zero: each of the
three calculations pins the G = 0 (cell-averaged) Coulomb potential to zero,
and although the neutral cells remove the monopole divergence, the *value* of
the average potential still depends on the second moment of the cell's charge
density -- which changes when an ion is replaced by a bare point charge. The
result is one rigid, k-independent offset per calculation (for ScF3 the Sc
column sits ~2 eV and the F3 column ~5 eV below the crystal column), exactly
the reference problem familiar from band-offset calculations.
The columns are therefore aligned the XPS way: a **chemically inert deep
level** -- the deepest fragment level that reappears in the crystal almost
unchanged (counterpart purity >= 80%) -- must have the same energy before and
after bond formation. The energy zero is kept at that anchor's *pre-bonding*
(fragment) value, i.e. the crystal column and the other fragment column are
shifted onto the frame of the deepest unhybridized sublattice level; the
constants are printed, and ``--no-align`` disables the whole step.
Irreducible representations
---------------------------
Every level -- fragment and crystal -- is labelled with the little-group irrep
of its Bloch state, using crystod's own machinery: the site-permutation
representation at k combined with the real-orbital Wigner matrices, projected
with the spgrep characters. The representation is verified against PySCF's
own overlap matrix (D+ S D = S) at every k point before it is used. Levels
of the same irrep are connected in the diagram, levels of different irreps are
not; the strength of each connection is the projection of the crystal orbital
onto the fragment orbital, and the fragment-fragment coupling matrix element
<phi_left| F(k) |phi_right> of the converged crystal Fock operator is reported
alongside, since two levels of the same irrep only mix appreciably when that
matrix element is large compared with their energy separation.
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
from dataclasses import dataclass
import numpy as np
from .crystal_orbital_diagram import (
CrystalOrbitalDiagram,
DiagramLevel,
_composition_string,
_format_kpoint,
assign_bond_characters,
assign_fragment_compositions,
parse_fragment_formula,
write_crystal_diagram_html,
)
from .operations import wigner_D_real
from .runtime_compat import get_chemical_symbols, get_scaled_positions
from .visualize_basis import SymmetryAdaptedOrbitalBasis
HARTREE_TO_EV = 27.211386245988
BOHR_TO_ANGSTROM = 0.52917721092
# Default target for the automatic SCF mesh: n_i = round(K / |a_i|) with the
# lattice constant in Angstrom. K = 8 reproduces 2x2x2 for the ~4 Angstrom
# perovskite-like cells this module is written for.
KMESH_TARGET_ANGSTROM = 8.0
# A fragment level whose Mulliken population sits mostly on the ghost basis of
# the *other* sublattice is a counterpoise artifact, not a state of the
# fragment, and is dropped (same rule as crystod-mol --pyscf).
GHOST_FRACTION_THRESHOLD = 0.35
# --onsite block diagonalization: sublattice Bloch combinations whose overlap
# eigenvalue falls below this are dropped (canonical orthogonalization), the
# same near-dependence guard the SCF itself applies to the full basis
ONSITE_OVERLAP_FLOOR = 1.0e-7
# Seed window in eV for clustering degenerate levels. The uniform FFT grid on
# which the Coulomb, the exchange-correlation and the point-charge potentials
# are evaluated does not respect the point group exactly, so levels that are
# degenerate by symmetry come out split -- by ~5e-4 eV for the occupied states
# of ScF3, but by tens of meV for the empty states of ZrO2 at a coarse cutoff.
# No fixed window can cover that range without merging genuinely distinct
# levels, so this is only the starting guess: groups are then merged until the
# irrep multiplicities come out integral (see _adaptive_groups), which is the
# statement that the group is a whole number of complete multiplets.
DEGENERACY_SEED_EV = 1.0e-3
# A group is never widened beyond this, so two levels that are really distinct
# can never be merged just because the projection is noisy.
DEGENERACY_MAX_WINDOW_EV = 0.30
# How far a multiplicity may sit from an integer and still count as one.
_MULTIPLICITY_TOL = 0.05
# The GTH basis sets PySCF ships, in the order --basis lists them (the
# testsuite checks this against pyscf.pbc.gto.basis.ALIAS). Apart from W in
# gth-dzvp, only the two molopt-sr sets cover the transition metals (and none
# of them the lanthanides); the rest stop at Ar or cover light elements only.
GTH_BASIS_SETS = (
"gth-szv-molopt-sr", "gth-dzvp-molopt-sr",
"gth-szv", "gth-dzv", "gth-dzvp", "gth-tzvp", "gth-tzv2p",
"gth-qzv2p", "gth-qzv3p",
"gth-szv-molopt", "gth-dzvp-molopt", "gth-tzvp-molopt",
"gth-tzv2p-molopt",
"gth-aug-dzvp", "gth-aug-tzvp", "gth-aug-tzv2p", "gth-aug-qzv2p",
"gth-aug-qzv3p",
"gth-cc-dzvp", "gth-cc-tzvp", "gth-cc-qzvp",
)
# The GTH pseudopotentials PySCF ships (same check in the testsuite).
GTH_PSEUDOPOTENTIALS = (
"gth-pbe", "gth-pade", "gth-lda", "gth-hfrev", "gth-blyp", "gth-bp",
"gth-hcth120", "gth-hcth407", "gth-hf", "gth-olyp", "gth-pbesol",
)
# The diffuse-richer basis the empty-shell caveats point at (light elements
# only -- offered only when it covers every element of the structure).
RICHER_BASIS = "gth-qzv2p"
# A fragment level counts as chemically inert -- usable as an alignment anchor
# -- when some crystal level consists of it to at least this fraction.
ALIGNMENT_PURITY = 0.80
# crystod orders the real orbitals as in complex_to_real_transform_orbital;
# PySCF uses m = -l..l except for p, which it orders x, y, z. The two agree
# for s, p and d and differ for f and above, so the AO rotation is permuted
# into PySCF's order.
_CRYSTOD_M_ORDER = {
0: [0],
1: [1, -1, 0],
2: [-2, -1, 0, 1, 2],
3: [3, -3, 2, -2, 1, -1, 0],
}
def _import_pyscf():
"""Import PySCF, or raise ``ImportError`` naming the ``[quantum]`` extra.
PySCF is an optional dependency (``pip install "CrystOD[quantum]"``);
the command-line front end reports the same condition as a one-line
``ERROR:`` before dispatching here, and the library API lets the
``ImportError`` propagate.
"""
from ._optional import require_pyscf
require_pyscf("crystod --diagram/--band/--dos/--visualize --pyscf "
"(PySCFCrystalOrbitalDiagram)")
from pyscf.pbc import dft, gto, scf, tools # noqa: F401
# some conda builds of pyscf default to a SINGLE OpenMP thread unless
# OMP_NUM_THREADS is exported; that turns a minutes-long diagram into an
# hours-long one. Use every core unless the user chose otherwise.
import os
from pyscf import lib
if not os.environ.get("OMP_NUM_THREADS"):
lib.num_threads(os.cpu_count())
def _pyscf_m_order(l: int) -> list[int]:
return [1, -1, 0] if l == 1 else list(range(-l, l + 1))
def _reorder_to_pyscf(l: int) -> np.ndarray:
"""Q with Q[pyscf_row, crystod_row] = 1 for the real orbitals of shell l."""
crystod_order = _CRYSTOD_M_ORDER.get(l, list(range(-l, l + 1)))
pyscf_order = _pyscf_m_order(l)
matrix = np.zeros((2 * l + 1, 2 * l + 1))
for column, m in enumerate(crystod_order):
matrix[pyscf_order.index(m), column] = 1.0
return matrix
def wigner_pyscf(l: int, rotation: np.ndarray) -> np.ndarray:
"""Real-orbital rotation matrix in PySCF's AO ordering."""
transform = _reorder_to_pyscf(l)
return transform @ wigner_D_real(l, rotation) @ transform.T
def default_kmesh(lattice: np.ndarray) -> list[int]:
"""n_i = round(8 Angstrom / |a_i|), at least 1 -- the user-facing rule of
thumb that a ~4 Angstrom cell wants a 2x2x2 mesh."""
lengths = np.linalg.norm(np.asarray(lattice, dtype=float), axis=1)
return [max(1, int(round(KMESH_TARGET_ANGSTROM / length))) for length in lengths]
@dataclass
class AOBlock:
"""One (element, shell) block of the PySCF AO space, for level labels.
Attributes:
element: Chemical symbol.
shell: Shell name such as ``"3d"``.
l: Azimuthal quantum number.
offset: First AO index of the block.
n_ao: Number of AOs in the block (``2l+1`` per site).
sites: Indices of the atoms carrying the shell.
column: ``"left"`` or ``"right"``.
radial: Radial amplitude of the contracted function at the sketch
radius ``r0``.
radial_profile: The same at every sketch radius, for the
multi-radius sign of the viewer.
"""
element: str
shell: str
l: int
offset: int
n_ao: int
sites: list[int]
column: str
radial: float = 1.0 # radial amplitude at the sketch radius r0
# amplitudes at SKETCH_RADII, for the multi-radius sign of the viewer:
# a single probe radius can sit right on the orthogonalization node of a
# semicore-carrying channel (Sc s of the valence R1+ of ScF3: -0.0013 at
# 2 bohr), where the drawn sign would be decided by noise
radial_profile: tuple = ()
[docs]
class PySCFCrystalOrbitalDiagram(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,
:func:`report_and_write`, is :meth:`run` (the SCFs, or a ``chk``
restart), :meth:`prepare_bands` for all :meth:`special_kpoints`,
:meth:`solve_at` and :meth:`site_symmetry_irreps` per k point,
:meth:`align_fragment_columns`, :meth:`atomic_ion_levels` with
:meth:`attach_atomic_columns`, then the shared HTML writer. PySCF is
an optional dependency (``pip install "CrystOD[quantum]"``).
Args:
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
:func:`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.
:func:`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``).
Attributes:
builder: The :class:`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 :meth:`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 :meth:`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 :meth:`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
"""
def __init__(self, cell, left_tokens, right_tokens, *, symprec=1e-5,
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-8, max_cycle=100,
max_memory=4000.0, verbose=0):
_import_pyscf()
# --conventional: draw the hover sketches in the conventional cell
# (display-only -- deliberately NOT part of _chk_params)
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)
self.cartesian = self.positions @ self.lattice
self.basis_name = basis
self.pseudo_name = pseudo
self.xc = xc
self.ke_cutoff = ke_cutoff
if ke_cutoff and ke_cutoff < 80.0:
# measured on ScF3 (dzvp, 2x2x2): below ~80 Hartree the FFT grid
# cannot resolve the compact Gaussians (semicore shells) and the
# SCF does not converge at all -- the grid here is the DENSITY
# grid, the analogue of VASP's augmentation grid (~4x ENCUT in
# energy), not of ENCUT itself.
print(f"WARNING: --ke-cutoff {ke_cutoff:g} Hartree is below the "
"~80 Hartree the GTH Gaussian basis needs; expect the SCF "
"to fail. 100-200 Hartree is the working range.")
self.sigma = float(sigma)
self.no_ghost = bool(no_ghost)
self.symmetrize = bool(symmetrize)
self.max_l = max_l
if projection not in ("lowdin", "mulliken"):
raise SystemExit(
"ERROR: --projection must be 'lowdin' or 'mulliken', "
f"not {projection!r}.")
self.projection = projection
# display label for the population rows ("Loewdin: ..."/"Mulliken: ...")
self.projection_label = ("Loewdin" if projection == "lowdin"
else "Mulliken")
self.chk_path = chk
# True when report_and_write derived CHK_{formula}.chk itself:
# an automatic checkpoint is a cache, so a stale or mismatched
# one is recomputed and overwritten instead of being an error
self.chk_auto = False
self.degeneracy_tol = (DEGENERACY_SEED_EV if degeneracy_tol is None
else float(degeneracy_tol))
self.degeneracy_window = DEGENERACY_MAX_WINDOW_EV
self.retry_sigma = 0.2
self.smeared: set[str] = set()
self.conv_tol = conv_tol
self.max_cycle = max_cycle
self.max_memory = max_memory
self.pyscf_verbose = verbose
self.sketch_specs = None # sketches always use all AO components
self.sketch_tokens = None
self.method_chip = f"PySCF {xc.upper()}/{basis}"
self._check_basis_coverage()
self._check_xc()
# a diffuse-richer set to recommend in the empty-shell caveats -- only
# if PySCF ships it for every element here (gth-qzv2p stops at Ar, so
# a transition-metal compound has no richer alternative at all)
self.richer_basis = (
RICHER_BASIS
if ("molopt" in self.basis_name.lower()
and self._basis_covers(RICHER_BASIS))
else "")
# footer sentence for the outermost isolated-ion columns (the shared
# writer's default describes the EHT VSIP columns instead)
self.ao_foot = (
" The outermost columns are each element's ISOLATED ion at its "
"formal charge -- one PySCF calculation per element with the "
"same basis, pseudopotential and functional (RKS/UKS; a cation "
"left with no pseudo-valence electrons keeps its bare-ion "
"one-electron spectrum) -- rigidly shifted per element so its "
"deepest shell with a fragment counterpart sits at that "
"shell's sublattice band center (the molecular vacuum and "
"periodic G=0 references share no common zero; the shift and "
"the raw vacuum levels are in each level's tooltip). The "
"connector lines into the sublattice columns show how the "
"ions' Bloch combinations split each shell at the chosen "
"k point. Shells the formal charge leaves EMPTY come with a "
"caveat (each affected tooltip says so): for a cation only "
"the lowest empty state of each l is trustworthy -- the "
"higher ones are finite-basis virtuals whose order and "
"energy follow the basis"
+ (", and the default condensed-phase basis has no diffuse "
"functions (a diffuse-richer basis such as --basis "
f"{self.richer_basis} removes that basis-side error, at "
"the cost of a heavier and possibly ill-conditioned "
"periodic SCF)"
if self.richer_basis else "")
+ "; every empty shell of an ANION is instead a discretized "
"continuum state that no basis makes physical (one more "
"electron on a free anion is vacuum-unbound -- in the "
"crystal it is the Madelung potential that binds)."
)
# --onsite: single-Hamiltonian mode. Only the crystal SCF runs; the
# fragment columns are the sublattice BLOCKS of its converged Fock,
# F[rows,rows] c = E S[rows,rows] c -- the pre-bonding sublattice
# states with the left-right mixing switched off. One operator for
# every column, so no reference alignment is needed and a level's
# rise/drop against its parents is purely the orbital interaction.
self.onsite = bool(onsite)
self.scf_columns = ("mo",) if self.onsite else ("mo", "left", "right")
if self.onsite and self.no_ghost:
print("NOTE: --onsite ignores --no-ghost (no fragment SCF runs; "
"the columns are crystal-Fock blocks).")
# canonicalize: the flag is numerically inert here, and keeping
# it would poison the chk parameter check for no reason
self.no_ghost = False
self.basis_chip = "(GTH valence basis)"
if self.onsite:
self.embedding_chip = "one Hamiltonian: crystal-Fock blocks"
self.foot_intro = (
"Crystal-orbital diagram (COD, --onsite): every column comes "
"from the ONE converged crystal Fock operator F(k). The "
"fragment columns are the per-(element, shell) ON-SITE "
"multiplets -- F(k) diagonalized within each shell's own "
"symmetry-adapted Bloch orbitals, the tight-binding on-site "
"energies <φ|F|φ> of the actual AO shells, one "
"level per induced irrep and no cross-shell mixing -- the "
"center its full eigenstates; states sharing an irrep of the "
"little group at k mix into bonding/antibonding crystal "
"orbitals. No fragment SCF and no reference alignment: a "
"crystal level's drop/rise against its parents is the "
"orbital interaction (level repulsion/hybridization) with "
"everything else, cross-shell and left–right alike. "
"Column occupations (electron arrows) are the formal ionic "
"counts from the oxidation states -- a display convention, "
"not an output of the Fock operator.")
else:
self.embedding_chip = ""
self.foot_intro = (
"Crystal-orbital diagram (COD): the fragment-sublattice "
"Bloch orbitals (columns) are the electronic states before "
"chemical bond formation -- each fragment is its own "
"periodic DFT calculation (formal-charge ions + ghost basis "
"+ point-charge lattice of the removed sublattice); states "
"sharing an irrep of the little group at k mix into "
"bonding/antibonding crystal orbitals (center). Energies: "
f"PySCF {xc.upper()} eigenvalues, the three calculations put "
"on one scale by deep-level (XPS-style) alignment. NOTE the "
"parent→crystal vertical offsets also carry the "
"point-charge-model-vs-crystal environment difference "
"(site-dependent, up to ~1.5 eV) -- read bonding from the "
"line colors (COOP), not from the offsets; --onsite removes "
"this by drawing every column from the crystal Fock itself.")
self._assign_fragments(left_tokens, right_tokens)
self._resolve_oxidation(oxidation)
if not self.onsite and not any(self.oxidation.values()):
# neutral sublattices (--oxidation El=0 ...): no point charges
self.foot_intro = self.foot_intro.replace(
"(formal-charge ions + ghost basis + point-charge lattice "
"of the removed sublattice)",
"(neutral sublattices + ghost basis; all oxidation states "
"are 0, so there is no point-charge lattice)").replace(
"point-charge-model-vs-crystal environment difference",
"fragment-model-vs-crystal environment difference")
self.kmesh = list(kmesh) if kmesh else default_kmesh(self.lattice)
if any(n < 1 for n in self.kmesh):
raise SystemExit("ERROR: every --kmesh entry must be at least 1.")
# the crystal cell first: with --no-ghost the fragment cells span only
# their own sublattice's AOs and are embedded into the crystal's AO
# space through the atom slices of the crystal cell
self.odd_electron: set[str] = set()
self.cells = {"mo": self._make_cell(None)}
self.n_ao = int(self.cells["mo"].nao_nr())
slices = self.cells["mo"].aoslice_by_atom()
self.embed_rows = {
column: np.concatenate([
np.arange(int(slices[atom, 2]), int(slices[atom, 3]))
for atom in self.side_atoms[column]
])
for column in ("left", "right")
}
self.cells["left"] = self._make_cell("left")
self.cells["right"] = self._make_cell("right")
self._check_shared_ao_space()
self._build_ao_blocks()
self.side_electrons = {
column: int(self.cells[column].nelectron) for column in ("left", "right")
}
crystal_electrons = int(self.cells["mo"].nelectron)
if electrons is None:
electrons = crystal_electrons
self.electrons = float(electrons)
self.crystal_electrons = crystal_electrons
self.mean_field = {}
self.density_matrix = {}
self.scf_energy = {}
self._pc_potential = {}
# ------------------------------------------------------------ fragments
def _assign_fragments(self, left_tokens, right_tokens) -> None:
formulas = {
"left": parse_fragment_formula(left_tokens, "--co-left"),
"right": parse_fragment_formula(right_tokens, "--co-right"),
}
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 column in ("left", "right"):
for element, count in formulas[column]:
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[column]} lists {element}{count} but the "
f"primitive cell has {composition[element]} {element} "
f"atom(s) (composition: {comp_str})."
)
assigned[element] = column
missing = [element for element in composition if element not in assigned]
if missing:
raise SystemExit(
f"ERROR: element(s) {', '.join(missing)} not assigned to "
f"--co-left/--co-right (primitive-cell composition: {comp_str}; "
"every atom must belong to one fragment)."
)
self.element_column = assigned
self.atom_column = [assigned[symbol] for symbol in self.symbols]
self.side_atoms = {
column: [index for index, side in enumerate(self.atom_column) if side == column]
for column in ("left", "right")
}
self.formula = {
column: _composition_string([self.symbols[i] for i in self.side_atoms[column]])
for column in ("left", "right")
}
def _resolve_oxidation(self, oxidation) -> None:
composition: dict[str, int] = {}
for symbol in self.symbols:
composition[symbol] = composition.get(symbol, 0) + 1
if oxidation is None:
from pymatgen.core import Composition
guesses = Composition(_composition_string(self.symbols)).oxi_state_guesses()
if not guesses:
raise SystemExit(
"ERROR: could not guess the oxidation states of "
f"{_composition_string(self.symbols)}; pass them explicitly, "
"e.g. --oxidation Sc=+3 F=-1."
)
oxidation = {element: float(q) for element, q in guesses[0].items()}
missing = [element for element in composition if element not in oxidation]
if missing:
raise SystemExit(f"ERROR: --oxidation misses element(s) {', '.join(missing)}.")
net = sum(oxidation[element] * count for element, count in composition.items())
if abs(net) > 1e-6:
raise SystemExit(
f"ERROR: oxidation states are not charge-neutral (net {net:+g} per cell)."
)
self.oxidation = oxidation
self.side_charge = {
column: sum(self.oxidation[self.symbols[i]] for i in self.side_atoms[column])
for column in ("left", "right")
}
# ---------------------------------------------------------- pyscf cells
def _make_cell(self, column):
"""Cell in which only ``column``'s atoms are real; None = the crystal.
By default every atom keeps its basis functions (ghosts elsewhere), so
all three cells span one AO space (counterpoise-consistent). With
``no_ghost`` the fragment cells contain no trace of the removed
sublattice at all -- its basis functions are excluded from the
variational space, which is the hard constraint that no fragment
wave function can sit on the removed atoms.
"""
from pyscf.pbc import gto
atoms, basis, pseudo = [], {}, {}
for index, symbol in enumerate(self.symbols):
real = column is None or self.atom_column[index] == column
if not real and self.no_ghost:
continue
tag = f"{symbol}{index}" if real else f"ghost-{symbol}{index}"
shells = gto.basis.load(self.basis_name, symbol)
if self.max_l is not None:
shells = [shell for shell in shells if shell[0] <= self.max_l]
if not shells:
raise SystemExit(
f"ERROR: --max-l {self.max_l} removes every shell of "
f"{symbol}.")
basis[tag] = shells
if real:
pseudo[tag] = self.pseudo_name
atoms.append((tag, tuple(self.cartesian[index])))
cell = gto.Cell()
cell.a = self.lattice
cell.atom = atoms
cell.unit = "Angstrom"
cell.basis = basis
cell.pseudo = pseudo
# formal-charge ions: the fragment carries the total oxidation state of
# its own sublattice, so together with the point charges of the removed
# one the cell is neutral
cell.charge = 0 if column is None else int(round(self.side_charge[column]))
cell.ke_cutoff = self.ke_cutoff
cell.max_memory = self.max_memory
cell.verbose = self.pyscf_verbose
# required by make_kpts(space_group_symmetry=True): the SCF then only
# solves the irreducible wedge of the mesh (2x2x2 of Pm-3m: 4 of 8)
cell.space_group_symmetry = True
cell.symmorphic = False
import warnings
with warnings.catch_warnings():
# neutral-sublattice runs (--oxidation Al=0 N=0) legitimately
# build odd-electron cells; the spin-0 inconsistency pyscf warns
# about is resolved below by forcing Fermi smearing
warnings.filterwarnings(
"ignore", message="Electron number .* not consistent")
cell.build()
# pin the PER-CELL electron count: pyscf's tot_electrons(nkpts) is
# atom_charges().sum()*nkpts - charge -- cell.charge is subtracted
# once for the whole Born-von-Karman supercell -- so a charged
# fragment on a k-mesh would keep charge*(nkpts-1)/nkpts spurious
# electrons per cell (Sc^+3 of ScF3 on 2x2x2 converged with 10.5
# instead of 8). With _nelectron set, tot_electrons(nkpts) is
# exactly nelectron*nkpts.
cell.nelectron = cell.nelectron
if cell.nelectron % 2:
# an odd count per cell needs fractional occupations in this
# spin-restricted driver: across the k-mesh pyscf's integer
# aufbau drops the unpaired electron when the BvK total is odd
# and hunts an unstable metallic degeneracy edge when it is
# even -- run() upgrades such a calculation to Fermi smearing,
# which conserves the count exactly
self.odd_electron.add("mo" if column is None else column)
return cell
def _check_shared_ao_space(self) -> None:
for column in ("left", "right"):
expected = (len(self.embed_rows[column]) if self.no_ghost
else self.n_ao)
actual = int(self.cells[column].nao_nr())
if actual != expected:
raise SystemExit(
f"ERROR: the {column} fragment spans {actual} AOs, expected "
f"{expected}; this is a bug, please report it."
)
def _embed(self, column, coefficients):
"""Zero-pad small-basis fragment eigenvectors into the shared AO space
(identity without --no-ghost)."""
if column == "mo" or not self.no_ghost:
return coefficients
full = np.zeros((self.n_ao, coefficients.shape[1]),
dtype=coefficients.dtype)
full[self.embed_rows[column]] = coefficients
return full
def _build_ao_blocks(self) -> None:
"""(element, shell) blocks of the PySCF AO space, in AO order."""
from pyscf.gto.mole import gto_norm
cell = self.cells["mo"]
slices = cell.aoslice_by_atom()
raw_labels = [label.split() for label in cell.ao_labels()]
blocks: list[AOBlock] = []
index = 0
sketch_r0 = 2.0 # bohr, the representative bonding-region radius
sketch_radii = (1.5, 2.0, 2.5, 3.0)
for shell in range(cell.nbas):
atom = int(cell.bas_atom(shell))
l = int(cell.bas_angular(shell))
exponents = np.asarray(cell.bas_exp(shell), dtype=float)
# bas_ctr_coeff is over unit-normalized primitives; the gto_norm
# factor makes the amplitudes below the actual AO values
contractions = (np.asarray(cell.bas_ctr_coeff(shell), dtype=float)
* gto_norm(l, exponents)[:, None])
# radial amplitude of each contracted function at r0, for the
# wave-function sketch: same-l shells of one atom accumulate with
# these weights, so the drawn lobe signs are those of the real
# wave function there (a bare coefficient of one shell would be
# wrong for semicore levels, whose orthogonalization tails invert)
radial_amplitudes = (
sketch_r0 ** l
* (contractions * np.exp(-exponents * sketch_r0 ** 2)[:, None]).sum(axis=0)
)
profiles = [
radius ** l
* (contractions * np.exp(-exponents * radius ** 2)[:, None]).sum(axis=0)
for radius in sketch_radii
]
for contraction in range(cell.bas_nctr(shell)):
# the AO label carries the shell name, e.g. '0 Sc 3dxy' -> '3d'
name = raw_labels[index][2]
shell_name = name[: len(name) - len(name.lstrip("0123456789"))] + "spdfgh"[l]
blocks.append(AOBlock(
element=self.symbols[atom], shell=shell_name, l=l,
offset=index, n_ao=2 * l + 1, sites=[atom],
column=self.atom_column[atom],
radial=float(radial_amplitudes[contraction]),
radial_profile=tuple(
float(profile[contraction]) for profile in profiles),
))
index += 2 * l + 1
self.ao_blocks = blocks
self.ao_offset_in_atom = {
block.offset: block.offset - int(slices[block.sites[0], 2]) for block in blocks
}
# merged (element, shell) specs for the level labels and the report
merged: dict[tuple[str, str], AOBlock] = {}
self.side_specs = {"left": [], "right": []}
for block in blocks:
key = (block.element, block.shell)
if key in merged:
merged[key].sites.append(block.sites[0])
continue
spec = AOBlock(block.element, block.shell, block.l, block.offset,
block.n_ao, list(block.sites), block.column)
merged[key] = spec
self.specs = list(merged.values())
for spec in self.specs:
self.side_specs[spec.column].append(spec)
# per-(element, shell) AO index lists, for Mulliken weights
self.spec_indices = {}
for (element, shell), spec in merged.items():
indices: list[int] = []
for block in blocks:
if block.element == element and block.shell == shell:
indices.extend(range(block.offset, block.offset + block.n_ao))
self.spec_indices[id(spec)] = np.array(indices, dtype=int)
# --------------------------------------------------- point-charge field
def _basis_covers(self, name) -> bool:
"""Does PySCF ship this GTH basis for every element of the cell?"""
from pyscf.pbc.gto import basis as pbc_basis
for element in set(self.symbols):
try:
pbc_basis.load(name, element)
except Exception:
return False
return True
def _check_basis_coverage(self) -> None:
"""Refuse a basis/pseudopotential that has no entry for one of the
elements, naming the sets that do -- PySCF's own BasisNotFoundError
arrives deep inside cell.build() and says nothing about the
alternatives (apart from W in gth-dzvp, only the two molopt-sr sets
reach the transition metals, and none of them the lanthanides)."""
from pyscf.pbc.gto import basis as pbc_basis
from pyscf.pbc.gto import pseudo as pbc_pseudo
for kind, name, loader, shipped in (
("basis", self.basis_name, pbc_basis.load, GTH_BASIS_SETS),
("pseudopotential", self.pseudo_name, pbc_pseudo.load,
GTH_PSEUDOPOTENTIALS)):
missing = []
for element in dict.fromkeys(self.symbols):
try:
loader(name, element)
except Exception:
missing.append(element)
if not missing:
continue
if name not in shipped:
raise SystemExit(
f"ERROR: PySCF has no {kind} '{name}' for "
f"{', '.join(missing)}.\n"
f" The GTH sets it ships: {', '.join(shipped)}")
usable = []
for candidate in shipped:
try:
for element in dict.fromkeys(self.symbols):
loader(candidate, element)
except Exception:
continue
usable.append(candidate)
option = "--basis" if kind == "basis" else "--pseudo"
advice = (f" {option} sets covering every element here: "
f"{', '.join(usable)}" if usable else
f" No GTH {kind} PySCF ships covers every "
"element of this structure.")
raise SystemExit(
f"ERROR: the {kind} '{name}' has no entry for "
f"{', '.join(missing)}.\n{advice}")
def _check_xc(self) -> None:
"""Reject an unusable --xc up front. libxc raises a bare KeyError
for a name it does not know (``pz``, ``vwn`` -- common shorthands
that are not libxc names), and PySCF's PERIODIC code has no
nonlocal-correlation path, so a VV10 functional (wb97m-v, b97m-v,
wb97x-v) dies inside get_veff with 'KNumInt has no attribute
nr_nlc_vxc' after the cells are already built."""
if self.xc.lower() in {"hf", "hartree-fock"}:
return
from pyscf.dft import libxc
try:
libxc.parse_xc(self.xc)
except KeyError:
raise SystemExit(
f"ERROR: '{self.xc}' is not a functional libxc knows.\n"
" crystod --help lists the verified names (note that "
"the LDA shorthands are 'lda' and 'svwn', not 'pz'/'vwn').")
if any(abs(b) > 0 or abs(c) > 0
for (b, c), _ in libxc.nlc_coeff(self.xc)):
raise SystemExit(
f"ERROR: '{self.xc}' carries VV10 nonlocal correlation, "
"which PySCF's periodic\n code does not implement "
"(only its molecular code does).\n"
" Use a functional without the -V suffix, e.g. "
"wb97x instead of wb97x-v.")
def _point_charge_potential(self, cell, column):
"""Grid values of -sum_i q_i/|r - R_i| for the removed sublattice.
Jellium-referenced (the G = 0 term is dropped), which is the same zero
of potential PySCF uses for the electrons and the nuclei, so all three
calculations share one energy reference. Validated against
:func:`crystod.point_charge_field.ewald_site_potential`.
"""
from pyscf.pbc import tools
other = {"left": "right", "right": "left"}[column]
sites = self.side_atoms[other]
if not sites:
return None
positions = np.array([self.cartesian[i] for i in sites]) / BOHR_TO_ANGSTROM
charges = np.array([self.oxidation[self.symbols[i]] for i in sites], dtype=float)
if not np.any(charges):
# all-zero oxidation states (neutral sublattices): no field, and
# no reason to pay the grid quadrature for exact zeros
return None
mesh = cell.mesh
Gv = cell.get_Gv(mesh)
coulG = tools.get_coulG(cell, mesh=mesh, Gv=Gv)
structure = np.exp(-1j * np.einsum("gx,ix->gi", Gv, positions))
# the electron charge is -1, hence the minus sign on the charges
potential_G = (structure @ (-charges)) * coulG
return tools.ifft(potential_G, mesh).real
def _point_charge_matrix(self, cell, column, kpts):
"""<phi_mu k| V_pointcharge |phi_nu k> for every k of ``kpts``.
The potential itself depends only on the cell and is cached; only the
AO quadrature is repeated, so this is cheap to call again for the
band k points of the diagram.
"""
from pyscf.pbc.dft import numint
if column not in self._pc_potential:
self._pc_potential[column] = self._point_charge_potential(cell, column)
potential = self._pc_potential[column]
if potential is None:
return None
coords = cell.get_uniform_grids(cell.mesh)
weight = cell.vol / len(coords)
matrices = []
for ao in numint.KNumInt().eval_ao(cell, coords, np.asarray(kpts).reshape(-1, 3)):
matrices.append(np.einsum("gi,g,gj->ij", ao.conj(), potential * weight, ao))
return np.asarray(matrices)
def _hcore_with_point_charges(self, cell, column, bare_get_hcore):
"""``get_hcore`` replacement that follows whatever k points it is asked
for -- the SCF mesh during the SCF, the special point during the band
step -- instead of freezing the mesh-sized matrix."""
cache: dict[bytes, np.ndarray] = {}
def get_hcore(cell_arg=None, kpts=None):
target = cell_arg if cell_arg is not None else cell
bare = np.asarray(bare_get_hcore(target, kpts))
if kpts is None:
requested = target.make_kpts([1, 1, 1])
elif hasattr(kpts, "kpts_ibz"):
# symmetry-adapted SCF: the Hamiltonian lives on the wedge
requested = np.asarray(kpts.kpts_ibz)
else:
requested = np.asarray(kpts)
requested = requested.reshape(-1, 3)
key = np.ascontiguousarray(requested).tobytes()
if key not in cache:
cache[key] = self._point_charge_matrix(target, column, requested)
correction = cache[key]
return bare if correction is None else bare + correction
return get_hcore
# ------------------------------------------------------------------ SCF
def _make_mean_field(self, cell, kpts, column, sigma=0.0, level_shift=0.0):
"""KRKS/KRHF with the point-charge field of the removed sublattice."""
from pyscf.pbc import dft, scf
if self.xc.lower() in {"hf", "hartree-fock"}:
mean_field = scf.KRHF(cell, kpts)
hybrid = True
else:
mean_field = dft.KRKS(cell, kpts)
mean_field.xc = self.xc
hybrid = abs(mean_field._numint.hybrid_coeff(self.xc)) > 1e-10
if hybrid and not getattr(self, "_hybrid_warned", False):
self._hybrid_warned = True
print("NOTE: exact exchange is evaluated on the FFT grid here; its grid error\n"
" is much larger than the semilocal one (raw degeneracy splittings\n"
" of a few 0.1 eV at ke_cutoff 80). The diagram re-diagonalizes the\n"
" group-averaged Fock, so the levels stay exactly symmetric, but for\n"
" hybrid energetics prefer --ke-cutoff 150 or more.")
mean_field.conv_tol = self.conv_tol
mean_field.max_cycle = self.max_cycle
mean_field.max_memory = self.max_memory
# MINAO/atom guesses assume all-electron shell structures and crash for
# compact GTH sets; hcore is the safe choice here
mean_field.init_guess = "hcore"
if column != "mo":
mean_field.get_hcore = self._hcore_with_point_charges(
cell, column, type(mean_field).get_hcore.__get__(mean_field)
)
if level_shift > 0:
# classic remedy for oscillating charged cells: bias the virtuals
# during the SCF (removed again before any eigenvalue is reported)
mean_field.level_shift = level_shift
if sigma > 0:
from pyscf.pbc.scf import addons
mean_field = addons.smearing_(mean_field, sigma=sigma / HARTREE_TO_EV,
method="fermi")
return mean_field
def _convergence_ladder(self):
"""Fallback (smearing eV, max cycles, level shift Hartree), in order."""
cycles = 3 * self.max_cycle
base = self.sigma if self.sigma > 0 else self.retry_sigma
return [
(base, cycles, 0.0),
# diffuse cation shells in a highly charged fragment cell make the
# aufbau occupations oscillate; a virtual-level shift stabilizes it
(base, cycles, 0.3),
(5.0 * base, cycles, 0.3),
]
[docs]
def run(self, report=print) -> None:
"""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 :attr:`mean_field`,
:attr:`density_matrix` and :attr:`scf_energy`, and are written to
:attr:`chk_path` when it is set. With ``onsite`` only the crystal
calculation runs.
Args:
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.
"""
import os
if self.chk_path and os.path.exists(self.chk_path):
if self._load_chk(report):
return
# The crystal is solved first: its converged density restricted to one
# sublattice's AO block (ghost rows/columns zeroed) is the best
# available guess for that fragment -- it differs from the pre-bonding
# state only by the bonding redistribution. The hcore guess that
# replaces it for the crystal itself is fine there, but for a highly
# charged cation fragment with diffuse shells (SrTi^6+ with dzvp) it
# starts so far away that the SCF never finds its way back.
slices = self.cells["mo"].aoslice_by_atom()
for column in self.scf_columns:
cell = self.cells[column]
# symmetry-reduced SCF mesh: only the irreducible wedge is solved
# (2x2x2 in Pm-3m: 4 of 8 k points), then to_khf() expands the
# result back to the full mesh for the band step and the guesses
try:
kpts = cell.make_kpts(self.kmesh, space_group_symmetry=True,
time_reversal_symmetry=True)
reduced = kpts.nkpts_ibz < kpts.nkpts
except Exception:
kpts = cell.make_kpts(self.kmesh)
reduced = False
if not reduced and not isinstance(kpts, np.ndarray):
kpts = cell.make_kpts(self.kmesh)
sigma = self.sigma
if column in self.odd_electron and sigma <= 0:
# an odd per-cell count needs fractional occupations (see
# _make_cell); Fermi smearing conserves it exactly
sigma = self.retry_sigma
report(f" {column:<5} has an odd electron count "
f"({cell.nelectron}): occupations use Fermi smearing, "
f"sigma {sigma:g} eV (set --sigma to change)")
self.smeared.add(column)
mean_field = self._make_mean_field(cell, kpts, column, sigma=sigma)
guess = None
if column != "mo" and "mo" in self.density_matrix:
if self.no_ghost:
# the fragment basis IS the sublattice block: take the
# crystal density restricted to it
rows = self.embed_rows[column]
guess = np.asarray(self.density_matrix["mo"])[
:, rows[:, None], rows[None, :]]
else:
mask = np.zeros(self.n_ao)
for atom in self.side_atoms[column]:
mask[int(slices[atom, 2]):int(slices[atom, 3])] = 1.0
guess = (np.asarray(self.density_matrix["mo"])
* mask[None, :, None] * mask[None, None, :])
if reduced:
# the stored crystal density lives on the full mesh; the
# symmetry-adapted SCF wants it on the irreducible wedge
guess = guess[np.asarray(kpts.ibz2bz)]
energy = mean_field.kernel(dm0=guess)
# A degenerate level straddling the Fermi energy makes the aufbau
# occupation jump between iterations (typical on a coarse mesh and
# for the highly charged cation fragments); fractional occupations
# and more cycles fix it. Escalate, and say which rung was used --
# smearing changes the level occupancies in the diagram.
for sigma, cycles, shift in self._convergence_ladder():
if mean_field.converged:
break
report(f" {column:<5} not converged; retrying with "
f"{sigma} eV Fermi smearing"
+ (f" and a {shift} Hartree virtual-level shift" if shift else "")
+ f" ({cycles} cycles)")
mean_field = self._make_mean_field(cell, kpts, column,
sigma=sigma, level_shift=shift)
mean_field.max_cycle = cycles
energy = mean_field.kernel(dm0=guess)
if sigma > 0:
self.smeared.add(column)
# the shift must not appear in any reported eigenvalue: get_bands
# rebuilds the Fock operator from the converged density, so with the
# attribute cleared the band energies are shift-free
mean_field.level_shift = 0.0
if not mean_field.converged:
raise SystemExit(
f"ERROR: the {column} SCF did not converge, even with Fermi "
"smearing. Try a different --kmesh, a larger --max-cycle or "
"--sigma, or a different --basis."
)
if reduced:
# expand the irreducible wedge back to the full mesh, exactly
# like the VASP-style band scripts; everything downstream
# (get_bands, the masked guesses) then sees a plain KRKS.
# to_khf() builds a fresh object, so the point-charge hcore
# override must be re-attached for the fragments.
mean_field = mean_field.to_khf()
if column != "mo":
mean_field.get_hcore = self._hcore_with_point_charges(
cell, column,
type(mean_field).get_hcore.__get__(mean_field),
)
self.mean_field[column] = mean_field
self.density_matrix[column] = mean_field.make_rdm1()
self.scf_energy[column] = float(energy)
report(f" {column:<5} {self.formula.get(column, 'crystal'):<8} "
f"E = {energy:16.8f} Hartree ({cell.nelectron} electrons, "
f"charge {cell.charge:+d})")
if self.chk_path:
self._save_chk(report)
# WAVECAR-style restart: the converged density matrices are all that the
# band step needs (get_bands rebuilds the Fock operator from them), so
# saving the three of them plus the defining parameters skips the SCFs
# entirely on the next run with the same --chk file.
def _chk_params(self) -> dict:
return {
"version": 1,
# densities written before the per-cell electron-count pin (see
# _make_cell) are wrong for charged fragments on a k-mesh; the
# loader waves the mismatch through when the bug could not bite
"nelectron_fix": 1,
"basis": self.basis_name, "pseudo": self.pseudo_name,
"xc": self.xc, "kmesh": list(self.kmesh),
"ke_cutoff": float(self.ke_cutoff), "max_l": self.max_l,
"no_ghost": self.no_ghost, "sigma": self.sigma,
# NOT compared: self.electrons only fills the diagram's arrows
# (_fill), so a different --electrons must not invalidate a
# perfectly good density. conv_tol/max_cycle DO shape it.
"conv_tol": float(self.conv_tol),
"max_cycle": int(self.max_cycle),
"oxidation": {el: float(q) for el, q in self.oxidation.items()},
"left": self.formula["left"], "right": self.formula["right"],
"symbols": list(self.symbols),
}
def _save_chk(self, report) -> None:
"""Write the checkpoint. Never fatal for an automatic file: the
cache is an optimization, and the three SCFs are already done --
losing the page because the working directory is read-only or full
would be absurd. Written to a temporary file and renamed, so an
interrupted write cannot leave a truncated checkpoint behind."""
import os
try:
self._write_chk()
except Exception as error:
if not self.chk_auto:
raise
report(f" WARNING: could not save {self.chk_path} ({error}); "
"continuing without the checkpoint")
for leftover in (self.chk_path + ".tmp",):
try:
os.remove(leftover)
except OSError:
pass
return
report(f" SCF saved to {self.chk_path} "
+ ("(reused automatically by the next --pyscf run on this "
"structure; delete the file to force a fresh SCF, or pass "
"--no-chk to skip it)"
if self.chk_auto else
"(reuse with --chk; delete the file to force a fresh "
"SCF)"))
def _write_chk(self) -> None:
import json
import os
temporary = self.chk_path + ".tmp"
with open(temporary, "wb") as handle:
# a file handle keeps the exact name (np.savez would append .npz)
# --onsite runs (and therefore saves) only the crystal SCF; the
# energy_columns array records which densities the file holds
np.savez_compressed(
handle,
params=json.dumps(self._chk_params(), sort_keys=True),
positions=self.positions, lattice=self.lattice,
smeared=np.array(sorted(self.smeared)),
energy_columns=np.array(list(self.scf_columns)),
energies=np.array([self.scf_energy[c]
for c in self.scf_columns]),
**{f"dm_{column}": np.asarray(self.density_matrix[column])
for column in self.scf_columns},
)
os.replace(temporary, self.chk_path)
def _chk_reject(self, reason, remedy, report) -> bool:
"""A checkpoint that cannot be reused.
A file the user named with --chk is an error -- they asked for that
density. The automatic CHK_{formula}.chk is a cache: say what
happened and recompute (the fresh SCF overwrites it)."""
if not self.chk_auto:
raise SystemExit(f"ERROR: {self.chk_path} {reason}\n"
f" {remedy}")
report(f" the automatic checkpoint {self.chk_path} {reason};")
report(" running the SCFs again and overwriting it")
return False
def _load_chk(self, report) -> bool:
"""True when the stored densities were adopted, False when an
automatic checkpoint was rejected (run() then recomputes)."""
import json
try:
data = np.load(self.chk_path, allow_pickle=False)
saved = json.loads(str(data["params"]))
except Exception as error:
# a truncated write, a foreign file that happens to match the
# automatic name, anything unreadable -- never a dead end
return self._chk_reject(
f"is not a readable CrystOD checkpoint ({error})",
"Delete the file (or point --chk somewhere else) and rerun.",
report)
current = self._chk_params()
# --onsite reads only the crystal density, which no_ghost never
# touches (it shapes the fragment SCFs) -- a full-run chk written
# with either setting is equally valid here
ignored = {"no_ghost"} if self.onsite else set()
mismatched = [key for key in current
if key not in ignored and saved.get(key) != current[key]]
if "nelectron_fix" in mismatched:
# checkpoint written before the per-cell electron-count pin (see
# _make_cell). Its stored densities are wrong exactly where the
# bug bit: a CHARGED fragment cell on a k-mesh larger than
# 1x1x1 converged with charge*(nkpts-1)/nkpts extra electrons
# per cell. Crystal-only (--onsite) and neutral-fragment
# checkpoints are unaffected and stay valid.
harmed = (int(np.prod(self.kmesh)) > 1
and any(int(round(self.side_charge[col])) != 0
for col in self.scf_columns if col != "mo"))
if harmed:
return self._chk_reject(
"predates the fragment electron-count fix: its charged "
"fragment densities converged with "
"charge*(nkpts-1)/nkpts spurious electrons per cell "
"(pyscf subtracts cell.charge once per "
"Born-von-Karman supercell, not per cell)",
"Delete the file and rerun to regenerate the densities.",
report)
mismatched.remove("nelectron_fix")
# shapes first: np.allclose RAISES on differently shaped arrays, and
# two polymorphs share one automatic file name (rocksalt and wurtzite
# AlN are both CHK_AlN.chk with 2 and 4 atoms per cell)
stored_positions = np.asarray(data["positions"])
stored_lattice = np.asarray(data["lattice"])
if (stored_positions.shape != self.positions.shape
or stored_lattice.shape != self.lattice.shape
or not np.allclose(stored_positions, self.positions,
atol=1e-6)
or not np.allclose(stored_lattice, self.lattice, atol=1e-6)):
mismatched.append("structure")
if mismatched:
return self._chk_reject(
"was written with different parameters "
f"({', '.join(sorted(mismatched))})",
"Delete the file or rerun with the matching options "
f"(crystod --chk-info {self.chk_path} shows the stored "
"conditions and a ready-to-paste option string).",
report)
self.smeared = set(str(s) for s in data["smeared"])
energies = np.asarray(data["energies"], dtype=float)
# pre-onsite checkpoints carry no energy_columns record; they always
# hold all three calculations in this fixed order
stored = ([str(c) for c in data["energy_columns"]]
if "energy_columns" in data.files
else ["mo", "left", "right"])
for column in self.scf_columns:
if f"dm_{column}" not in data.files:
return self._chk_reject(
"was written with --onsite and holds only the crystal "
f"density, not the {column} fragment column",
"Rerun without --chk (or with a checkpoint from a full "
"three-SCF run) to build the fragment columns.",
report)
cell = self.cells[column]
# a mean-field object is still needed for get_bands, but with the
# density matrix passed explicitly it is never iterated
kpts = cell.make_kpts(self.kmesh)
self.mean_field[column] = self._make_mean_field(
cell, kpts, column, sigma=self.sigma)
self.density_matrix[column] = np.asarray(data[f"dm_{column}"])
energy = float(energies[stored.index(column)])
self.scf_energy[column] = energy
report(f" {column:<5} {self.formula.get(column, 'crystal'):<8} "
f"E = {energy:16.8f} Hartree "
f"({cell.nelectron} electrons, charge {cell.charge:+d}) "
f"[read from {self.chk_path}]")
# only name the calculations this run actually uses (--onsite reads
# the crystal density alone; a smeared fragment SCF never enters it)
relevant = self.smeared & set(self.scf_columns)
if relevant:
report(f" ({', '.join(sorted(relevant))} carried Fermi "
"smearing when the file was written)")
return True
[docs]
def prepare_bands(self, kpoints) -> None:
"""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
:meth:`run`; :meth:`solve_at` then reads the cache.
Args:
kpoints: List of three-component primitive reciprocal
coordinates.
Returns:
``None``.
"""
self._band_cache = {}
keys = [self._kpoint_key(kpoint) for kpoint in kpoints]
for column in self.scf_columns:
cell = self.cells[column]
kpts_band = cell.get_abs_kpts(np.array(kpoints, dtype=float))
energies, coefficients = self.mean_field[column].get_bands(
kpts_band, cell=cell, dm_kpts=self.density_matrix[column]
)
for key, energy, coefficient in zip(keys, energies, coefficients):
self._band_cache[(column, key)] = (
np.asarray(energy).real * HARTREE_TO_EV,
self._embed(column, np.asarray(coefficient)),
)
@staticmethod
def _kpoint_key(kpoint):
return tuple(np.round(np.asarray(kpoint, dtype=float), 9))
def _bands_at(self, column, kpoint):
"""(energies in eV, coefficients) of one calculation at one k point."""
key = (column, self._kpoint_key(kpoint))
cached = getattr(self, "_band_cache", {}).get(key)
if cached is not None:
return cached
cell = self.cells[column]
kpts_band = cell.get_abs_kpts(np.array([kpoint], dtype=float))
energies, coefficients = self.mean_field[column].get_bands(
kpts_band, cell=cell, dm_kpts=self.density_matrix[column]
)
return (np.asarray(energies[0]).real * HARTREE_TO_EV,
self._embed(column, np.asarray(coefficients[0])))
def _population_shares(self, vectors, degeneracy, S, S_half, specs=None):
"""Sorted per-(element, shell) AO populations of one level space in
the selected projection (the displayed composition measure).
``specs`` restricts the listed shells (fragment columns list their
own sublattice only)."""
if self.projection == "mulliken":
gross = (vectors.conj() * (S @ vectors)
).real.sum(axis=1) / degeneracy
else:
gross = (np.abs(S_half @ vectors) ** 2).sum(axis=1) / degeneracy
shares = []
for spec in (self.specs if specs is None else specs):
value = float(gross[self.spec_indices[id(spec)]].sum())
if abs(value) >= 0.001:
shares.append((value, f"{spec.element} {spec.shell}"))
shares.sort(key=lambda item: -item[0])
return shares
def _block_bands(self, column, S):
"""Per-shell on-site multiplets of the crystal Fock (--onsite).
For every (element, shell) of the sublattice, diagonalizes (F, S)
restricted to THAT SHELL's own Bloch AOs -- the tight-binding
on-site energies <phi|F|phi> of the actual shells, one level per
induced irrep, matching the site-symmetry table exactly. No
cross-shell mixing on purpose: diagonalizing a whole sublattice
block instead is variationally unstable with this diffuse basis --
the raw AO block lets diffuse cation functions fall into the
removed side's potential wells (a "Ti 3d GM3+" at O-2p depth with
28% Loewdin weight on O: a disguised anion state, caught by the
user on SrTiO3), while Loewdin- or projection-orthogonalized
blocks load strongly overlapped shells with huge orthogonalization
penalties (the O 2s on-site swung from -6.9 to +0.5 to +12.6 eV
across those three conventions; per-shell Rayleigh quotients have
no such freedom). Shell AO sets map onto themselves under the
space group, so exact multiplets of the group-averaged Fock carry
over (raw FFT-grid splittings under --no-symmetrize).
"""
fock = self.fock["mo"]
chunks = []
dropped = 0
for spec in self.side_specs[column]:
rows = np.asarray(self.spec_indices[id(spec)])
overlap = S[np.ix_(rows, rows)]
values, vectors = np.linalg.eigh(overlap)
keep = values > ONSITE_OVERLAP_FLOOR
dropped += int(values.size - int(keep.sum()))
basis = vectors[:, keep] / np.sqrt(values[keep])
energies, mixing = np.linalg.eigh(
basis.conj().T @ fock[np.ix_(rows, rows)] @ basis)
states = basis @ mixing
full = np.zeros((self.n_ao, states.shape[1]), dtype=complex)
full[rows] = states
chunks.append((energies.real, full))
energies = np.concatenate([e for e, _ in chunks])
coefficients = np.hstack([c for _, c in chunks])
order = np.argsort(energies)
return energies[order], coefficients[:, order], dropped
[docs]
def overlap_at(self, kpoint) -> np.ndarray:
"""PySCF AO overlap matrix ``S(k)`` of the crystal cell.
Args:
kpoint: Three primitive reciprocal coordinates.
Returns:
Hermitian complex array of shape ``(n_ao, n_ao)``.
"""
cell = self.cells["mo"]
kpts_band = cell.get_abs_kpts(np.array([kpoint], dtype=float))
return np.asarray(cell.pbc_intor("int1e_ovlp", hermi=1, kpts=kpts_band))[0]
@staticmethod
def _sqrt_overlap(overlap) -> np.ndarray:
"""S^(1/2), for Loewdin populations (|coefficient|^2 in the
symmetrically orthogonalized basis -- the orthonormal set closest to
the atomic orbitals, i.e. the site-bound attribution)."""
values, vectors = np.linalg.eigh(overlap)
return (vectors * np.sqrt(np.clip(values, 0.0, None))) @ vectors.conj().T
# -------------------------------------------------------------- symmetry
[docs]
def little_group_data(self, kpoint):
"""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]``.
Args:
kpoint: Three primitive reciprocal coordinates.
Returns:
``(irreps, mapping, labels, representation)`` as in
:meth:`CrystalOrbitalDiagram.little_group_data`.
"""
from spgrep.core import get_spacegroup_irreps_from_primitive_symmetry
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)
permutations = self.builder.get_permutation_reps_at_k(
little_rotations=self.builder.rotations[mapping],
little_translations=self.builder.translations[mapping],
kpoint=kpoint,
)
# crystod builds the permutation representation in the *periodic* gauge,
# chi_mu k = sum_T exp(i k (T + tau_mu)) phi_mu(r - tau_mu - T), while
# PySCF puts the phase on the lattice translation alone (the atomic
# gauge). The two bases differ by Lambda_mu = exp(2 pi i k . tau_mu),
# so the representation must be conjugated with it. The difference is
# invisible whenever k . tau is a multiple of 1/2 for every site (ScF3,
# SrTiO3) and essential when it is not (the 1/4, 3/4 sites of fluorite).
slices = self.cells["mo"].aoslice_by_atom()
site_phase = np.ones(self.n_ao, dtype=complex)
for atom in range(len(self.symbols)):
value = np.exp(2j * np.pi * float(np.dot(kpoint, self.positions[atom])))
site_phase[int(slices[atom, 2]):int(slices[atom, 3])] = value
blocks = self.ao_blocks
representation = []
for op_index, op in enumerate(mapping):
rotation = np.real(self.builder.rotations_cartesian[op])
permutation = permutations[op_index]
wigners: dict[int, np.ndarray] = {}
matrix = np.zeros((self.n_ao, self.n_ao), dtype=complex)
for block in blocks:
if block.l not in wigners:
wigners[block.l] = wigner_pyscf(block.l, rotation)
wigner = wigners[block.l]
atom = block.sites[0]
for image in blocks:
if image.l != block.l:
continue
if self.ao_offset_in_atom[image.offset] != self.ao_offset_in_atom[block.offset]:
continue
phase = permutation[image.sites[0], atom]
if abs(phase) < 1e-12:
continue
matrix[image.offset:image.offset + image.n_ao,
block.offset:block.offset + block.n_ao] = phase * wigner
representation.append(site_phase[:, None] * matrix / site_phase[None, :])
return irreps, mapping, labels, representation
[docs]
def site_symmetry_irreps(self, kpoint, representation, irreps, labels):
"""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.
Args:
kpoint: Three primitive reciprocal coordinates (not used by the
computation; kept for symmetry with :meth:`little_group_data`).
representation: The AO representation matrices from
:meth:`little_group_data`.
irreps: The spgrep irreps from the same call.
labels: Their ISO-IR labels.
Returns:
``{(element, shell): ["GM1+", "2GM4-", ...]}``: the irreps the
shell's Bloch orbitals span, multiplicities as prefixes.
"""
from .runtime_compat import get_character
order = len(representation)
content = {}
for spec in self.specs:
indices = self.spec_indices[id(spec)]
characters = np.array([
np.trace(D[np.ix_(indices, indices)]) for D in representation
])
parts = []
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:
name = label.split("(")[0]
parts.append(name if count == 1 else f"{count}{name}")
content[(spec.element, spec.shell)] = parts
return content
# --------------------------------------------------------------- solving
def _group_levels(self, energies, vectors):
"""Cluster eigenvalues into degenerate groups (eV, grid-noise aware)."""
groups = []
start = 0
for index in range(1, len(energies) + 1):
if (index == len(energies)
or energies[index] - energies[start] > self.degeneracy_tol):
groups.append((float(np.mean(energies[start:index])),
vectors[:, start:index]))
start = index
return groups
def _multiplicities(self, space, overlap, representation, irreps):
"""Irrep content of an eigenspace (may be fractional if it is not a
whole number of complete multiplets)."""
from .runtime_compat import get_character
gram = space.conj().T @ overlap @ space
values, basis = np.linalg.eigh(gram)
keep = values > 1e-10
orthonormal = space @ (basis[:, keep] / np.sqrt(values[keep]))
characters = np.array([
np.trace(orthonormal.conj().T @ overlap @ D @ orthonormal)
for D in representation
])
order = len(representation)
return [
float(np.real(np.sum(characters * np.conj(
np.array(get_character(irrep), dtype=complex))) / order))
for irrep in irreps
]
def _is_complete_multiplet(self, space, overlap, representation, irreps) -> bool:
multiplicities = self._multiplicities(space, overlap, representation, irreps)
if not any(value > 0.5 for value in multiplicities):
return False
return all(abs(value - round(value)) < _MULTIPLICITY_TOL
for value in multiplicities)
def _adaptive_groups(self, energies, vectors, overlap, representation, irreps):
"""Cluster levels into complete multiplets.
The seed window only separates obviously distinct levels; adjacent
seeds are then merged until the irrep multiplicities of the group are
integral, which is exactly the condition that no multiplet is cut in
half. Grid noise splits symmetry-degenerate levels by anything from
1e-4 to several 1e-2 eV depending on the system and the cutoff, so this
is far more robust than any fixed tolerance -- and the merging can
never run away, because it stops at DEGENERACY_MAX_WINDOW_EV.
"""
seeds = self._group_levels(energies, vectors)
groups = []
index = 0
while index < len(seeds):
start_energy = seeds[index][0]
space = seeds[index][1]
weights = [space.shape[1]]
centres = [seeds[index][0]]
last = index
while not self._is_complete_multiplet(space, overlap, representation, irreps):
if (last + 1 >= len(seeds)
or seeds[last + 1][0] - start_energy > self.degeneracy_window):
break
last += 1
space = np.hstack([space, seeds[last][1]])
weights.append(seeds[last][1].shape[1])
centres.append(seeds[last][0])
energy = float(np.average(centres, weights=weights))
groups.append((energy, space))
index = last + 1
return groups
[docs]
def align_fragment_columns(self, records):
"""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.
Args:
records: The per-k-point list built by :func:`report_and_write`:
dicts with ``"name"``, ``"kpoint"`` and ``"levels"`` (the
:meth:`solve_at` result), one per diagram k point. The
level energies are shifted in place.
Returns:
``(shifts, anchors)``: ``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.
"""
anchors = {}
for column in ("left", "right"):
pairs = [] # (fragment_energy, delta, fragment_label, k_name, purity)
for record in records:
fragment_levels = {
level.level_id: level for level in record["levels"][column]
}
best: dict[str, tuple[float, object]] = {}
for crystal in record["levels"]["mo"]:
# absolute projections: the renormalized composition can
# show ~100% for a level whose true overlap with the
# retained fragment states is tiny
weights = getattr(crystal, "absolute_composition",
crystal.composition)
for level_id, weight in weights:
if not level_id.startswith(column):
continue
if weight > best.get(level_id, (0.0, None))[0]:
best[level_id] = (weight, crystal)
for level_id, (purity, crystal) in best.items():
fragment = fragment_levels[level_id]
pairs.append((fragment.energy, crystal.energy - fragment.energy,
fragment.label, record["name"], purity))
inert = [pair for pair in pairs if pair[4] >= ALIGNMENT_PURITY]
pool = inert or pairs
if not pool:
anchors[column] = None
continue
deepest = min(pool, key=lambda pair: pair[0])
key = tuple(deepest[2].split()[:2]) # (element, shell), any k
# the same shell across the k points -- but only levels of the same
# band (within 1 eV of the anchor), since one shell can span several
# irreps of quite different chemistry at a single k point
same = [pair for pair in pool
if tuple(pair[2].split()[:2]) == key
and abs(pair[0] - deepest[0]) < 1.0]
# the reference offset is one constant, but the anchor shell may
# hybridize more at some k than at others; trust only the pairs
# within 2% of the best purity, where the chemistry is smallest
best_purity = max(pair[4] for pair in same)
trusted = [pair for pair in same if pair[4] >= best_purity - 0.02]
values = [pair[1] for pair in trusted]
purest = max(trusted, key=lambda pair: pair[4])
anchors[column] = {
"label": purest[2],
"fragment_energy": purest[0],
"delta": float(np.mean(values)),
"spread": float(np.max(values) - np.min(values)) if len(values) > 1 else 0.0,
"n_k": len(trusted),
"purity": best_purity,
"fallback": not inert,
}
available = [c for c in ("left", "right") if anchors.get(c)]
if not available:
return None, anchors
reference = min(available, key=lambda c: anchors[c]["fragment_energy"])
delta_ref = anchors[reference]["delta"]
shifts = {reference: 0.0, "mo": -delta_ref}
other = {"left": "right", "right": "left"}[reference]
if anchors.get(other):
shifts[other] = anchors[other]["delta"] - delta_ref
else:
shifts[other] = -delta_ref # no anchor: move with the crystal
shifts["reference"] = reference
for record in records:
for column in ("left", "right", "mo"):
for level in record["levels"][column]:
level.energy += shifts[column]
return shifts, anchors
# ------------------------------------------- isolated formal-charge ions
[docs]
def atomic_ion_levels(self):
"""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:
``None``. Fills :attr:`atomic_ions` with one entry per element:
``{"charge", "nelec", "method", "shells": [(shell_name, l,
energy_eV), ...], "caveats"}``, one shell per (element, shell)
spec of the AO basis.
"""
from pyscf import gto as mol_gto
from pyscf import dft as mol_dft
from pyscf import scf as mol_scf
from pyscf.data.nist import HARTREE2EV
self.atomic_ions = {}
element_shells: dict[str, list] = {}
for column in ("left", "right"):
for spec in self.side_specs[column]:
element_shells.setdefault(spec.element, []).append(
(spec.shell, spec.l))
for element, shells in element_shells.items():
charge = int(round(self.oxidation[element]))
basis_shells = mol_gto.basis.load(self.basis_name, element)
if self.max_l is not None:
basis_shells = [shell for shell in basis_shells
if shell[0] <= self.max_l]
def build_ion(spin):
return mol_gto.M(
atom=f"{element} 0.0 0.0 0.0",
basis={element: basis_shells},
pseudo={element: self.pseudo_name},
charge=charge, spin=spin,
verbose=0, max_memory=self.max_memory,
)
try:
mol = build_ion(0)
except RuntimeError:
# odd electron count: parity forces an odd spin
mol = build_ion(1)
nelec = mol.nelectron
if nelec == 0:
# no electrons left: the one-electron spectrum of the bare
# pseudo-ion (canonically orthogonalized hcore)
hcore = mol_scf.RHF(mol).get_hcore()
overlap = mol.intor("int1e_ovlp")
values, vectors = np.linalg.eigh(overlap)
keep = values > 1e-10
X = vectors[:, keep] / np.sqrt(values[keep])
energies, rotation = np.linalg.eigh(X.conj().T @ hcore @ X)
coeff = X @ rotation
occupations = np.zeros(len(energies))
method = "0 electrons: bare-ion one-electron levels (hcore)"
else:
# ground spin state by trial, not by assumption: scan the
# parity-consistent spins (up to 8 unpaired electrons --
# Gd/Cm reach 4f7 5d1 / 5f7 6d1) and keep the lowest
# converged energy. A neutral N atom (--oxidation N=0)
# takes spin 3 (Hund's 4S term), not the spin-1 doublet a
# bare parity fallback would pick; closed shells still win
# at spin 0. Atomic calculations cost well under a second.
def solve_ion(spin):
trial = mol if spin == mol.spin else build_ion(spin)
mean_field = (mol_dft.RKS(trial) if spin == 0
else mol_dft.UKS(trial))
mean_field.xc = self.xc
mean_field.conv_tol = self.conv_tol
mean_field.max_cycle = self.max_cycle
mean_field.kernel()
if not mean_field.converged:
# DIIS stalls on near-degenerate open d shells
# (neutral Ni, spin 0/2); second-order SCF usually
# lands them. Without this the scan would settle
# on a converged EXCITED spin state and say nothing.
mean_field = mean_field.newton()
mean_field.kernel()
return mean_field
solutions = []
for trial_spin in range(mol.spin, min(nelec, 8) + 1, 2):
try:
solutions.append(solve_ion(trial_spin))
except Exception as error:
# e.g. more alpha electrons than the (--max-l
# truncated) basis has orbitals: skip the trial,
# do not kill the diagram
print(f" note: the {element}^{charge:+d} spin-"
f"{trial_spin} trial failed ({error})")
continue
if not solutions:
raise SystemExit(
f"ERROR: no spin state of the isolated "
f"{element}^{charge:+d} ion could be solved in "
"this basis (is --max-l too small for its "
"electron count?)")
converged = [sol for sol in solutions if sol.converged]
pick = min(converged or solutions,
key=lambda sol: float(sol.e_tot))
lowest = min(solutions, key=lambda sol: float(sol.e_tot))
if float(lowest.e_tot) < float(pick.e_tot) - 1e-6:
print(f" WARNING: an unconverged spin-"
f"{lowest.mol.spin} state of {element}^"
f"{charge:+d} lies {(float(pick.e_tot) - float(lowest.e_tot)) * HARTREE2EV:.2f} eV "
f"below the chosen spin-{pick.mol.spin} solution "
"-- the atomic column may show an excited spin "
"state")
mol = pick.mol
spin = mol.spin
if spin == 0:
energies, coeff = pick.mo_energy, pick.mo_coeff
occupations = pick.mo_occ
method = f"RKS {self.xc.upper()}, {nelec} electrons"
else:
# open-shell ion or neutral atom: alpha channel,
# flagged in the tooltip
energies, coeff = (pick.mo_energy[0],
pick.mo_coeff[0])
occupations = pick.mo_occ[0]
method = (f"UKS {self.xc.upper()}, {nelec} electrons, "
f"spin {spin} (alpha levels)")
if not pick.converged:
method += "; WARNING: SCF not converged"
# dominant l of each MO (Loewdin), then per-l pairing of the MO
# multiplets (energy order) with the element's shells (2s < 3s)
overlap = mol.intor("int1e_ovlp")
values, vectors = np.linalg.eigh(overlap)
sqrt_overlap = (vectors * np.sqrt(np.clip(values, 0.0, None))
) @ vectors.conj().T
l_of_ao = np.repeat(
[mol.bas_angular(shell) for shell in range(mol.nbas)
for _ in range(mol.bas_nctr(shell))],
[2 * mol.bas_angular(shell) + 1 for shell in range(mol.nbas)
for _ in range(mol.bas_nctr(shell))])
gross = np.abs(sqrt_overlap @ coeff) ** 2
mo_l = np.array([
int(np.argmax([gross[l_of_ao == l, i].sum()
for l in range(int(l_of_ao.max()) + 1)]))
for i in range(coeff.shape[1])
])
shell_levels = []
# empty-shell caveats: for a CATION the only spectroscopically
# trustworthy empty level per l is the lowest EMPTY one --
# occupied semicore shells below it do not count (Na+ 3s, the
# first empty s above the occupied 2s, is basis-converged to
# 0.2 eV; validated on Al^3+ against NIST Al III: 3s/3p/3d
# fine, 4s off by +7 eV in gth-dzvp-molopt-sr, which carries
# no diffuse functions). Higher empty multiplets are
# finite-basis virtuals; for an anion EVERY empty level is a
# discretized continuum state (one more electron is
# vacuum-unbound at any basis size), and a formally NEUTRAL
# atom has no Coulomb tail to guarantee bound empty states,
# so all of its empty shells are flagged too.
caveats: dict[str, str] = {}
for l in sorted({l for _, l in shells}):
names = sorted((name for name, shell_l in shells
if shell_l == l),
key=lambda name: int(name[:-1]))
indices = np.where(mo_l == l)[0]
ordered = indices[np.argsort(energies[indices])]
degeneracy = 2 * l + 1
multiplets = [ordered[i:i + degeneracy]
for i in range(0, len(ordered), degeneracy)]
if len(multiplets) != len(names):
print(f" WARNING: {element}^{charge:+d}: "
f"{len(multiplets)} l={l} ion multiplets for "
f"{len(names)} shells; pairing the lowest ones")
empty_rank = 0
for name, multiplet in zip(names, multiplets):
if float(occupations[multiplet].sum()) < 1e-6:
if charge < 0:
caveats[name] = "continuum"
elif charge == 0 or empty_rank > 0:
caveats[name] = "virtual"
empty_rank += 1
spread = float(energies[multiplet].max()
- energies[multiplet].min()) * HARTREE2EV
if spread > 0.1:
# a torn multiplet: the energy-ordered chunking
# assumed exact degeneracy (symmetry-broken UKS
# solutions of open-shell ions can violate it)
print(f" WARNING: the {element}^{charge:+d} "
f"{name} ion multiplet is not degenerate "
f"(spread {spread:.2f} eV; symmetry-broken "
"open-shell solution?) -- its mean is shown")
shell_levels.append((
name, l,
float(np.mean(energies[multiplet])) * HARTREE2EV))
self.atomic_ions[element] = {
"charge": charge, "nelec": nelec, "method": method,
"shells": shell_levels, "caveats": caveats,
}
[docs]
def attach_atomic_columns(self, records):
"""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.
Args:
records: The per-k-point list of :meth:`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.
"""
site_counts: dict[str, set] = {}
for column in ("left", "right"):
for spec in self.side_specs[column]:
site_counts.setdefault(spec.element, set()).update(spec.sites)
# per-(element, shell) fragment band centers over all k points
sums: dict[tuple[str, str], float] = {}
weights: dict[tuple[str, str], float] = {}
for record in records:
for column in ("left", "right"):
for level in record["levels"][column]:
parts = level.label.split()
if len(parts) < 3:
continue
key = (parts[0], parts[1])
sums[key] = sums.get(key, 0.0) \
+ level.energy * level.degeneracy
weights[key] = weights.get(key, 0.0) + level.degeneracy
anchors: dict[str, tuple[str, float]] = {}
shifts: dict[str, float] = {}
for element, ion in self.atomic_ions.items():
anchor = None
for name, _, energy in sorted(ion["shells"],
key=lambda item: item[2]):
if (element, name) in weights:
anchor = (name, energy)
break
if anchor is None:
shifts[element] = 0.0
anchors[element] = ("none", 0.0)
continue
center = (sums[(element, anchor[0])]
/ weights[(element, anchor[0])])
shifts[element] = center - anchor[1]
anchors[element] = (anchor[0], shifts[element])
for record in records:
levels = record["levels"]
for column in ("left", "right"):
ao_column = f"{column}-ao"
levels[ao_column] = []
ao_ids: dict[tuple[str, str], str] = {}
elements = list(dict.fromkeys(
spec.element for spec in self.side_specs[column]))
entries = []
for element in elements:
ion = self.atomic_ions[element]
count = len(site_counts[element])
for name, l, energy in ion["shells"]:
entries.append((energy + shifts[element], element,
name, l, count, ion))
entries.sort(key=lambda item: item[0])
equivalent = self.builder.spglib_dataset["equivalent_atoms"]
for energy, element, name, l, count, ion in entries:
prefix = f"{count}" if count > 1 else ""
level = DiagramLevel(
level_id=f"{ao_column}{len(levels[ao_column])}",
column=ao_column,
energy=energy,
degeneracy=2 * l + 1,
irrep="",
label=f"{prefix}{element} {name}",
vectors=None,
)
level.electrons = None
level.display_composition = [(f"{element} {name}", 1.0)]
anchor_name, shift = anchors[element]
if anchor_name == "none":
anchor_line = (
"WARNING: no fragment level is dominated by a "
f"{element} shell -- drawn on the raw molecular "
"vacuum reference (no anchor)")
else:
anchor_line = (
f"vacuum level {energy - shift:+.2f} eV, "
f"shifted {shift:+.2f} eV so the deepest shell "
f"with a fragment counterpart ({element} "
f"{anchor_name}) sits at its "
f"{self.formula[column]}-column band center")
orbits = len({int(equivalent[site])
for site in site_counts[element]})
if count > 1:
inequivalent = ("" if orbits == 1 else
f" ({orbits} inequivalent positions)")
tail = (f"the {self.formula[column]} column shows "
f"how the {count} {element} "
f"ions'{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")
caveat = ion.get("caveats", {}).get(name, "")
if caveat == "virtual":
if ion["charge"] > 0:
opening = (
"empty shell beyond the lowest EMPTY state "
"of its l -- a finite-basis VIRTUAL, not a "
"physical Rydberg level")
else:
opening = (
"empty shell of a formally NEUTRAL atom -- "
"a finite-basis VIRTUAL (no ionic Coulomb "
"tail guarantees a bound counterpart)")
if self.richer_basis:
remedy = (
f" ({self.basis_name} has no diffuse "
"functions, so its energy and even its s/p "
"order follow the contraction; a "
"diffuse-richer basis such as --basis "
f"{self.richer_basis} removes the "
"basis-side error -- bare-ion Al^+3 4s/4p "
"land within 0.1 eV of NIST -- though "
"Kohn-Sham virtuals of ions that keep "
"electrons retain the functional's own "
"eV-scale attachment error)")
elif "molopt" in self.basis_name.lower():
remedy = (
f" ({self.basis_name} has no diffuse "
"functions, so its energy and even its s/p "
"order follow the contraction; PySCF ships "
"no diffuse-richer GTH basis for these "
"elements -- the molopt-sr sets are the "
"only ones reaching beyond Ar)")
else:
remedy = (
" (its position is limited by the basis's "
"diffuse coverage and, for ions that keep "
"electrons, by the functional's Kohn-Sham "
"virtual error)")
caveat = "\nNOTE: " + opening + remedy
elif caveat == "continuum":
caveat = (
"\nNOTE: one more electron on this anion is "
"vacuum-unbound, so this empty level is a "
"discretized CONTINUUM state of the finite "
"basis, with no physical counterpart at any "
"basis size; in the crystal it is the Madelung "
"potential that provides the binding")
species = "atom" if ion["charge"] == 0 else "ion"
level.detail = (
f"isolated {element}^{ion['charge']:+d} {species} "
f"({ion['method']}; PySCF {self.basis_name}, same "
"pseudopotential as the crystal)\n"
f"{anchor_line}\n{tail}{caveat}")
ao_ids[(element, name)] = level.level_id
levels[ao_column].append(level)
# splitting connector lines, weighted like the tooltip rows
for level in levels[column]:
linked = []
for name_irrep, value in getattr(
level, "display_composition", []):
parts = name_irrep.split()
if len(parts) >= 2 \
and (parts[0], parts[1]) in ao_ids \
and value >= 0.001:
linked.append(
(ao_ids[(parts[0], parts[1])], value))
level.composition = linked
return anchors
# ---------------------------------------------------- wave-function sketch
[docs]
def sketch_partners(self, level, kpoint, sites):
"""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
(:meth:`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.
Args:
level: A ``DiagramLevel`` from :meth:`solve_at`.
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 .visualize_basis import realify_basis_space
rows, _ = realify_basis_space(level.vectors.T)
rows = np.asarray(rows)
if self.projection == "mulliken":
overlap = getattr(level, "overlap", None)
if overlap is None:
overlap = self.overlap_at(kpoint)
else:
overlap = getattr(level, "sqrt_overlap", None)
if overlap is None:
overlap = self._sqrt_overlap(self.overlap_at(kpoint))
n_prim = len(self.symbols)
width = 9
slot_of = {0: 0, 1: 1, 2: 4}
angular = {0: 0.28209479, 1: 0.48860251, 2: 0.63078313}
# a fragment level is drawn from its own sublattice's components only:
# amplitude on the ghost basis of the removed sublattice is the
# variational tail toward the point charges (symmetry-allowed BSSE
# borrowing), and drawing it as an atom-centred orbital on the empty
# site misreads it; the ghost fraction is reported numerically instead
column = getattr(level, "column", "mo")
blocks = [b for b in self.ao_blocks
if column == "mo" or b.column == column]
# one representative display row per primitive atom (|amp| is the
# same on every translate of an atom, the Bloch phase is unimodular)
representative: dict[int, int] = {}
for row_index, (atom, _) in enumerate(sites):
representative.setdefault(atom, row_index)
raw_partners = []
channel_pop = np.zeros((n_prim, width)) # multiplet-summed Loewdin
channel_amp2 = np.zeros((n_prim, width)) # multiplet-summed |amp|^2
for vector in rows:
v = np.asarray(vector, dtype=complex)
# per-AO population of this partner in the selected projection
# (--projection): Loewdin |S^(1/2)c|^2 -- squared coefficients in
# the symmetrically orthogonalized basis, non-negative and summing
# to 1 -- or Mulliken gross populations Re[c* (S c)], whose
# overlap cross terms can go negative on diffuse empty levels
if self.projection == "mulliken":
gross = (v.conj() * (overlap @ v)).real
else:
gross = np.abs(overlap @ v) ** 2
amp_re = np.zeros((len(sites), width))
amp_im = np.zeros_like(amp_re)
for block in blocks:
if block.l not in slot_of:
continue
slot = slot_of[block.l]
site = block.sites[0]
channel_pop[site, slot] += float(
gross[block.offset:block.offset + block.n_ao].sum())
coefficients = np.asarray(
vector[block.offset:block.offset + block.n_ao]
) * (angular[block.l] * block.radial)
for row_index, (atom, translation) in enumerate(sites):
if atom != site:
continue
phase = np.exp(2j * np.pi * float(np.dot(kpoint, translation)))
values = coefficients * phase
amp_re[row_index, slot:slot + block.n_ao] += values.real
amp_im[row_index, slot:slot + block.n_ao] += values.imag
for site in range(n_prim):
row0 = representative[site]
for l, slot in slot_of.items():
n_m = 2 * l + 1
channel = (amp_re[row0, slot:slot + n_m]
+ 1j * amp_im[row0, slot:slot + n_m])
channel_amp2[site, slot] += float(np.linalg.norm(channel)) ** 2
raw_partners.append((amp_re, amp_im))
# one calibration factor per (atom, l) channel, from the MULTIPLET
# sums: within a degenerate multiplet the population/amplitude^2 ratio
# of a channel is partner-independent by symmetry, so a shared factor
# keeps symmetry-equivalent atoms exactly equal and preserves the
# sigma/pi contrast between partners, while lobe areas become
# proportional to the Loewdin electron weight on the site instead of
# the diffuse-function amplitude at r0 (a gth Sc 4p is ~5x an F 2p
# there; per-partner rescaling instead broke the F-site equivalence
# slightly through the RREF canonicalization below)
col_factor = np.zeros((n_prim, width))
for site in range(n_prim):
for l, slot in slot_of.items():
n_m = 2 * l + 1
amp2 = channel_amp2[site, slot]
if amp2 > 1e-24:
col_factor[site, slot:slot + n_m] = np.sqrt(
max(channel_pop[site, slot], 0.0) / amp2)
tiled = col_factor[[atom for atom, _ in sites], :]
partner_rows = []
for amp_re, amp_im in raw_partners:
amp_re = amp_re * tiled
amp_im = amp_im * tiled
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
def _symmetrize_fock(self, fock, overlap, representation):
"""Group-average an operator so its eigenspaces carry complete irreps
EXACTLY. The raw SCF has no point-group constraint -- the FFT grid
(and, much worse, grid-evaluated exact exchange of hybrid functionals)
breaks degeneracies by up to a few 0.1 eV -- so the displayed levels
are re-diagonalized from F_avg = (1/|G|) sum_g D(g)+ F D(g), which is
invariant under the verified AO representation by construction.
Returns (energies, coefficients, max |F_avg - F|)."""
average = sum(D.conj().T @ fock @ D for D in representation)
average = average / len(representation)
average = 0.5 * (average + average.conj().T)
deviation = float(np.max(np.abs(average - fock)))
s_values, s_vectors = np.linalg.eigh(overlap)
keep = s_values > 1e-9 * float(s_values.max())
basis = s_vectors[:, keep] / np.sqrt(s_values[keep])
reduced = basis.conj().T @ average @ basis
reduced = 0.5 * (reduced + reduced.conj().T)
energies, rotation = np.linalg.eigh(reduced)
return energies, basis @ rotation, deviation
def _symmetrized_bands(self, column, energies, coefficients, overlap,
representation):
"""Re-derive one calculation's levels from the group-averaged Fock.
With --no-ghost the fragment spans only its own AO block, so the
averaging runs in that subspace (the space-group operations never mix
the sublattices, hence D is block-diagonal in them).
"""
if self.no_ghost and column != "mo":
rows = self.embed_rows[column]
sub = np.ix_(rows, rows)
small_overlap = overlap[sub]
small_representation = [D[sub] for D in representation]
small_c = coefficients[rows]
fock = (small_overlap @ small_c @ np.diag(energies)
@ small_c.conj().T @ small_overlap)
e, c, deviation = self._symmetrize_fock(
fock, small_overlap, small_representation)
return e, self._embed(column, c), deviation
fock = (overlap @ coefficients @ np.diag(energies)
@ coefficients.conj().T @ overlap)
return self._symmetrize_fock(fock, overlap, representation)
def _ghost_fraction(self, column, space, overlap) -> float:
"""Mean Mulliken population of a level space on the ghost basis of the
*other* sublattice -- a counterpoise artifact when it dominates."""
other = {"left": "right", "right": "left"}[column]
slices = self.cells["mo"].aoslice_by_atom()
rows: list[int] = []
for atom in self.side_atoms[other]:
rows.extend(range(int(slices[atom, 2]), int(slices[atom, 3])))
if not rows:
return 0.0
rows = np.array(rows, dtype=int)
gross = (space.conj() * (overlap @ space)).real
total = np.clip(gross.sum(axis=0), 1e-12, None)
return float(np.mean(gross[rows, :].sum(axis=0) / total))
def _dominant_spec(self, space, S, column):
"""(element, shell) block with the largest Mulliken population."""
weighted = S @ space
best, best_weight = None, -np.inf
for spec in self.side_specs[column]:
indices = self.spec_indices[id(spec)]
weight = float(np.real(np.sum(np.conj(space[indices]) * weighted[indices])))
if weight > best_weight:
best, best_weight = spec, weight
return best
[docs]
def solve_at(self, kpoint):
"""Solve the fragment and crystal levels at one k point from the SCFs.
Reads the band energies and coefficients cached by
:meth:`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 :attr:`last_coupling`.
Args:
kpoint: Three primitive reciprocal coordinates.
Returns:
``(levels, labels)`` with the columns ``"left"``, ``"mo"`` and
``"right"`` as in :meth:`CrystalOrbitalDiagram.solve_at` (the
``-ao`` columns are added later by :meth:`attach_atomic_columns`);
energies in eV on the raw per-calculation references until
:meth:`align_fragment_columns` shifts them.
Raises:
SystemExit: The AO representation does not leave the PySCF
overlap invariant (please report the case).
"""
S = self.overlap_at(kpoint)
S_half = self._sqrt_overlap(S)
irreps, mapping, labels, representation = self.little_group_data(kpoint)
worst = max(float(np.max(np.abs(D.conj().T @ S @ D - S))) for D in representation)
if worst > 1e-6:
raise SystemExit(
"ERROR: the AO representation does not leave the PySCF overlap invariant "
f"(residual {worst:.2e}); please report this case."
)
self.last_gauge_residual = worst
self.last_dropped = 0
def strip(label):
return label.split("(")[0]
levels = {"left": [], "mo": [], "right": []}
self.fock = {}
self.last_symbreak = {}
# --onsite solves the crystal first: the fragment columns are the
# sublattice blocks of its Fock operator and need it in hand
order = (("mo", "left", "right") if self.onsite
else ("left", "right", "mo"))
for column in order:
if self.onsite and column != "mo":
energies, coefficients, dropped = self._block_bands(column, S)
self.last_dropped += dropped
else:
energies, coefficients = self._bands_at(column, kpoint)
if self.symmetrize:
energies, coefficients, deviation = self._symmetrized_bands(
column, energies, coefficients, S, representation)
self.last_symbreak[column] = deviation
if column == "mo":
# F(k) = S C E C+ S, exact for the eigenvectors of F with metric S
self.fock["mo"] = (
S @ coefficients @ np.diag(energies) @ coefficients.conj().T @ S
)
counts: dict[str, int] = {}
for energy, group in self._adaptive_groups(
energies, coefficients, S, representation, irreps):
if (column != "mo"
and self._ghost_fraction(column, group, S)
> GHOST_FRACTION_THRESHOLD):
continue
for irrep_label, space in self._irrep_split(
group, S, representation, irreps, labels
):
name = strip(irrep_label)
detail = ""
ghost_fraction = 0.0
if column == "mo":
# "GM4- #2" = second GM4- multiplet from the bottom;
# "GM4-(2)" read like a degeneracy count
counts[name] = counts.get(name, 0) + 1
label = f"{name} #{counts[name]}"
else:
spec = self._dominant_spec(space, S, column)
label = f"{spec.element} {spec.shell} {name}"
if not self.no_ghost:
ghost_fraction = self._ghost_fraction(column, space, S)
if ghost_fraction >= 0.005:
detail = ("weight on the removed sublattice's "
"(ghost) basis: "
f"{100 * ghost_fraction:.0f}% "
"(variational tail toward the point "
"charges; not drawn in the sketch)")
level = DiagramLevel(
level_id=f"{column}{len(levels[column])}",
column=column,
energy=float(energy),
degeneracy=space.shape[1],
irrep=name,
label=label,
vectors=space,
detail=detail,
)
level.ghost_fraction = ghost_fraction
# shared references, for the population-scaled sketches
level.overlap = S
level.sqrt_overlap = S_half
levels[column].append(level)
if column != "mo":
# 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]}"
# fragment projection (level.composition): crystal levels in the
# Loewdin-orthogonalized fragment-level basis. Internal only -- it
# positions the connector lines and feeds the alignment anchors
# (absolute_composition keeps the unnormalized weights) but is NOT
# displayed as percentages: the fragment eigenstates mix AO shells
# among themselves (the "Sc 3d R5+" fragment level carries 4d AO
# character), so its weights disagree with the AO populations below
# and showing both confused more than it explained.
assign_fragment_compositions(levels, S)
for crystal in levels["mo"]:
# the ONE displayed composition (panel bars, hover tooltip and
# the terminal all quote this list): per-(element, shell) AO
# populations of the PySCF eigenvector in the selected projection
# (--projection lowdin/mulliken; Loewdin sums to exactly 100%
# with non-negative entries) -- the same partial-charge measure
# as `crystod --dos --pyscf`. Symmetry does the orbital
# selection: a state of irrep G only picks up the G-adapted
# combination of each shell (forbidden shells project to ~0), so
# every entry is labeled with the crystal irrep.
shares = self._population_shares(
crystal.vectors, crystal.degeneracy, S, S_half)
crystal.display_composition = [
(f"{name} {crystal.irrep}", value) for value, name in shares
]
populations = ", ".join(f"{name} {100 * value:.1f}%"
for value, name in shares)
crystal.detail = f"{self.projection_label}: {populations}"
# the fragment/onsite columns carry the same displayed composition,
# but list ONLY their own sublattice's shells -- a fragment level is
# a sublattice state, and per-shell entries of the other element
# read as contamination ("O states inside a Ti-only level"). The
# weight that does sit beyond the own shells (the counterpoise
# ghost tail in the fragment-SCF mode, and in every mode the
# Loewdin attribution of the overlap density) is aggregated into
# one closing note instead.
for column in ("left", "right"):
own = self.side_specs[column]
for level in levels[column]:
shares = self._population_shares(
level.vectors, level.degeneracy, S, S_half, specs=own)
level.display_composition = [
(f"{name} {level.irrep}", value) for value, name in shares
]
row = self.projection_label + ": " + ", ".join(
f"{name} {100 * value:.1f}%" for value, name in shares)
remainder = 1.0 - sum(value for value, _ in shares)
if remainder >= 0.005:
if self.onsite:
row += (f" (+{100 * remainder:.1f}% "
f"{self.projection_label}-attributed to the "
"other sublattice's basis: overlap density; "
"the state has no coefficients there)")
else:
row += (f" (+{100 * remainder:.1f}% on the removed "
"sublattice's basis: ghost tail + "
f"{self.projection_label} attribution)")
level.detail = (row if not level.detail
else f"{row}\n{level.detail}")
# left-right coupling <phi_left| F(k) |phi_right> of the crystal Fock
# operator: same irrep and a large matrix element compared with the
# level separation is what actually splits bonding from antibonding
self.last_coupling = []
fock = self.fock["mo"]
for left in levels["left"]:
for right in levels["right"]:
if left.irrep != right.irrep:
continue
# the fragment states are NOT mutually orthogonal, so the bare
# matrix element h = <phi_L|F|phi_R> carries an
# overlap-times-mean-energy part s*(e_L+e_R)/2 that mimics a
# coupling even between non-interacting states (diffuse F 3s/3d
# fragment virtuals at +50 eV showed |h| of 5-15 eV against
# semicore levels this way). The reported strength is the
# first-order Loewdin-orthogonalized coupling
# H~ = h - s (e_L + e_R)/2
# with e_X = <phi_X|F|phi_X> the crystal-Fock expectations --
# invariant under the G=0 reference (F -> F + V S shifts h by
# V s and both e_X by V).
overlap_block = left.vectors.conj().T @ S @ right.vectors
raw_block = left.vectors.conj().T @ fock @ right.vectors
mean_energy = 0.5 * (
float(np.trace(left.vectors.conj().T @ fock
@ left.vectors).real) / left.degeneracy
+ float(np.trace(right.vectors.conj().T @ fock
@ right.vectors).real) / right.degeneracy)
block = raw_block - mean_energy * overlap_block
strength = float(np.sqrt(np.sum(np.abs(block) ** 2) / left.degeneracy))
raw_strength = float(np.sqrt(
np.sum(np.abs(raw_block) ** 2) / left.degeneracy))
overlap_norm = float(np.sqrt(
np.sum(np.abs(overlap_block) ** 2) / left.degeneracy))
gap = abs(left.energy - right.energy)
# two-level mixing: tan(2 theta) = 2|H| / dE, so the minority
# weight of each mixed orbital is sin^2(theta)
angle = 0.5 * np.arctan2(2.0 * strength, gap)
self.last_coupling.append(
(left, right, strength, gap, float(np.sin(angle) ** 2),
overlap_norm, raw_strength)
)
self.last_coupling.sort(key=lambda item: -item[4])
# COOP bonding character of every crystal level (see
# assign_bond_characters); the occupations must be filled first
self._fill(levels["mo"], self.electrons)
for column in ("left", "right"):
self._fill(levels[column], self.side_electrons[column])
spec_ranges = {(spec.element, spec.shell): self.spec_indices[id(spec)]
for spec in self.specs}
assign_bond_characters(
levels, S, self.embed_rows["left"], self.embed_rows["right"],
spec_ranges, sqrt_overlap=S_half, hamiltonian=fock,
)
return levels, labels
# --------------------------------------------------------------------- report
def describe_chk(path: str) -> None:
"""Print the calculation conditions stored in a --chk checkpoint.
The file is a compressed npz (binary for size and load speed); this is
the human-readable window into it: the defining parameters, what it
holds, and a ready-to-paste option string that reproduces them
(``crystod --chk-info FILE``).
"""
import json
import os
if not os.path.exists(path):
raise SystemExit(f"ERROR: {path} does not exist.")
try:
data = np.load(path, allow_pickle=False)
params = json.loads(str(data["params"]))
except Exception as error:
raise SystemExit(
f"ERROR: {path} is not a CrystOD --chk checkpoint ({error}).")
stored = ([str(c) for c in data["energy_columns"]]
if "energy_columns" in data.files
else ["mo", "left", "right"])
energies = np.asarray(data["energies"], dtype=float)
lattice = np.asarray(data["lattice"], dtype=float)
lengths = np.linalg.norm(lattice, axis=1)
n_kpts, n_ao = np.asarray(data["dm_mo"]).shape[:2]
smeared = sorted(str(s) for s in data["smeared"])
kind = ("full three-SCF run"
if all(f"dm_{c}" in data.files for c in ("mo", "left", "right"))
else "crystal density only (--onsite run)")
oxidation = " ".join(f"{el}={q:+g}"
for el, q in sorted(params["oxidation"].items()))
size = os.path.getsize(path) / 1e6
print(f" * {path} -- CrystOD SCF checkpoint "
f"(WAVECAR-style restart, {size:.1f} MB) *")
print(f" structure : {params['left']} + {params['right']}, "
f"{len(params['symbols'])} atoms/cell, "
f"a = {lengths[0]:.4f} / {lengths[1]:.4f} / {lengths[2]:.4f} A")
print(f" method : {params['xc'].upper()} / {params['basis']} / "
f"{params['pseudo']}, ke_cutoff {params['ke_cutoff']:g} Ha, "
f"k-mesh {'x'.join(map(str, params['kmesh']))}"
+ (f", max_l {params['max_l']}"
if params.get("max_l") is not None else "")
+ (f", smearing sigma {params['sigma']:g} eV"
if params.get("sigma") else ""))
print(f" electrons : {params['electrons']:g} per cell "
f"(oxidation {oxidation})")
print(f" fragments : left {params['left']} | right {params['right']}"
+ (", own-sublattice basis (--no-ghost)"
if params.get("no_ghost") else ", counterpoise ghosts"))
print(f" contents : {kind}; densities on {n_kpts} k points x "
f"{n_ao} AOs")
print(" energies : " + " | ".join(
f"{column} {energies[index]:.8f} Ha"
for index, column in enumerate(stored)))
if smeared:
print(f" smearing : {', '.join(smeared)} carried Fermi smearing "
"when written")
reuse = (f"--co-left {params['left']} --co-right {params['right']} "
f"--xc {params['xc']} --basis {params['basis']} "
f"--pseudo {params['pseudo']} "
f"--kmesh {' '.join(map(str, params['kmesh']))} "
f"--ke-cutoff {params['ke_cutoff']:g}"
+ (f" --max-l {params['max_l']}"
if params.get("max_l") is not None else "")
+ (" --no-ghost" if params.get("no_ghost") else "")
+ (f" --sigma {params['sigma']:g}"
if params.get("sigma") else "")
+ f" --oxidation {oxidation}"
# a crystal-only checkpoint can only feed --onsite runs
+ ("" if kind.startswith("full") else " --onsite"))
print(f" reuse with: {reuse} --chk {path}")
def report_and_write(cell, *, left, right, symprec, electrons, kpoint_filter,
output_path, structure_label, oxidation=None, basis=None,
pseudo=None, xc="pbe", kmesh=None, ke_cutoff=200.0,
sigma=0.0, degeneracy_tol=None, align=True, no_ghost=False,
symmetrize=True, max_l=None, projection="lowdin",
chk=None, no_chk=False, onsite=False,
conventional=False, verbose=0):
"""Terminal report + HTML for the PySCF crystal-orbital diagram."""
diagram = PySCFCrystalOrbitalDiagram(
cell, left, right, symprec=symprec, electrons=electrons,
oxidation=oxidation, basis=basis or "gth-dzvp-molopt-sr",
pseudo=pseudo or "gth-pbe", xc=xc, kmesh=kmesh, ke_cutoff=ke_cutoff, sigma=sigma,
degeneracy_tol=degeneracy_tol, no_ghost=no_ghost, symmetrize=symmetrize,
max_l=max_l, projection=projection, chk=chk, onsite=onsite,
conventional=conventional, verbose=verbose,
)
if chk is None and not no_chk:
# --pyscf runs are expensive, so the three converged densities are
# cached by default under the compound's name; a run whose options
# do not match the stored ones simply recomputes and overwrites
# (an explicit --chk file stays strict -- see _chk_reject)
diagram.chk_path = f"CHK_{diagram.compound_formula()}.chk"
diagram.chk_auto = True
dataset = diagram.builder.spglib_dataset
print("\n * Space group *")
print(f" {dataset['international']} ({dataset['number']})\n")
# ---- stage 1: the valence orbitals the pseudopotential leaves ----------
print(" * Valence basis (pseudopotential valence shells + polarization) *")
for column in ("left", "right"):
for element in dict.fromkeys(
spec.element for spec in diagram.side_specs[column]
):
shells = [spec.shell for spec in diagram.side_specs[column]
if spec.element == element]
sites = len({site for spec in diagram.side_specs[column]
if spec.element == element for site in spec.sites})
print(f" {element:<3} {' '.join(shells):<32} x{sites} site(s)")
print(f" basis {diagram.basis_name} / pseudo {diagram.pseudo_name} / "
f"functional {diagram.xc.upper()}")
print(f" AO space shared by all three calculations: {diagram.n_ao} orbitals\n")
# ---- the fragments ----------------------------------------------------
if diagram.onsite:
print(" * Fragment columns (--onsite): sublattice blocks of the "
"crystal Fock *")
for column in ("left", "right"):
print(f" {column:<5} {diagram.formula[column]:<8} "
f"{diagram.side_electrons[column]} electrons (formal "
"count), levels = per-shell on-site multiplets "
"<phi|F(k)|phi> of its own shells")
print(f" crystal {diagram.formula['left'] + diagram.formula['right']:<7} "
f"{diagram.crystal_electrons} electrons "
f"= {diagram.side_electrons['left']} + "
f"{diagram.side_electrons['right']}")
print(" ONE Hamiltonian for every column (no fragment SCF, no point\n"
" charges, no reference alignment): a crystal level's rise or\n"
" drop against its parents is purely the left-right orbital\n"
" interaction\n")
else:
neutral = not any(diagram.oxidation.values())
print(" * Fragments ("
+ ("neutral sublattices" if neutral else "formal-charge ions")
+ " + ghost basis of the removed sublattice) *")
other = {"left": "right", "right": "left"}
for column in ("left", "right"):
if neutral:
felt = " (all oxidation states 0: no point-charge lattice)"
else:
felt = ", in the " + " + ".join(
f"{element}^{diagram.oxidation[element]:+g}"
for element in dict.fromkeys(
diagram.symbols[i]
for i in diagram.side_atoms[other[column]]
)
) + " point-charge lattice"
print(f" {column:<5} {diagram.formula[column]:<8} charge "
f"{diagram.side_charge[column]:+g}, "
f"{diagram.side_electrons[column]} electrons{felt}")
print(f" crystal {diagram.formula['left'] + diagram.formula['right']:<7} "
f"{diagram.crystal_electrons} electrons "
f"= {diagram.side_electrons['left']} + "
f"{diagram.side_electrons['right']}")
print(" every cell is neutral (no monopole divergence), but each calculation\n"
" still pins its own G=0 average potential to zero, so the raw columns\n"
" are offset by one constant each -- removed below by deep-level alignment")
if diagram.onsite:
pass
elif diagram.no_ghost:
print(" fragment basis: OWN sublattice only (--no-ghost) -- the removed\n"
" sublattice's functions are excluded from the variational space,\n"
" so no fragment wave function can sit on the removed atoms\n")
else:
print(" fragment basis: counterpoise (ghost functions of the removed\n"
" sublattice kept); sketches draw own-sublattice components only,\n"
" the ghost weight of each level is reported numerically\n")
from pyscf import lib as _pyscf_lib
print(f" * Self-consistent calculations ({'x'.join(map(str, diagram.kmesh))} "
f"k-mesh, ke_cutoff {diagram.ke_cutoff} Hartree, "
f"{_pyscf_lib.num_threads()} OpenMP threads) *")
diagram.run()
kpoints = diagram.special_kpoints()
if kpoint_filter is not None:
available = [name for name, _ in kpoints]
kpoints = [(name, k) for name, k in kpoints if name == kpoint_filter]
if not kpoints:
raise SystemExit(
f"ERROR: k point '{kpoint_filter}' is not a special point of this "
f"space group (available: {', '.join(available)})."
)
diagram.prepare_bands([kpoint for _, kpoint in kpoints])
# ---- solve every k point first; printing follows the alignment ---------
records = []
for name, kpoint in kpoints:
levels, _ = diagram.solve_at(kpoint)
irreps, mapping, labels, representation = diagram.little_group_data(kpoint)
records.append({
"name": name,
"kpoint": kpoint,
"levels": levels,
"coupling": list(diagram.last_coupling),
"residual": diagram.last_gauge_residual,
"symbreak": dict(getattr(diagram, "last_symbreak", {})),
"content": diagram.site_symmetry_irreps(
kpoint, representation, irreps, labels),
})
# ---- deep-level alignment ----------------------------------------------
if diagram.onsite:
align = False
print("\n * Single-Hamiltonian mode (--onsite): every column is "
"measured on the\n crystal Fock's own scale -- no alignment "
"step needed *")
if diagram.onsite:
pass
elif align:
shifts, anchors = diagram.align_fragment_columns(records)
print("\n * Deep-level alignment (pre-bonding reference) *")
if shifts is None:
print(" no chemically inert fragment level found; columns left "
"on their raw G=0 references")
else:
reference = shifts["reference"]
anchor = anchors[reference]
print(f" anchor : {anchor['label']} of the {reference} fragment "
f"({anchor['fragment_energy']:.2f} eV, counterpart purity "
f"{100 * anchor['purity']:.0f}%, "
f"k-spread {anchor['spread']:.3f} eV over {anchor['n_k']} k)")
print(f" shifts : left {shifts['left']:+.3f} eV | "
f"crystal {shifts['mo']:+.3f} eV | "
f"right {shifts['right']:+.3f} eV")
other = {"left": "right", "right": "left"}[reference]
if anchors.get(other):
info = anchors[other]
print(f" ({other} anchored via {info['label']}, purity "
f"{100 * info['purity']:.0f}%, k-spread {info['spread']:.3f} eV"
+ ("; WARNING: purity below "
f"{100 * ALIGNMENT_PURITY:.0f}%, the anchor mixes and the "
"column offset inherits that error" if info["fallback"] else "")
+ ")")
print(" each calculation pins its own average (G=0) potential to zero;"
" the columns\n are shifted so this inert deep level keeps its"
" pre-bonding energy -- XPS-style")
else:
print("\n * Deep-level alignment disabled (--no-align): each column keeps "
"its own G=0 reference *")
# ---- stage 0: isolated formal-charge ions (outermost columns) ----------
if not diagram.onsite:
diagram.atomic_ion_levels()
ao_anchors = diagram.attach_atomic_columns(records)
print("\n * Isolated formal-charge ions (outermost columns) *")
for element, ion in diagram.atomic_ions.items():
shells = ", ".join(f"{name} {energy:+.2f}"
for name, _, energy in ion["shells"])
print(f" {element}^{ion['charge']:+d} ({ion['method']}): "
f"{shells} eV (vacuum)")
anchor_name, shift = ao_anchors[element]
if anchor_name == "none":
print(f" -> WARNING: no fragment level is dominated "
f"by a {element} shell; the ion column keeps the "
"raw molecular vacuum reference (no anchor)")
else:
print(f" -> column shifted {shift:+.2f} eV: the "
"deepest shell with a fragment counterpart "
f"({element} {anchor_name}) is anchored to its "
"sublattice band center (molecular vacuum and "
"periodic G=0 references share no common zero)")
caveats = ion.get("caveats", {})
virtuals = [n for n, kind in caveats.items()
if kind == "virtual"]
continuum = [n for n, kind in caveats.items()
if kind == "continuum"]
if virtuals:
if diagram.richer_basis:
remedy = (f"{diagram.basis_name} has no diffuse "
"functions; a richer basis, e.g. --basis "
f"{diagram.richer_basis}, removes the "
"basis-side error")
elif "molopt" in diagram.basis_name.lower():
remedy = (f"{diagram.basis_name} has no diffuse "
"functions, and PySCF ships no richer GTH "
"set for these elements")
else:
remedy = ("limited by the basis's diffuse coverage "
"and the Kohn-Sham virtual error")
print(f" note: {', '.join(virtuals)}: finite-basis "
"virtuals, not physical Rydberg levels "
f"({remedy})")
if continuum:
print(f" note: {', '.join(continuum)}: discretized "
"continuum -- one more electron on the free anion "
"is unbound; in the crystal the Madelung potential "
"provides the binding")
# ---- per-k-point report -------------------------------------------------
for record in records:
name, kpoint, levels = record["name"], record["kpoint"], record["levels"]
print(f"\n * k point {name} {_format_kpoint(kpoint)} *")
print(f" (AO representation verified against the PySCF overlap: "
f"max |D+SD - S| = {record['residual']:.1e})")
if record.get("symbreak"):
parts = ", ".join(f"{col} {1000 * dev:.0f} meV"
for col, dev in record["symbreak"].items())
print(" raw-SCF point-group breaking, removed by Fock "
f"group-averaging: {parts}")
for col_name, column in (("crystal", "mo"),
(diagram.formula["left"], "left"),
(diagram.formula["right"], "right")):
for lv in levels.get(column, []):
if getattr(lv, "partial", False):
print(f" WARNING: aufbau leaves the {col_name} "
f"column's {lv.irrep} level at {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 "
"of the model")
# ---- stage 2: site-symmetry induced irreps -----------------------
print(" site-symmetry induced representations:")
for (element, shell), parts in record["content"].items():
if parts:
print(f" {element} {shell:<4} = {' + '.join(parts)}")
for column in ("left", "right"):
parts = ", ".join(
f"{lv.label} ({lv.energy:.2f}"
+ (f", ghost {100 * getattr(lv, 'ghost_fraction', 0.0):.0f}%"
if getattr(lv, "ghost_fraction", 0.0) >= 0.05 else "")
+ ")"
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 " "
print(f" {lv.label:<10} {lv.energy:9.2f} eV x{lv.degeneracy}"
f" {occupancy:<4} {composition}")
if record["coupling"]:
# dE and the mixing fraction on the ALIGNED scale -- the raw
# separations carried the per-calculation reference offsets
rescored = []
for (lv_left, lv_right, strength, _gap, _mix,
overlap_norm, raw_strength) in record["coupling"]:
gap = abs(lv_left.energy - lv_right.energy)
mixing = float(np.sin(0.5 * np.arctan2(2.0 * strength, gap)) ** 2)
rescored.append((lv_left, lv_right, strength, gap, mixing,
overlap_norm, raw_strength))
rescored.sort(key=lambda item: -item[4])
record["coupling_aligned"] = rescored
print(" same-irrep couplings, Loewdin-corrected "
"|H~| = |<L|F|R> - S (e_L+e_R)/2| "
"(mix = sin^2 of the two-level mixing angle):")
for (lv_left, lv_right, strength, gap, mixing,
overlap_norm, _raw) in rescored[:8]:
print(f" {lv_left.label:<16} x {lv_right.label:<16} "
f"|H~| = {strength:6.2f} eV |S| = {overlap_norm:5.3f}"
f" dE = {gap:7.2f} eV mix = {100 * mixing:4.1f}%")
# the terminal shows only the top-8 couplings per k point; the full list
# is a result worth keeping, so it is written next to the HTML
coupling_path = output_path
if coupling_path.endswith(".html"):
coupling_path = coupling_path[: -len(".html")]
coupling_path += "_coupling.txt"
with open(coupling_path, "w") as handle:
scale_note = ("energies on the crystal Fock's own scale (--onsite)"
if diagram.onsite else "energies on the aligned scale")
handle.write(
"# same-irrep couplings between the fragment levels, from the\n"
f"# converged crystal Fock operator F(k); {scale_note}\n"
"# |S| = |<phi_L| S |phi_R>|: the fragment states are NOT\n"
"# mutually orthogonal\n"
"# |H|raw = |<phi_L| F |phi_R>|: carries an unphysical\n"
"# overlap-times-mean-energy part |S|*(e_L+e_R)/2\n"
"# |H~| = |<phi_L|F|phi_R> - S (e_L+e_R)/2|, e_X = <phi_X|F|phi_X>:\n"
"# first-order Loewdin-orthogonalized coupling, the\n"
"# resonance integral to reason with\n"
"# mix = sin^2(theta) with tan(2 theta) = 2|H~| / dE "
"(two-level mixing fraction)\n"
f"# columns: irrep left_level E_left(eV) right_level "
"E_right(eV) |S| |H|raw(eV) |H~|(eV) dE(eV) mix(%)\n")
for record in records:
handle.write(f"\n# k point {record['name']} "
f"{_format_kpoint(record['kpoint'])}\n")
for (lv_left, lv_right, strength, gap, mixing,
overlap_norm, raw_strength) in record.get(
"coupling_aligned", []):
handle.write(
f"{lv_left.irrep:<6} {lv_left.label:<18} "
f"{lv_left.energy:9.3f} {lv_right.label:<18} "
f"{lv_right.energy:9.3f} {overlap_norm:6.3f} "
f"{raw_strength:8.3f} {strength:8.3f} {gap:8.3f} "
f"{100 * mixing:7.2f}\n")
entries = [(r["name"], r["kpoint"], r["levels"]) for r in records]
write_crystal_diagram_html(diagram, entries, output_path, structure_label)
print(f"\nCrystal-orbital diagram written to {output_path}")
print(f"Same-irrep couplings written to {coupling_path}")
def main(argv: list[str] | None = None) -> None:
import argparse
from pathlib import Path
from .crystal_orbital_diagram import parse_oxidation_tokens
parser = argparse.ArgumentParser(
description="Crystal-orbital diagram from three periodic PySCF calculations."
)
parser.add_argument("--poscar", default="POSCAR")
parser.add_argument("--co-left", nargs="+", required=True, metavar="FORMULA")
parser.add_argument("--co-right", nargs="+", required=True, metavar="FORMULA")
parser.add_argument("--oxidation", nargs="+", default=None, metavar="EL=Q")
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)
parser.add_argument("--basis", default="gth-dzvp-molopt-sr")
parser.add_argument("--pseudo", default="gth-pbe")
parser.add_argument("--xc", default="pbe")
parser.add_argument("--kmesh", type=int, nargs=3, default=None,
metavar=("N1", "N2", "N3"))
parser.add_argument("--ke-cutoff", type=float, default=200.0)
parser.add_argument("--sigma", type=float, default=0.0,
help="Fermi smearing width in eV (0 = integer "
"occupations; a cell with an odd electron "
"count always smears, 0.2 eV unless --sigma "
"is positive)")
parser.add_argument("--no-symmetrize", action="store_true",
help="do not re-diagonalize the group-averaged Fock (debug: shows "
"the raw grid-broken degeneracies)")
parser.add_argument("--max-l", type=int, default=None,
help="drop basis shells with l above this from every element "
"(e.g. 2 removes the f polarization functions)")
parser.add_argument("--no-ghost", action="store_true",
help="exclude the removed sublattice's basis functions from the "
"fragment calculations entirely (hard constraint: no fragment "
"wave function on the removed atoms; loses counterpoise "
"consistency)")
parser.add_argument("--no-align", action="store_true",
help="keep each calculation's own G=0 reference instead of "
"the deep-level (XPS-style) column alignment")
parser.add_argument("--degeneracy-tol", type=float, default=None,
help="seed window in eV for clustering degenerate levels "
f"(default {DEGENERACY_SEED_EV}; groups are then merged "
"until the irrep multiplicities are integral)")
parser.add_argument("--projection", choices=("lowdin", "mulliken"),
default="lowdin",
help="population measure for the sketch lobe sizes and "
"the per-(element, shell) rows: Loewdin |S^(1/2)c|^2 "
"(default; non-negative, sums to 100%%) or Mulliken "
"gross populations Re[c*(Sc)]")
parser.add_argument("--no-chk", action="store_true",
help="do not cache the SCFs in CHK_{formula}.chk")
parser.add_argument("--chk", default=None, metavar="FILE",
help="WAVECAR-style restart file: written after the "
"SCFs if missing, read (skipping all three SCFs) if "
"present; parameters are verified before reuse")
parser.add_argument("--onsite", action="store_true",
help="single-Hamiltonian mode: only the crystal SCF "
"runs, and the fragment columns are the per-shell "
"on-site multiplets <phi|F|phi> of the crystal Fock "
"(tight-binding on-site energies, one level per "
"induced irrep) -- no point charges, no reference "
"alignment; a full-run --chk is reused, only the "
"crystal density is read")
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 (display-only; compatible "
"with any --chk)")
parser.add_argument("--output", default=None)
parser.add_argument("--tolerance", type=float, default=1e-5)
parser.add_argument("--verbose", type=int, default=0)
args = parser.parse_args(argv)
from .star_of_k import read_poscar_or_exit
cell = read_poscar_or_exit(args.poscar)
stem = Path(args.poscar).name
for extension in (".vasp", ".poscar"):
if stem.lower().endswith(extension):
stem = stem[: -len(extension)]
report_and_write(
cell,
left=args.co_left,
right=args.co_right,
symprec=args.tolerance,
electrons=args.electrons,
kpoint_filter=args.kpoint,
output_path=args.output or f"CrystOD_{stem}_pyscf.html",
structure_label=stem,
oxidation=(parse_oxidation_tokens(args.oxidation) if args.oxidation else None),
basis=args.basis,
pseudo=args.pseudo,
xc=args.xc,
kmesh=args.kmesh,
ke_cutoff=args.ke_cutoff,
sigma=args.sigma,
degeneracy_tol=args.degeneracy_tol,
align=not args.no_align,
no_ghost=args.no_ghost,
onsite=args.onsite,
conventional=args.conventional,
symmetrize=not args.no_symmetrize,
max_l=args.max_l,
projection=args.projection,
chk=args.chk, no_chk=args.no_chk,
verbose=args.verbose,
)
if __name__ == "__main__":
main()
