Parent lattices and sites

SymmetryOp{D}

A space-group operation in lattice coordinates: a rotation R (D×D integer matrix) plus a fractional translation t (length-D float vector). For a single Bravais lattice, all t == 0. For a multilattice (more than one position in the dset), the dset induces fractional translations (screws / glides) that must be tracked alongside the rotation.

SymmetryOp{3} is the type of every element of ParentLattice{3}.space_group. The parametric D keeps the type ready for 2D or 4D extensions.

Thin wrapper around Spacey's SpacegroupOp. The fractional translation t is canonicalized to [0,1)^D at construction (delegated to Spacey._canonicalize_τ); we don't need to wrap again.

Examples

The identity element of any 3D point group has R = I and t = 0:

julia> id = SymmetryOp{3}([1 0 0; 0 1 0; 0 0 1], [0.0, 0.0, 0.0])
SymmetryOp{3}([1 0 0; 0 1 0; 0 0 1], [0.0, 0.0, 0.0])

julia> id.R
3×3 Matrix{Int64}:
 1  0  0
 0  1  0
 0  0  1

ParentLattice{D}

The geometric description of the parent multilattice for an enumeration: basis vectors A, dset dset (basis sites in fractional coordinates), and the cached space_group of the multilattice (rotation + fractional-translation pairs).

The dset captures the multilattice basis — for a Bravais lattice, length==1; for HCP, length==2; for perovskite ABO₃, length==5. The dset does not need to contain the origin — placing the origin where it makes physical sense (e.g., at the inversion center for diamond) is a user choice and the enumeration math doesn't require the dset to include the origin.

What the constructor canonicalizes silently

1. Periodic-coordinate wrap. Each dset position is folded into [0,1)^D via mod(., 1). So [1.5, -0.5, 0.5] becomes [0.5, 0.5, 0.5] — mathematically equivalent under lattice translation. (Matches the convention in ASE / pymatgen.)
2. Bravais origin shift (only when length(dset) == 1). A single-site dset has a degenerate choice of origin — there's no geometric structure picking one position over another. We shift the lone dset entry to the origin, so the resulting space_group doesn't carry artifact t-translations introduced by the user's choice of origin. For multilattices (length(dset) ≥ 2) we never shift — the relative positions encode physically meaningful structure (placing the origin at diamond's inversion center is the canonical example).

Numerical scale check

The basis is rejected as singular if |det(A)| / prod(‖aⱼ‖) ≤ 1e-12. This is the Hadamard ratio — 1 for orthogonal columns, 0 for linearly dependent columns, dimensionless and therefore unit-independent. So the check works whether the user works in Ångströms, nm, or meters; what it catches is geometric near-singularity, not absolute determinant magnitude.

Space group

Computed once at construction by calling Spacey.spacegroup(c::Crystal) with a uniform-species Crystal built from (A, dset). Cached for ergonomics — call sites read parent.space_group without re-invoking Spacey.

dset-permutation precompute (multilattice support)

For each space-group operation (N, t), the constructor precomputes the dset permutation π(i) and the integer lattice shifts v_i satisfying N·d_i + t = d_{π(i)} + v_i. Stored as dset_perms::Vector{Vector{Int}} and dset_shifts::Vector{Vector{Vector{Int}}}, both indexed by op. The multilattice supercell-permutation-group construction consumes these to build the correct n_D · n-site permutation group; for single-lattice parents (n_D = 1) π is trivially [1] and v is [zeros(Int, D)], and the multilattice path degenerates to the single-lattice path. Tolerance eps_dset (default 1e-6) is configurable via the constructor.

The parametric D is the spatial dimension. Almost all uses are D=3; D=2 is reserved for the future surface/2D extension. The constructor infers D from size(A,1).

basis(p::ParentLattice{D}) -> Matrix{Float64}

Return the basis matrix of the parent lattice. Columns are basis vectors in Cartesian coordinates; the matrix is D×D.

Examples

julia> p = ParentLattice([0.0 0.5 0.5; 0.5 0.0 0.5; 0.5 0.5 0.0]);

julia> basis(p)
3×3 Matrix{Float64}:
 0.0  0.5  0.5
 0.5  0.0  0.5
 0.5  0.5  0.0

dset(p::ParentLattice{D}) -> Vector{Vector{Float64}}

Return the dset (basis sites in fractional coordinates, canonicalized to [0,1)^D). Length 1 for a Bravais lattice (one site, shifted to the origin), length ≥ 2 for a multilattice. See ParentLattice for canonicalization details.

Examples

julia> p = ParentLattice([1.0 0 0; 0 1 0; 0 0 1], [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]]);

julia> dset(p)
2-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0]
 [0.5, 0.5, 0.5]

space_group(p::ParentLattice{D}) -> Vector{SymmetryOp{D}}

Return the cached multilattice space group (rotation + fractional-translation pairs) of the parent. Computed once at construction by Spacey.

Examples

julia> p = ParentLattice([0.0 0.5 0.5; 0.5 0.0 0.5; 0.5 0.5 0.0]);

julia> length(space_group(p))
48

ndset(p::ParentLattice) -> Int

Return the number of sites in the dset — i.e., length(dset(p)).

Examples

julia> p = ParentLattice([1.0 0 0; 0 1 0; 0 0 1], [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]]);

julia> ndset(p)
2

n_nonzero_translations(p::ParentLattice; tol::Real = 1e-9) -> Int

Count symmetry operations whose fractional translation t has any component with |t_i| > tol. Distinguishes symmorphic space groups (all t == 0, so this returns 0) from non-symmorphic ones (screw axes / glide planes, where some operations carry an intrinsic translation).

Examples

julia> p_sc = ParentLattice([1.0 0 0; 0 1 0; 0 0 1]);  # simple cubic — Pm-3m, symmorphic

julia> n_nonzero_translations(p_sc)
0

julia> A_hcp = [1.0 -0.5 0.0; 0.0 sqrt(3)/2 0.0; 0.0 0.0 sqrt(8/3)];

julia> p_hcp = ParentLattice(A_hcp, [[0.0, 0.0, 0.0], [1/3, 2/3, 1/2]]);  # HCP — P6_3/mmc

julia> n_nonzero_translations(p_hcp)
12

lattice_rotations(p::ParentLattice{D}) -> Vector{Matrix{Int}}

Project out just the rotation parts of p.space_group, dropping the fractional translations. This is what callers like getSymInequivHNFs and Supercell's constructor want — the HNF symmetry equivalence and the supercell stabilizer detection are pure-rotation tests; fractional translations don't enter.

Returns a fresh Vector{Matrix{Int}} (one allocation per call). Cheap enough for typical use; if profiling motivates it, we can later cache this on ParentLattice itself.

Examples

julia> p = ParentLattice([1.0 0 0; 0 1 0; 0 0 1]);  # simple cubic, point group order 48

julia> rots = lattice_rotations(p);

julia> length(rots)
48

julia> rots[1]  # the identity is first
3×3 Matrix{Int64}:
 1  0  0
 0  1  0
 0  0  1

Site{D}

A single position in the parent cell where atomic substitution can happen. Carries the position itself plus the set of species labels allowed at this position.

The position is in fractional coordinates of the parent lattice (ParentLattice{D}.A), matching the convention used by ParentLattice{D}.dset. Allowed labels are stored as a BitSet of integers from 0:k-1, where k is the number of species in the problem. (D is typically 3, for 3D crystals.)

A site is inactive if length(allowed_labels) == 1 — only one species can occupy it, so it has no configurational freedom and gets stripped from the labeling space during enumeration. A site is active otherwise.

Site does not validate the position against any specific ParentLattice. Cross-validation happens at the enumerate(parent, sites) boundary — keeps Site parametric on D only.

Examples

A binary substitution site (active — two allowed labels) and a fixed-species site (inactive — one allowed label):

julia> Site([0.0, 0.0, 0.0], [0, 1])
Site{3}([0.0, 0.0, 0.0], species {0, 1})

julia> Site([0.5, 0.5, 0.5], [0])
Site{3}([0.5, 0.5, 0.5], species {0}  [inactive])

SymbolSite{D}

Intermediate type for constructing Sites with atomic-symbol labels. Holds a position and a Vector{Symbol} of allowed species (e.g. [:Al, :Ga]). The integer↔symbol mapping is established when the SymbolSite is folded into a Sites collection — the Sites constructor walks the list in order and assigns the first-seen symbol to label 0, the second-seen to label 1, etc.

Users typically don't construct a SymbolSite directly. It's the return type of Site(position, syms::AbstractVector{Symbol}):

Site([0.0, 0.0, 0.0], [:Al, :Ga])    # → SymbolSite{3}, not Site{3}

See Sites for how a Vector{SymbolSite{D}} becomes a fully- populated Sites with species_symbols.

Sites{D}

A collection of Site{D}s plus an equivalence relation declaring which sites must carry the same label across configurations.

Equivalencies are user-declared — they are NOT derived from the parent's space group. Use cases include slab geometries (mirror-image layers must share composition; the slab vacuum breaks the parent's 3D periodicity that would otherwise tie them) and any other physical constraint the user knows but the symmetry analysis can't see.

The equivalence relation is stored as a Union-Find (IntDisjointSets) over the site indices 1:length(list). This makes transitivity automatic by data structure: if you declare site i ↔ j and j ↔ k, the data structure correctly returns i, j, k as a single class without any further work on the user's part.

Two-variant constructor

• Incremental: build the Sites with no equivalencies first via Sites([Site(...), Site(...)]), then declare equivalencies one-by-one via equate!(sites, i, j).
• Upfront partition: Sites([Site(...), Site(...)], [[1,2], [3,4]]) validates the partition and builds the Union-Find in one shot. Suitable when the user already knows the equivalence classes from problem setup.

The two variants produce the same internal state and can be mixed (start with the upfront variant and call equate! later).

Examples

Upfront-partition: three sites with sites 1 and 2 tied (e.g., mirror-image slab layers); site 3 inactive.

julia> list = [Site([0.0, 0.0, 0.0], [0, 1]),
               Site([0.5, 0.5, 0.5], [0, 1]),
               Site([0.25, 0.25, 0.25], [0])];

julia> s = Sites(list, [[1, 2]]);

julia> n_active(s)
2

julia> n_canonical(s)
2

julia> canonical(s, 2)  # site 2 collapses to site 1's root
1

julia> n_effective(s)  # active AND canonical = 1
1

is_active(s::Site) -> Bool

A site is active iff it has more than one allowed label. Active sites are the ones the enumeration algorithm will assign labels to.

Examples

julia> is_active(Site([0.0, 0.0, 0.0], [0, 1]))
true

julia> is_active(Site([0.5, 0.5, 0.5], [0]))
false

is_inactive(s::Site) -> Bool

A site is inactive iff its allowed_labels has exactly one element — only one species can occupy it, so it contributes no configurational freedom to the enumeration.

Examples

julia> is_inactive(Site([0.0, 0.0, 0.0], [0, 1]))
false

julia> is_inactive(Site([0.5, 0.5, 0.5], [0]))
true

equate!(sites::Sites, i::Integer, j::Integer) -> Sites

Declare sites i and j equivalent in the user-supplied partition. Idempotent (equating already-equated sites is a no-op) and transitive (equating (i,j) then (j,k) puts all three in one class). Returns sites for chainability.

Examples

julia> s = Sites([Site([0.0, 0.0, 0.0], [0, 1]),
                  Site([0.5, 0.5, 0.5], [0, 1]),
                  Site([0.25, 0.25, 0.25], [0, 1])]);

julia> equate!(s, 1, 2);  # tie sites 1 and 2

julia> equate!(s, 2, 3);  # transitivity collapses all three into one class

julia> n_canonical(s)
1

canonical(sites::Sites, i::Integer) -> Int

Return the canonical (root) site index of the equivalence class containing site i. Sites in the same equivalence class share their root; the dispatcher uses one root per class as the labeling-space representative.

Examples

julia> s = Sites([Site([0.0, 0.0, 0.0], [0, 1]),
                  Site([0.5, 0.5, 0.5], [0, 1])], [[1, 2]]);

julia> canonical(s, 1)
1

julia> canonical(s, 2)  # same equivalence class → same root
1

active_canonical_sites(sites::Sites{D}) -> Vector{Tuple{Int, Site{D}}}

Return a vector of (index, Site{D}) pairs for sites that are both active (more than one allowed label) and canonical (the root of their equivalence class). This is the labeling space the enumeration algorithm sees: stripped of inactive sites and collapsed across equivalencies.

Examples

julia> list = [Site([0.0, 0.0, 0.0], [0, 1]),
               Site([0.5, 0.5, 0.5], [0, 1]),
               Site([0.25, 0.25, 0.25], [0])];

julia> s = Sites(list, [[1, 2]]);

julia> [i for (i, _) in active_canonical_sites(s)]  # site 1 is the only active canonical
1-element Vector{Int64}:
 1

n_active(sites::Sites)

Count of active sites (those with more than one allowed label).

n_canonical(sites::Sites)

Count of equivalence classes (number of distinct canonical roots).

n_effective(sites::Sites)

Count of active canonical sites — the actual dimension of the labeling space after stripping inactive sites and collapsing equivalencies. This is what enumerate(...) will use when it asks "how many free configurational variables does this problem have?"

species_symbols(sites::Sites) -> Union{Nothing, Vector{Symbol}}

Return the atomic-symbol mapping carried by sites, or nothing if the Sites was constructed with integer labels and no explicit mapping.

The returned vector is dense and 0-indexed by label: species_symbols(s)[i + 1] is the symbol for integer label i. Length equals k, the number of distinct labels across the sites.

julia> sites = Sites([
           Site([0.0, 0.0, 0.0], [:Al, :Ga]),
           Site([0.25, 0.25, 0.25], [:As]),
       ]);

julia> species_symbols(sites)
3-element Vector{Symbol}:
 :Al
 :Ga
 :As

julia> integer_sites = Sites([Site([0.0, 0.0, 0.0], [0, 1])]);

julia> species_symbols(integer_sites) === nothing
true

to_atom_labeling(structure, sites::Sites) -> Vector{Symbol}

Translate an integer labeling (to_labeling(structure) :: Vector{Int8}) to its atomic-symbol equivalent using sites.species_symbols. Throws ArgumentError if sites has no symbol mapping.

This is the natural read-back for users who constructed Sites with atomic symbols and want the post-enumeration labeling in those symbols rather than as raw integer labels.

julia> using Enumlib

julia> p = ParentLattice([0.0 0.5 0.5; 0.5 0.0 0.5; 0.5 0.5 0.0]);

julia> sites = Sites(p, [:Al, :Ga]);

julia> e = enumerate(p, sites; supercells = VolumeRange(1:1));

julia> to_atom_labeling(e[1], sites)
1-element Vector{Symbol}:
 :Ga