diff --git a/BandGraphs/Project.toml b/BandGraphs/Project.toml index 53676424..794f8fb9 100644 --- a/BandGraphs/Project.toml +++ b/BandGraphs/Project.toml @@ -10,7 +10,9 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Crystalline = "ae5e2be0-a263-11e9-351e-f94dad1eb351" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" GraphsMatching = "c3af3a8c-b79e-4b01-bf44-c718d7e0e0d6" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377" diff --git a/BandGraphs/build/connections_2d.jl b/BandGraphs/build/connections_2d.jl index 2406f878..a53d3810 100644 --- a/BandGraphs/build/connections_2d.jl +++ b/BandGraphs/build/connections_2d.jl @@ -114,6 +114,7 @@ CONNECTIONSD_2D[9] = C2[ C2(LK2(:Γ, "0, 0"), LK2(:Σ, "2α, 0")), # Γ = [0, 0] ↓ Σ = [2α, 0] C2(LK2(:Y, "1, 0"), LK2(:Σ, "2α, 0")), # Y = [1, 0] ↓ Σ = [2α, 0] C2(LK2(:S, "1, 2"), LK2(:Ω, "α, β")), # S = [1, 2] ↓ Ω = [α, β] (manual addition) + C2(LK2(:Γ, "0, 0"), LK2(:Ω, "α, β")), # Γ = [0, 0] ↓ Ω = [α, β] (manual addition) ] # Plane group 10 (p4) @@ -157,6 +158,7 @@ CONNECTIONSD_2D[14] = C2[ C2(LK2(:Γ, "0, 0"), LK2(:Σ, "α, 0")), # Γ = [0, 0] ↓ Σ = [α, 0] C2(LK2(:M, "1/2, 0"), LK2(:Σ, "α, 0")), # M = [1/2, 0] ↓ Σ = [α, 0] C2(LK2(:K, "1/3, 1/3"), LK2(:Ω, "α, β")), # K = [1/3, 1/3] ↓ Ω = [α, β] (manual addition) + C2(LK2(:Γ, "0, 0"), LK2(:Ω, "α, β")), # Γ = [0, 0] ↓ Ω = [α, β] (manual addition) ] # Plane group 15 (p31m) @@ -206,7 +208,7 @@ for sgnum in 1:MAX_SGNUM[2] lgirsd = lgirreps(sgnum, Val(2)) sg = spacegroup(sgnum, Val(2)) for timereversal in (false, true) - timereversal && (lgirsd = Dict(klab => realify(lgirs) for (klab, lgirs) in lgirsd)) + timereversal && (lgirsd = realify(lgirsd)) subts = SubductionTable{2}[] for c in cs push!(subts, SubductionTable(c, sg, lgirsd)) diff --git a/BandGraphs/data/connections/2d/connections-tr.jld2 b/BandGraphs/data/connections/2d/connections-tr.jld2 index 14e301ce..9129a0ef 100644 Binary files a/BandGraphs/data/connections/2d/connections-tr.jld2 and b/BandGraphs/data/connections/2d/connections-tr.jld2 differ diff --git a/BandGraphs/data/connections/2d/connections.jld2 b/BandGraphs/data/connections/2d/connections.jld2 index 14e301ce..9129a0ef 100644 Binary files a/BandGraphs/data/connections/2d/connections.jld2 and b/BandGraphs/data/connections/2d/connections.jld2 differ diff --git a/BandGraphs/data/connections/2d/subductions-tr.jld2 b/BandGraphs/data/connections/2d/subductions-tr.jld2 index c0b11a4d..20446fc2 100644 Binary files a/BandGraphs/data/connections/2d/subductions-tr.jld2 and b/BandGraphs/data/connections/2d/subductions-tr.jld2 differ diff --git a/BandGraphs/data/connections/2d/subductions.jld2 b/BandGraphs/data/connections/2d/subductions.jld2 index 5b71c715..6ca77cb2 100644 Binary files a/BandGraphs/data/connections/2d/subductions.jld2 and b/BandGraphs/data/connections/2d/subductions.jld2 differ diff --git a/BandGraphs/data/connections/3d/connections-tr.jld2 b/BandGraphs/data/connections/3d/connections-tr.jld2 index 9a050cb4..aff67059 100644 Binary files a/BandGraphs/data/connections/3d/connections-tr.jld2 and b/BandGraphs/data/connections/3d/connections-tr.jld2 differ diff --git a/BandGraphs/data/connections/3d/connections.jld2 b/BandGraphs/data/connections/3d/connections.jld2 index 765920c3..15db63af 100644 Binary files a/BandGraphs/data/connections/3d/connections.jld2 and b/BandGraphs/data/connections/3d/connections.jld2 differ diff --git a/BandGraphs/data/connections/3d/subductions-tr.jld2 b/BandGraphs/data/connections/3d/subductions-tr.jld2 index 68c52918..6ed659bf 100644 Binary files a/BandGraphs/data/connections/3d/subductions-tr.jld2 and b/BandGraphs/data/connections/3d/subductions-tr.jld2 differ diff --git a/BandGraphs/data/connections/3d/subductions.jld2 b/BandGraphs/data/connections/3d/subductions.jld2 index d2d6901d..d914f3e9 100644 Binary files a/BandGraphs/data/connections/3d/subductions.jld2 and b/BandGraphs/data/connections/3d/subductions.jld2 differ diff --git a/BandGraphs/ext/BandGraphsGraphMakieExt.jl b/BandGraphs/ext/BandGraphsGraphMakieExt.jl index f74de673..06ff5647 100644 --- a/BandGraphs/ext/BandGraphsGraphMakieExt.jl +++ b/BandGraphs/ext/BandGraphsGraphMakieExt.jl @@ -23,7 +23,7 @@ import BandGraphs: plot_flattened_bandgraph, make_vertices_dragable! # exported ## Unfolding a band graph `g` via `kg_trail` function plot_flattened_bandgraph( n::SymVector{D}, - lgirsd::Dict{String, Vector{LGIrrep{D}}}; + lgirsd::Dict{String, <:AbstractVector{LGIrrep{D}}}; timereversal=true ) where D # TODO: take a `SubductionTable` as input? @@ -66,7 +66,8 @@ function plot_flattened_bandgraph( # manifolds (lines/planes, etc) if isnothing(xys) force_equal_layers = find_equal_maximal_layers(klabs_trail, bandg.partitions) - xs, ys, _ = solve_positions(Zarate(), g_trail; force_equal_layers) + force_layer = assign_layer_indices(g_trail) + xs, ys, _ = solve_positions(Zarate(), g_trail; force_layer, force_equal_layers) ys = linearize_nonmaximal_y_coordinates(g_trail, ys) else xs, ys = xys @@ -75,6 +76,12 @@ function plot_flattened_bandgraph( plot_flattened_bandgraph(g_trail, klabs_trail, xs, ys) end +function assign_layer_indices(g_trail) + # the x-layer of the vertices should equal the trail index of the vertex + vs = vertices(g_trail) + return vs .=> [g_trail[label_for(g_trail, i)].trailidx for i in vs] +end + function plot_flattened_bandgraph( g_trail :: MetaGraph, # ::IsDirected klabs_trail :: AbstractVector{<:AbstractString}, @@ -158,6 +165,11 @@ function find_equal_maximal_layers(klabs_trail, partitions) # don't need to process k-label if already discovered i ∈ seen && continue + # if this is a transition point between disconnected components (`kidx == -1` which + # maps to `klabs_trail[trailidx] = klab = ""`), there's nothing to enforce and we + # just go to the next vertex + isempty(klab) && continue + # also, check if this is a maximal manifold: if not, continue pidx = something(findfirst(p->p.klab == klab, partitions)) partitions[pidx].maximal || continue @@ -260,7 +272,12 @@ function make_vertices_dragable_bandgraph!(ax, p, g_trail) # find related irreps, their vertex indices, and change their y-coordinates as well vert_id = label_for(g_trail, idx)[1:2] - related_vert_idxs = findall(l->l[1:2]==(vert_id), g_trail.vertex_labels) + irlab, irmul = vert_id + irlab_normalized = replace(irlab, "′"=>"") # to also move monodromy-related verts + related_vert_idxs = findall(g_trail.vertex_labels) do l + irlab′, irmul′ = l + irlab_normalized == replace(irlab′, "′"=>"") && irmul == irmul′ + end for idx′ in related_vert_idxs prev_r′ = p[:node_pos][][idx′] p[:node_pos][][idx′] = Makie.Point(prev_r′[1], mouse_r[2]) diff --git a/BandGraphs/src/BandGraphs.jl b/BandGraphs/src/BandGraphs.jl index d9c3d5e3..f83570d6 100644 --- a/BandGraphs/src/BandGraphs.jl +++ b/BandGraphs/src/BandGraphs.jl @@ -13,6 +13,8 @@ using BlockArrays: BlockArray, Block # for `assemble_adjacency` # EXPORTS export + Connection, + SubductionTable, SymVector, Partition, SubGraph, @@ -28,43 +30,76 @@ export split_nonmaximal_nodes, permute_subgraphs, subduction_tables, - plot_flattened_bandgraph + plot_flattened_bandgraph, + findall_separable_vertices # ---------------------------------------------------------------------------------------- # # TODO: improve sharing of types between `BandGraphs` & `build/crawl_and_write_bandpaths.jl` include("subduction-types.jl") +include("subduction-table.jl") # ---------------------------------------------------------------------------------------- # # DATA -const SUBDUCTIONSD_TR = Dict{Int, Vector{SubductionTable{3}}}() # `timereversal = true` -const SUBDUCTIONSD = Dict{Int, Vector{SubductionTable{3}}}() # `timereversal = false` +const SUBDUCTIONSD_TR_3D = Dict{Int, Vector{SubductionTable{3}}}() # `timereversal = true` +const SUBDUCTIONSD_3D = Dict{Int, Vector{SubductionTable{3}}}() # `timereversal = false` +const SUBDUCTIONSD_TR_2D = Dict{Int, Vector{SubductionTable{2}}}() # `timereversal = true` +const SUBDUCTIONSD_2D = Dict{Int, Vector{SubductionTable{2}}}() # `timereversal = false` function __init__() # TODO: Fix the awfulness here: the loading below only works if the dataset is created # with the types from BandGraphs first - not a good circular thing to have... datapath = joinpath(dirname(@__DIR__), "data/connections/3d/subductions-tr.jld2") JLD2.jldopen(datapath, "r") do jldfile for (k,v) in jldfile["subductionsd"] - SUBDUCTIONSD_TR[k] = v + SUBDUCTIONSD_TR_3D[k] = v end end datapath = joinpath(dirname(@__DIR__), "data/connections/3d/subductions.jld2") JLD2.jldopen(datapath, "r") do jldfile for (k,v) in jldfile["subductionsd"] - SUBDUCTIONSD[k] = v + SUBDUCTIONSD_3D[k] = v + end + end + datapath = joinpath(dirname(@__DIR__), "data/connections/2d/subductions-tr.jld2") + JLD2.jldopen(datapath, "r") do jldfile + for (k,v) in jldfile["subductionsd"] + SUBDUCTIONSD_TR_2D[k] = v + end + end + datapath = joinpath(dirname(@__DIR__), "data/connections/2d/subductions.jld2") + JLD2.jldopen(datapath, "r") do jldfile + for (k,v) in jldfile["subductionsd"] + SUBDUCTIONSD_2D[k] = v end end end """ - subduction_tables(sgnum; timereversal=true) --> Vector{SubductionTable{3}} + subduction_tables(sgnum, ::Val{D}; timereversal=true) --> Vector{SubductionTable{D}} Return a vector of `SubductionTable`s from stored tabulations, with (`timereversal = true`) or without (`timereversal = false`) time-reversal symmetry in space group `sgnum`. """ -function subduction_tables(sgnum; timereversal::Bool=true) - (timereversal ? SUBDUCTIONSD_TR[sgnum] : SUBDUCTIONSD[sgnum])::Vector{SubductionTable{3}} +function subduction_tables(sgnum::Integer, ::Val{D}=Val(3); timereversal::Bool=true) where D + if D == 3 + if timereversal + return SUBDUCTIONSD_TR_3D[sgnum]::Vector{SubductionTable{3}} + else + return SUBDUCTIONSD_3D[sgnum]::Vector{SubductionTable{3}} + end + elseif D == 2 + if timereversal + return SUBDUCTIONSD_TR_2D[sgnum]::Vector{SubductionTable{2}} + else + return SUBDUCTIONSD_2D[sgnum]::Vector{SubductionTable{2}} + end + else + throw(DomainError(D, "dimension must be 1, 2, or 3")) + end +end +function subduction_tables(sgnum::Integer, D::Integer; timereversal::Bool=true) + return subduction_tables(sgnum, Val(D); timereversal) end # ---------------------------------------------------------------------------------------- # @@ -76,6 +111,9 @@ include("graphs.jl") include("subgraph-permutations.jl") include("graph_routing.jl") include("unfold.jl") +include("complete-split.jl") +include("subsetsum.jl") +include("separable.jl") # ---------------------------------------------------------------------------------------- # # EXTENSIONS: loaded via Requires.jl on Julia versions klab * "′") + irmul′ = irmul + tmp′ = find_vertex_in_partitions(partitions, irlab′, irmul′) + if !isnothing(tmp′) + kidx′, iridx′, midx′ = tmp′ + # we must also split up the monodromy-related vertex - and in general, we can + # just explore all subgraphs for the affected subgraphs as well + splits_bandg = _complete_split(splits_bandg, irmul′, kidx′, iridx′, midx′, irdim_v, + #=monodromy=# true) + end + + return splits_bandg +end +function complete_split(bandg :: BandGraph{D}, v::Tuple{LGIrrep{D}, Int}) where D + lgir, irmul = v + return complete_split(bandg, (label(lgir), irmul)) +end + +function _complete_split( + bandg :: Union{BandGraph{D}, BandGraphPermutations{D}}, + irmul :: Int, + kidx :: Int, + iridx :: Int, + midx :: Int, + irdim :: Int, + monodromy :: Bool = false + ) where D + + partitions = bandg.partitions + + # modify the specific partition that `v = (kidx, iridx, midx)` belongs to by removing + # it and adding two new "split-off" surrogate vertices in its place + split_p = build_split_partition(partitions, kidx, iridx, midx, irdim) + + # all subsequent partitions in `bandg` have their `iridxs` incremented by 1 since the + # degree of the graph is now one larger + split_partitions = Vector{Partition{D}}(undef, length(partitions)) + for (kidx′, p′) in enumerate(partitions) + if kidx′ < kidx + split_partitions[kidx′] = p′ + elseif kidx′ == kidx + split_partitions[kidx′] = split_p + else # kidx′ > kidx + split_partitions[kidx′] = Partition(p′.klab, p′.lgirs, p′.multiples, p′.maximal, + p′.kidx, p′.iridxs .+ (irdim-1)) + end + end + + # `split_partitions` is now the new "store" of vertex data for the graph; but the actual + # graph structure is specified in `subgraphs`, so this must be updated accordingly; + # this is the real work of the splitting - and where we must create the different kinds + # of admissible vertex splits by reconnecting the edges in a meaningful and + # comprehensive way + subgraphs = bandg isa BandGraphPermutations ? bandg.subgraphs_ps : bandg.subgraphs + + # first, find the parts of `subgraphs` that are affected by the split + affected_subgraphs_idxs = findall(s -> s.p_max.kidx == kidx, subgraphs) + affected_subgraphs = subgraphs[affected_subgraphs_idxs] + + row_in_A = partitions[kidx].multiples[midx][irmul] + affected_subgraphs_As = [Matrix{Int}[] for _ in 1:length(affected_subgraphs)] + for (i, s) in enumerate(affected_subgraphs) + A = if bandg isa BandGraphPermutations + As = s.As + length(As) == 1 || error("cannot split a subgraph which already has multiple permutations") + As[1] + else + s.A + end + + cols_in_A = findall(!iszero, A[row_in_A,:]) + A_part = A[row_in_A, cols_in_A] + if !(length(A_part) == irdim && all(isone, A_part)) + # unexpected nonzero parts of subgraph A: expected `fill(1, irdim)`; vertex is + # not fully splittable at `v`; this likely indicates that `v` is connected to + # a nonmaximal irrep which is not 1-dimensional; return `nothing` as sentinel + return nothing + # FIXME: Consider if this really is meaningful behavior in the context of + # keyword argument `separable_degree` to `findall_separable_vertices`: + # it doesn't seem that meaningful - it's not clear that this is really + # even necessary: maybe the vertex should just not be split completely + # in this case. Or we should, synthetically, allow a block adjacency + # matrix like [1 0; 0 1; 0 1] for a case where `A_part` = [1,2] and + # we are splitting a 3-dimensional irrep. I wonder if that breaks any + # of our "sum rules" for band graphs + end + + # now, we replace A = [1,1,...] by a permutation of the identity matrix (e.g., for + # irdim=2, we replace A = [1,1] by either [1 0; 0 1] or [0 1; 1 0]); if this is the + # first of the affected subgraphs, we only include the identity matrix, since the + # inclusion of additional permutations would give graphs that are isomorphic under + # edge contraction with respect to the split-off surrogate vertices; if it's a + # monodromy-related vertex, we just permute in any case to be safe + for (q, perm) in enumerate(permutations(1:irdim)) + if !monodromy + (i == 1 && q ≠ 1) && continue # skip [0 1; 1 0] cf. above + end + if first(s.p_nonmax.klab) == 'Ω' + # we want specifically to retain a mapping relationship between monodromy- + # paired irreps that are tied together via Ω<...> since these edges signify + # that we physically require monodromy-paired irreps to have the same energy + q ≠ 1 && continue # skip [0 1; 1 0] + end + + top = A[1:row_in_A-1, :] + bottom = A[row_in_A+1:end, :] + A_part = I(irdim)[:, perm] # TODO: optimize to avoid explicit creation of A_part + vcat_to_bottom = zeros(Int, irdim, size(A, 2)) + vcat_to_bottom[:, cols_in_A] .= A_part + split_A = vcat(top, bottom, vcat_to_bottom) + + push!(affected_subgraphs_As[i], split_A) + end + end + + # build a new `split_subgraphs_ps`, using `split_partitions` & `affected_subgraphs_As` + i′ = 0 + split_subgraphs_ps = Vector{SubGraphPermutations{D}}(undef, length(subgraphs)) + for (i, s) in enumerate(subgraphs) + As = if i ∉ affected_subgraphs_idxs + bandg isa BandGraphPermutations ? s.As : [s.A] + else + affected_subgraphs_As[i′+=1] + end + p_max = split_partitions[s.p_max.kidx] + p_nonmax = split_partitions[s.p_nonmax.kidx] + split_s = SubGraphPermutations(p_max, p_nonmax, As, s.pinned) + split_subgraphs_ps[i] = split_s + end + + split_bandg = BandGraphPermutations(split_partitions, split_subgraphs_ps) + return split_bandg +end + +function find_vertex_in_partitions( + partitions :: Vector{<:Partition}, irlab :: AbstractString, irmul :: Int) + klab = klabel(irlab) + for (kidx, p) in enumerate(partitions) + klab == p.klab || continue + + idxs = findall(lgir -> label(lgir) == irlab, p.lgirs) + + isempty(idxs) && error(lazy"could not find irrep label $irlab at $klab ($p.lgirs)") + irmul > length(idxs) && error(lazy"irrep index of $irmul exceeds cardinality of irrep $irlab") + kidx == p.kidx || error("unexpected misalignment of partition `kidx` and its index in the `bandg.partitions` vector") + + iridx = idxs[irmul] + midx = something(findfirst(∋(iridx), p.multiples)) + + # return indices corresponding to `v = (irlab, irmul)` into: + # - kidx: partitions[kidx] + # - iridx: partitions[kidx].lgirs[iridx] + # - midx: partitions[kidx].multiples[midx] + return (kidx, iridx, midx) + end + return nothing +end + +function build_split_partition( + partitions :: AbstractVector{Partition{D}}, kidx, iridx, midx, irdim) where D + + # we want to build a new, augmented partition, where the vertex corresponding to + # `(kidx, iridx, midx)` is split-up into two new "split-surrogate" vertices; this + # new augmented partition differs from the original partition in `lgirs`, `multiples`, + # and `iridxs` + + # original partition + p₀ = partitions[kidx] + + # augmented `multiples` + multiples₀ = p₀.multiples + complete_removal = length(multiples₀[midx]) == 1 + # ↑ if `true`, there will be no more irreps of the removed kind in the partition after + # removal; if `false`, one or more irreps of the removed kind will remain in the + # partition, after removal + + split_multiples = Vector{UnitRange{Int}}(undef, length(multiples₀) + !complete_removal) + split_multiples[1:midx-1] = multiples₀[1:midx-1] # entries before `midx` + for i in (midx+!complete_removal):length(split_multiples)-1 # entries after `midx` + split_multiples[i] = multiples₀[i+complete_removal] .- 1 + end + # insert augmented-irrep multiples + split_multiples[end] = (last(multiples₀[end])) : (last(multiples₀[end]) + (irdim-1)) + if !complete_removal # retain some of removed multiples + split_multiples[midx] = first(multiples₀[midx]) : (last(multiples₀[midx]) - 1) + end + + # augmented irreps + lgirs₀ = p₀.lgirs + split_lgirs = [lgir for (i, lgir) in enumerate(lgirs₀) if i != iridx] + + # insert two dummy irreps corresponding to the split-off irreps; we give them "trivial" + # irrep matrices, since all we care about is their labels - for the labels, we give + # them the original label, appended with an "ˣ" suffix + lgir₀ = lgirs₀[iridx] + split_irlab = label(lgir₀) * 'ˣ' + split_lgir = LGIrrep( + split_irlab, + group(lgir₀), + [ones(ComplexF64, 1, 1) for _ in 1:length(group(lgir₀))], + [zeros(Float64, D) for _ in 1:length(group(lgir₀))], + REAL, false) + append!(split_lgirs, fill(split_lgir, irdim)) # insert `irdim` times + + # augmented indices into irreps + iridxs₀ = p₀.iridxs + split_iridxs = first(iridxs₀) : (last(iridxs₀) + (irdim-1)) + + return Partition(p₀.klab, split_lgirs, split_multiples, p₀.maximal, p₀.kidx, + split_iridxs) +end + + +# a split can be invalid if it involves a monodromy-paired vertex: in particular, if the +# split happens to disconnect the monodromy-paired vertices, then the split is invalid since +# these vertices must always be connected (i.e., live together in a connected component); +# equivalently, all monodromy-related vertices v and v′ must be connected by a path in the +# graph; this function checks that and is intended to be used to "filter out" invalid +# splits that might be generated by the combinatorial splitting in `complete_split` +function is_valid_split(bandg :: BandGraph) + partitions = bandg.partitions + any(p->last(p.klab) == '′', partitions) || return true + + g = assemble_simple_graph(bandg) + for p′ in Iterators.reverse(partitions) + klab′ = p′.klab + last(klab′) == '′' || continue + + klab = SubString(klab′, firstindex(klab′), prevind(klab′, lastindex(klab′))) # faster `chop(klab′)` + i = findfirst(p -> p.klab == klab, partitions) + if isnothing(i) + klabs = [p.klab for p in partitions] + error(lazy"could not find a monodromy-related partition for $klab′ (trying to match to $klab); available partitions are $klabs") + end + p = partitions[i] # "original" partition (not monodromy-translated) + + iridxs′ = p′.iridxs + iridxs = p.iridxs + + # now check that there is a path in `g` from `iridx` to `iridx′` + bool = all(zip(iridxs, iridxs′)) do (iridx, iridx′) + has_path(g, iridx, iridx′) + end + if !bool + return false + end + end + return true +end \ No newline at end of file diff --git a/BandGraphs/src/graph_routing.jl b/BandGraphs/src/graph_routing.jl index 033826fd..1639a078 100644 --- a/BandGraphs/src/graph_routing.jl +++ b/BandGraphs/src/graph_routing.jl @@ -26,7 +26,7 @@ function chinese_postman( check_connectedness::Bool=true ) where T # if disconnected, split into connected components and run on each component separately, - # concatenating trails and indicating disconnections between by -1 elements + # concatenating trails and indicating disconnections between by a sentinel value of `-1` if check_connectedness && !is_connected(g) return _chinese_postman_split_disconnected_components(g) end diff --git a/BandGraphs/src/graphs.jl b/BandGraphs/src/graphs.jl index 68a04ff4..33cb647d 100644 --- a/BandGraphs/src/graphs.jl +++ b/BandGraphs/src/graphs.jl @@ -39,6 +39,28 @@ function assemble_graph(subgraphs, partitions) return g end assemble_graph(bandg::BandGraph) = assemble_graph(bandg.subgraphs, bandg.partitions) + + +function assemble_simple_graph(subgraphs, partitions) + # does not have the right "weights" (and hence, not the right adjacency matrix either), + # since a simple graph must have unit weights, but can still be used for connectivity + # analysis for instance + g = Graph(last(partitions[end].iridxs)) + + # add edges to graph + for subgraph in subgraphs + max_iridxs, nonmax_iridxs = subgraph.p_max.iridxs, subgraph.p_nonmax.iridxs + for (j,a) in enumerate(eachcol(subgraph.A)) # j: local column index in block/subgraph + j′ = nonmax_iridxs[j] # global column index in overall graph + i = something(findfirst(≠(0), a)) # local row index in block/subgraph + i′ = max_iridxs[i] # global column index in overall graph + add_edge!(g, i′, j′) + end + end + return g +end +assemble_simple_graph(bandg::BandGraph) = assemble_simple_graph(bandg.subgraphs, bandg.partitions) + # ---------------------------------------------------------------------------------------- # function partition_graph(subgraphs, partitions) @@ -73,7 +95,7 @@ function split_nonmaximal_nodes(kg::MetaGraph) vertex_data_type = @NamedTuple{klab::String, code::Int, maximal::Bool}) for i in vertices(kg) klab = label_for(kg, i) - kg[klab].maximal || break # assume that all non-maximal k-points come in sequence + kg[klab].maximal || continue add_vertex!(kg′, (klab, 1), (; klab=klab, code=i, maximal=true)) end for i in vertices(kg) @@ -81,16 +103,22 @@ function split_nonmaximal_nodes(kg::MetaGraph) kg[klab_nonmax].maximal && continue neighbor_codes = neighbors(kg, i) neighbor_labels = label_for.(Ref(kg), neighbor_codes) - j = 0 - for (n, klabⁿ) in enumerate(neighbor_labels) - for m in n+1:length(neighbor_labels) - j += 1 - klabᵐ = neighbor_labels[m] - nonmax_vertex_label = (klab_nonmax, j) - add_vertex!(kg′, nonmax_vertex_label, - (; klab=klab_nonmax, code=i, maximal=false)) - add_edge!(kg′, (klabⁿ, 1), nonmax_vertex_label, nothing) - add_edge!(kg′, (klabᵐ, 1), nonmax_vertex_label, nothing) + if length(neighbor_labels) == 1 + nonmax_vertex_label = (klab_nonmax, 1) + add_vertex!(kg′, nonmax_vertex_label, (; klab=klab_nonmax, code=i, maximal=false)) + add_edge!(kg′, (neighbor_labels[1], 1), nonmax_vertex_label, nothing) + else + j = 0 + for (n, klabⁿ) in enumerate(neighbor_labels) + for m in n+1:length(neighbor_labels) + j += 1 + klabᵐ = neighbor_labels[m] + nonmax_vertex_label = (klab_nonmax, j) + add_vertex!(kg′, nonmax_vertex_label, + (; klab=klab_nonmax, code=i, maximal=false)) + add_edge!(kg′, (klabⁿ, 1), nonmax_vertex_label, nothing) + add_edge!(kg′, (klabᵐ, 1), nonmax_vertex_label, nothing) + end end end end diff --git a/BandGraphs/src/separable.jl b/BandGraphs/src/separable.jl new file mode 100644 index 00000000..963f9151 --- /dev/null +++ b/BandGraphs/src/separable.jl @@ -0,0 +1,370 @@ +""" +Find all the vertices in a band graph `bandg` whose associated irrep `lgir` for which +`criterion(lgir)` is true. + +Return a vector of 2-tuples, whose first element is a valid `lgir`, and whose second element +is a linear index into the multiplicity of this irrep in `bandg`. +""" +function findall_vertices(criterion, bandg :: BandGraph{D}) where D + vs = Vector{Tuple{LGIrrep{D}, Int}}() + for p in bandg.partitions + # only include non-monodromy partitions, since monodromy irreps will be split + # with their parents anyway; i.e., don't double count + is_monodromy = endswith(p.klab, "′") + is_monodromy && continue + + # only consider splitting at maximal vertices + p.maximal || continue + + for ms in p.multiples + for (i, m) in enumerate(ms) + lgir = p.lgirs[m] + if criterion(lgir) :: Bool + push!(vs, (lgir, i)) + end + end + end + end + return vs +end + +function findall_separable_vertices( + criterion, n::SymVector{D}, subts, lgirsd; kws...) where D + bandg = build_subgraphs(n, subts, lgirsd) + return findall_separable_vertices(criterion, bandg; kws...) +end +function findall_separable_vertices( + criterion, + bandg :: BandGraph{D}; + max_permutations::Union{Nothing, <:Real} = 1e6, + separable_degree :: Union{Nothing, <: Integer} = nothing) where D + + separable = Tuple{BandGraph{D}, BandGraph{D}, LGIrrep{D}}[] + + # find which vertices fulfil criterion; i.e., which vertices to check seperability for + vs = findall_vertices(lgir -> isspecial(lgir) && criterion(lgir), bandg) + isempty(vs) && return false, separable + + # fast-path "necessary conditions" check which does not require permutation construction + unfeasible_vs_idxs = fast_path_separability_irdim_infeasibility(bandg, vs, separable_degree) + deleteat!(vs, unfeasible_vs_idxs) + isempty(vs) && return false, separable + + # move on to additional conditions which require us to look at individual permutations + subgraphs_ps = permute_subgraphs(bandg.subgraphs) + bandgp = BandGraphPermutations(bandg.partitions, subgraphs_ps) + + # if there are too many permutations to consider, we bail out (returning `missing`) + if !isnothing(max_permutations) && length(bandgp) > max_permutations + printstyled("band permutations (", length(bandgp), + ") exceed `max_permutations=", max_permutations, "`; aborting ...\n"; + color=:yellow) + return missing + end + + # okay, now we actually have to do the "full" checks for any remaining vertices by + # explicitly considering each possible permutation explicitly + for v in vs + lgir, irmul = v + for bandg′ in bandgp # iterate over graph permutations + split_bandg = is_separable_at_vertex(bandg′, (lgir, irmul); + verbose=true, separable_degree) + if !isnothing(split_bandg) + # constrain ourselves to return only a single representative seperable band + # permutation if any such permutation exists + push!(separable, (bandg′, split_bandg, lgir)) + break + end + end + end + + return !isempty(separable), separable +end + +function fast_path_separability_irdim_infeasibility( + bandg::BandGraph{D}, + vs::Vector{Tuple{LGIrrep{D}, Int}}, + separable_degree :: Union{Nothing, <:Integer} = nothing + ) where D + # - for `v` to be splittable, it must be possible to separate the irreps of every + # partition into groups that contain a single split-off surrogate vertex `vˣ`, with + # each such group having a consistent number of bands; this is precisely the job + # of `solve_subset_sum_variant` to decide if this is possible (but there, across + # bands over subsets of partitions, not over individual partitions, as here). + # - this is equivalent to ignoring the additional constraints from compatibility + # relations; but since they are _additional_ this is a still a necessary condition. + # - if `v` has a monodromy partner, we simply remove that partition from the + # consideration, since the monodromy partner would be split up in exactly the same + # way + # - we return a list of indices `infeasible_idxs` into `vs`, with each index signifying + # a vertex that failed this test + + partitions_irdims = [irdim.(p.lgirs) for p in bandg.partitions] + infeasible_idxs = Int[] + for (i, v) in enumerate(vs) + lgir, irmul = v + irlab = label(lgir) + kidx, iridx, _ = something(find_vertex_in_partitions(bandg.partitions, irlab, irmul)) + v_irdim = partitions_irdims[kidx][iridx] + + # check if there are any monodromy-related partner to `v`; if so, don't include the + # monodromy-related partition in the set of "exterior" partitions + klab′ = klabel(irlab) * '′' + kidx′ = findfirst(p->p.klab == klab′, bandg.partitions) + + # now see if there is a possible irrep-grouping across partitions exists that allow + # for a single split-off surrogate vertex in each group; collect irrep dimensions + # "outside" the one belong to `v` in `Sʲs`, and those belong to `v` in `T` and `R` + # with `T` denoting the irrep dimensions of the split-off surrogate vertices `vˣ` + # and `R` denoting the remaining irreps in the partition where `v` lives + Sʲs = if isnothing(kidx′) # no monodromy partner + [Sʲ for (j, Sʲ) in enumerate(partitions_irdims) if j ≠ kidx] + else + [Sʲ for (j, Sʲ) in enumerate(partitions_irdims) if j ∉ (kidx, kidx′)] + end + T = fill(1, v_irdim) # TODO: maybe someday a more general thing based on non-maximal irrep multiplicities incident on `v`? + R = Int[v′_irdim for (iridx′, v′_irdim) in enumerate(partitions_irdims[kidx]) if iridx′≠iridx] + + feasible, _ = if isnothing(separable_degree) || separable_degree == irdim(lgir) + solve_subset_sum_variant(Sʲs, T, R) + else + solve_subset_sum_variant_flexibleT(Sʲs, T, R, separable_degree) + end + if !feasible + push!(infeasible_idxs, i) + end + end + + infeasible_idxs +end + +function is_separable_at_vertex( + bandg :: BandGraph{D}, + v :: Tuple{LGIrrep{D}, Int}; + verbose :: Bool = true, + separable_degree :: Union{Nothing, <:Integer} = nothing) where D + + # The function determines whether a graph G is separable at a vertex v, as described in + # our Overleaf document; before doing "heavy-lifting", we do a fast-path "necessary + # conditions" check, to potentially avoid actually creating any explicit band splits + # if the vertex is clearly not separable + + + # fast path checks + # TODO: articulation point check (must be articulation point); but this only works + # if `v` does not have a monodromy partner vertex + + # "slow-path" checks + splits_bandg = complete_split(bandg, v) + isnothing(splits_bandg) && return nothing # no valid splits at `v` + + # set how many components the band graph should be separated into after the split, in + # order for us to classify it is a "succesful" split + split_lgir = first(v) + split_irdim = dim(split_lgir) + if isnothing(separable_degree) + # "full" separability: into as many distinct band graphs as the split-irrep's dimen. + separable_degree = split_irdim + else + # "custom": any number less than or equal to `split_irdim` + if split_irdim < separable_degree + error(lazy"`separable_degree` (=$separable_degree) cannot exceed the splitted irrep's dimension (=$split_irdim)") + end + end + + # obtain the connected components of the split-up graph, across both partitions and + # energy-separation: `bandg_cs[partitions-grouping][energy-grouping]` + bandg_cs = group_partition_connected_components(bandg) # components + + if isnothing(bandg_cs) + # CASE 1: original band graph is fully connected - easy to check separability + # (all distinct components after the split will cover the same set of + # partitions) + for split_bandg in splits_bandg # iterate over split permutations + if !is_valid_split(split_bandg) + verbose && printstyled(" ... skipping an invalid split\n"; color=:yellow) + continue + end + split_g = assemble_graph(split_bandg) + split_N = length(connected_components(split_g)) + if split_N == separable_degree + return split_bandg + end + end + else + # CASE 2: original band graph is not fully connected - harder to check separability + # (the distinct components may involve different partitions) + split_surrogate_irlab = label(split_lgir) * "ˣ" # label for irrep after split + for split_bandg in splits_bandg # iterate over split permutations + if !is_valid_split(split_bandg) + verbose && printstyled(" ... skipping an invalid split\n"; color=:yellow) + continue + end + split_bandg_cs = group_partition_connected_components(split_bandg) + ijs = findall_vertices_in_components(split_bandg_cs, split_surrogate_irlab) + + # after split, there must be at least `separable_degree` components that + # each contain a "split surrogate" vertex; if not, no need to look further + length(ijs) == separable_degree || continue + + # now we must check if it is possible to group the "split-off" components into + # isolated band graphs, each with a consistent occupation number + i = first(first(ijs)) # index into `partition-group` in `split_bandg_cs` + @assert all(ij -> ij[1] == i, ijs) + T = [occupation(split_bandg_cs[i][j]) for (i,j) in ijs] + i_components = split_bandg_cs[i] + js = getindex.(ijs, 2) + R = [occupation(i_components[j′]) for j′ in 1:length(i_components) if j′ ∉ js] + Sʲs = [[occupation(split_bandg_c) for split_bandg_c in split_bandg_cs[i′]] + for i′ in 1:length(split_bandg_cs) if i′≠i] + feasible, _ = solve_subset_sum_variant(Sʲs, T, R; verbose=false) + if feasible + # it's feasible to group the split-off vertices into their own distinct and + # isolated band graphs + return split_bandg + end + end + + end + + return nothing +end + +function Graphs.induced_subgraph(bandg :: BandGraph{D}, vlist) where D + partitions = bandg.partitions + issorted(vlist) || (vlist = sort(vlist)) + + # determine which partitions are correspond to the vertex indices in `vlist` + # (this aligns with the vertex-indexing ordering from `assemble_graph`) + kidx = 1 + irmax = 0 + code = 0 + partitions′ = empty(partitions) + vmap = Int[] + lgirs_idxs = Int[] + for p in partitions + for j in eachindex(p.lgirs) + code += 1 + if insorted(code, vlist) + push!(lgirs_idxs, j) + push!(vmap, code) + end + end + isempty(lgirs_idxs) && continue + lgirs = p.lgirs[lgirs_idxs] + multiples = _identify_multiples(lgirs, first(lgirs)) + p′ = Partition{D}(p.klab, lgirs, multiples, p.maximal, kidx, irmax+1:irmax+length(lgirs)) + push!(partitions′, p′) + + kidx += 1 + irmax = last(p′.iridxs) + empty!(lgirs_idxs) + end + + # determine the subgraphs connected by elements in vlist (exploit that this is a + # bipartite graph) + subgraphs = bandg.subgraphs + subgraphs′ = empty(subgraphs) + for subgraph in subgraphs + klab_max = subgraph.p_max.klab + klab_nonmax = subgraph.p_nonmax.klab + idx_max = findfirst(p-> p.klab==klab_max, partitions′) + isnothing(idx_max) && continue + idx_nonmax = findfirst(p -> p.klab==klab_nonmax, partitions′) + if isnothing(idx_nonmax) + error("could not find expected nonmaximal partition from k-point $klab_nonmax") + end + + # determine include parts of the subgraph adjacency matrix + p_max′ = partitions′[idx_max] + p_nonmax′ = partitions′[idx_nonmax] + A = subgraph.A + A′ = Matrix{Int}(undef, length(p_max′.lgirs), length(p_nonmax′.lgirs)) + for (local_idx_max′, iridx_max′) in enumerate(p_max′.iridxs) + for (local_idx_nonmax′, iridx_nonmax′) in enumerate(p_nonmax′.iridxs) + iridx_max = vmap[iridx_max′] + iridx_nonmax = vmap[iridx_nonmax′] + local_idx_max = something(findfirst(==(iridx_max), subgraph.p_max.iridxs)) + local_idx_nonmax = something(findfirst(==(iridx_nonmax), subgraph.p_nonmax.iridxs)) + A′[local_idx_max′, local_idx_nonmax′] = A[local_idx_max, local_idx_nonmax] + end + end + subgraph′ = SubGraph{D}(p_max′, p_nonmax′, A′, subgraph.pinned) + push!(subgraphs′, subgraph′) + end + + bandg′ = BandGraph(subgraphs′, partitions′) + return bandg′, vmap +end + +function _identify_multiples!(multiples :: Vector{UnitRange{Int}}, lgirs, lgir_target, start) + # we assume that `lgirs` is sorted in the sense that identical `lgirs` are adjacent + stop = start + for j in start+1:length(lgirs) + lgir = lgirs[j] + if label(lgir) == label(lgir_target) + stop = j + else + push!(multiples, start:stop) + return _identify_multiples!(multiples, lgirs, lgirs[j], j) + end + end + push!(multiples, start:stop) + return multiples +end +function _identify_multiples(lgirs, lgir_target, start::Int=1) + _identify_multiples!(UnitRange{Int}[], lgirs, lgir_target, start) +end + +# return the connected components of `bandg`; returned components are a +# `Vector{Vector{BandGraph}}` with indexing `bandg_cs[partitions-grouping][energy-grouping]` +# if there is only one connected component (fully connected), returns `nothing` as sentinel +function group_partition_connected_components(bandg :: BandGraph{D}) where D + g = assemble_simple_graph(bandg) + cs = connected_components(g) + if length(cs) == 1 + # `nothing` as sentinel for being fully connected (i.e., 1 connected component) + return nothing + end + + klabs_cs = Set{String}[] + bandg_cs = Vector{BandGraph{D}}[] + for c in cs + bandg_c, vmap = induced_subgraph(bandg, c) + klabs_c = Set(p.klab for p in bandg_c.partitions) + idx = findfirst(==(klabs_c), klabs_cs) + if isnothing(idx) + push!(klabs_cs, klabs_c) + push!(bandg_cs, [bandg_c]) + else + push!(bandg_cs[idx], bandg_c) + end + end + + return bandg_cs +end + +function findall_vertices_in_components( + bandg_cs :: Vector{Vector{BandGraph{D}}}, irlab :: String) where D + klab = klabel(irlab) + contains_vertices_idxs = Tuple{Int, Int}[] + for (i, bandgs) in enumerate(bandg_cs) # over each k-connected partition + for (j, bandg_c) in enumerate(bandgs) # over (energy-separated) components of a k-connected partition + is_matching_partition_group = false + partitions_c = bandg_c.partitions + for partition_c in partitions_c + if partition_c.klab == klab + is_matching_partition_group = true + if any(lgir->label(lgir)==(irlab), partition_c.lgirs) + push!(contains_vertices_idxs, (i,j)) + end + break # found what we looked for; go to next graph + else + continue # keep looking + end + end + is_matching_partition_group || break # wrong partitions; skip related graphs + end + end + return contains_vertices_idxs +end \ No newline at end of file diff --git a/BandGraphs/src/subduction-table.jl b/BandGraphs/src/subduction-table.jl index dba0e517..a4ac6bb4 100644 --- a/BandGraphs/src/subduction-table.jl +++ b/BandGraphs/src/subduction-table.jl @@ -9,15 +9,50 @@ function SubductionTable( lgirsᴳ = lgirsd[bare_kᴳ_label] # determine the free parameters for which kᴳ and kᴴ are compatible/"intersect" + # NB: `Gp` below is a primitive reciprocal lattice vector s.t. kᴴₚ(αβγ′) + Gp = kᴳₚ + # with kᴴₚ and kᴳₚ denoting the primitive settings of kᴴ and kᴳ, respectively. cntr = centering(sg) - bool, αβγ, αβγ′, G = Crystalline.can_intersect(primitivize(kᴳ.kv, cntr), - primitivize(kᴴ.kv, cntr)) + bool, αβγ, αβγ′, Gp = Crystalline.can_intersect(primitivize(kᴳ.kv, cntr), + primitivize(kᴴ.kv, cntr)) bool || error("failed to prove compatibility of $(kᴳ) and $(kᴴ)") iszero(αβγ) || error("implementation cannot handle nonzero free parameters for G-group") - lgirsᴴ = lgirsd[string(kᴴ.label)] - cosetsᴴ = cosets(sg, group(first(lgirsᴴ))) - lgirsᴴ = Crystalline.remap_lgirreps_to_point_in_kstar(lgirsᴴ, kᴴ.kv, cosetsᴴ) + # `Gp` (primitive setting here) could in principle contain parts that are actually + # along `kᴴ`, which would give a different value of αβγ′: we want the reciprocal vector + # to have no parts along `kᴴ` (i.e., the reciprocal vector should lie in a line/plane + # orthogonal to the plane/line spanned by `kᴴ`; we fix it up here. + # [a motivating example is space group 122 (M ↓ Λ)] + if !iszero(Gp) + kᴴₚ = primitivize(kᴴ.kv, cntr) + Q = pseudo_inverse_smith(Matrix(kᴴₚ.free)) * Gp + if norm(Q) > Crystalline.DEFAULT_ATOL + ΔGp = kᴴₚ.free * Q + if all(v -> abs(v - round(v)) < Crystalline.DEFAULT_ATOL, ΔGp) + # we can & should change `Gp` by a reciprocal lattice vector, if, in fact, + # such a change would make `Gp` orthogonal to `kᴴₚ` + Gp_tentative = Gp - ΔGp + αβγ′_tentative = αβγ′ + Q + if norm(kᴴₚ(αβγ′_tentative) ⋅ Gp_tentative) < Crystalline.DEFAULT_ATOL + Gp = Gp_tentative + αβγ′ = αβγ′_tentative + else + # give up ¯\_(ツ)_/¯ + # can e.g., happen in space group 220 for H ↓ Λ = [1,1,1] ↓ [-α, α, -α] + # but apparently does not change the subduction table, so giving up + # seems okay, although undeniably iffy - the whole approach here is iffy + end + else + error("there be dragons here!") + # there be dragons... the `Gp` vector is not actually orthogonal to `kᴴ` + # ought to be verified that it cannot change the irrep regardless + end + end + # NB: we currently don't use `Gp` for anything - the idea is we shouldn't need to + end + + lgirsᴴ_tmp = lgirsd[string(kᴴ.label)] + cosetsᴴ = cosets(sg, group(first(lgirsᴴ_tmp))) + lgirsᴴ = Crystalline.remap_lgirreps_to_point_in_kstar(lgirsᴴ_tmp, kᴴ.kv, cosetsᴴ) # compute subduction table of G- and H-group irreps at their αβ′-intersection table = [subduction_count(lgirᴳ, lgirᴴ, αβγ′) for lgirᴳ in lgirsᴳ, lgirᴴ in lgirsᴴ] @@ -43,4 +78,12 @@ function SubductionTable( end return SubductionTable{D}(num(sg), c, irlabsᴳ, irlabsᴴ, table, monodromy) -end \ No newline at end of file +end + +function pseudo_inverse_smith(A) + # compute the pseudo-inverse of a matrix `A` using the Smith normal form + # (https://en.wikipedia.org/wiki/Pseudoinverse#Using_the_Smith_normal_form) + F = Crystalline.smith(A) + A⁺ = F.Tinv * pinv(diagm(F)) * F.Sinv + return A⁺ +end diff --git a/BandGraphs/src/subgraph-permutations.jl b/BandGraphs/src/subgraph-permutations.jl index 7b312619..c9045c81 100644 --- a/BandGraphs/src/subgraph-permutations.jl +++ b/BandGraphs/src/subgraph-permutations.jl @@ -61,13 +61,35 @@ function subgraph_permutations(subgraph) # for each set of same-irrep nodes, we can perform all possible permutations of the # columns in the adjacency matrix associated with those nodes; if there are multiple # same-irrep node sets, we must perform all permutations across all sets - multiples = subgraph.p_nonmax.multiples - pss = collect.(permutations.(multiples)) - Np = prod(length, pss) # aggregate number of same-irrep permutations, across irreps + multiples = subgraph.p_nonmax.multiples + pss = permutations.(multiples) + + # if the permuted columns of the subgraph adjacency matrix are identical, however, there + # is no point in including them, since the permutation would leave the associated + # subgraph adjacency matrix unchanged; we filter them out below + pss′ = [Vector{Int}[] for same_ir_idx in 1:length(pss)] + for (same_ir_idx, ps) in enumerate(pss) + As′_part = Matrix{Int}[] + ps′ = pss′[same_ir_idx] + for p in ps + A′_part = subgraph.A[:,p] + A′_part ∈ As′_part && continue + push!(As′_part, A′_part) + push!(ps′, p) + end + end + + # finally, collect the distinct subgraph adjacency matrices for these permutations + Np = prod(length, pss′) # aggregate number of same-irrep permutations, across irreps + if Np > 1_000_000 + # TODO: we probably should not be materializing these band permutations but should + # store this lazily via a permutation structure + error(lazy"impossibly many subgraph permutations ($Np); bailing out..!") + end As = Vector{Matrix{Int}}(undef, Np) - for (j,is) in enumerate(CartesianIndices(Tuple(Base.OneTo.(length.(pss))))) - # TODO: do this type-stably and without `collect`ing `pss` - cols = reduce(vcat, [pss[k][i] for (k,i) in enumerate(Tuple(is))]) + for (j,is) in enumerate(CartesianIndices(Tuple(Base.OneTo.(length.(pss′))))) + # TODO: do this type-stably + cols = reduce(vcat, [pss′[k][i] for (k,i) in enumerate(Tuple(is))]) As[j] = subgraph.A[:,cols] end return SubGraphPermutations(subgraph.p_max, subgraph.p_nonmax, As, #=pinned=# false) diff --git a/BandGraphs/src/subgraphs.jl b/BandGraphs/src/subgraphs.jl index 4a6d5228..515efaec 100644 --- a/BandGraphs/src/subgraphs.jl +++ b/BandGraphs/src/subgraphs.jl @@ -1,7 +1,7 @@ function build_subgraphs( n::SymVector{D}, subductions, - lgirsd::Dict{String, Vector{LGIrrep{D}}} + lgirsd::Dict{String, <:AbstractVector{LGIrrep{D}}} ) where D seen_irs = 0 # current number of irreps aggregated across all partitions @@ -18,6 +18,7 @@ function build_subgraphs( # monodromy additions (same irrep population as their "non-monodromy"/"original" irrep # variants, but connect differently to non-maximal k-manifolds); PRE 96, 023310 (2017) + monodromy_pairs = Tuple{Int, Int}[] for s in subductions s.monodromy || continue klab_max′ = string(s.c.kᴳ.label) # monodromy k-label (with a '′') @@ -47,6 +48,9 @@ function build_subgraphs( p_max = build_partition(klab_max′, n.mults[idx], lgirs_max′, kidx, seen_irs) push!(partitions_max, p_max) seen_irs = last(p_max.iridxs) + + kidx_original = something(findfirst(p->p.klab==klab_max, partitions_max)) + push!(monodromy_pairs, (kidx_original, kidx)) end # non-maximal k-point partitions @@ -103,7 +107,7 @@ function build_subgraphs( push!(partitions_nonmax, p_nonmax) end end - + # TODO: sort `partitions_max` and `partitions_nonmax` so that we benefit maximally from # from the pinning order when we consider permutations: basically, want to sort # in such a way that subgraphs with most permutations get pinned @@ -111,11 +115,11 @@ function build_subgraphs( # build all the subgraphs that make up the full graph subgraphs = SubGraph{D}[] for p_max in partitions_max # maximal k-manifolds - i = p_max.kidx + # i = p_max.kidx klab_max = p_max.klab lgirs_max = p_max.lgirs for p_nonmax in partitions_nonmax # non-maximal k-manifolds - j = p_nonmax.kidx + # j = p_nonmax.kidx klab_nonmax = p_nonmax.klab lgirs_nonmax = p_nonmax.lgirs @@ -148,12 +152,50 @@ function build_subgraphs( end end + # if we had any monodromy additions, we also add a non-maximal partition associated + # with the generic point Ω and its trivial irrep Ω₁, in order to tie together every + # pair of monodromy-related irreps via Ω₁: we do this to ensure that band graphs + # which are connected in energy are also connected as graphs (via Ω₁, at least) + if !isempty(monodromy_pairs) + lgir_Ω = only(lgirsd["Ω"]) + lg_Ω = group(lgir_Ω) + for (original_kidx, monodromy_kidx) in monodromy_pairs + original_p_max = partitions_max[original_kidx] + monodromy_p_max = partitions_max[monodromy_kidx] + + # prepare to add new Ω-partition to graph; but one that is specific to the + # connection between the two monodromy points + klab_Ω′ = "Ω" * klabel(original_p_max.klab) + irlab_Ω′ = klab_Ω′ * "₁" + + kidx += 1 + global_iridxs = seen_irs .+ (1:n.μ) + seen_irs = last(global_iridxs) + + lg_Ω′ = LittleGroup{D}(lg_Ω.num, lg_Ω.kv, klab_Ω′, lg_Ω.operations) + lgir_Ω′ = LGIrrep{D}(irlab_Ω′, lg_Ω′, lgir_Ω.matrices, lgir_Ω.translations, + lgir_Ω.reality, lgir_Ω.iscorep) + p_Ω′ = Partition(klab_Ω′, fill(lgir_Ω′, n.μ), [1:n.μ], false, kidx, global_iridxs) + push!(partitions_nonmax, p_Ω′) + + A = zeros(Int, length(original_p_max.lgirs), n.μ) + col_start = 1 + for (row, lgir) in enumerate(original_p_max.lgirs) + cols = col_start:(col_start+irdim(lgir)-1) + A[row,cols] .= 1 + col_start = last(cols)+1 + end + push!(subgraphs, SubGraph(original_p_max, p_Ω′, A, #=pinned=# true)) + push!(subgraphs, SubGraph(monodromy_p_max, p_Ω′, A, #=pinned=# true)) + end + end + partitions = vcat(partitions_max, partitions_nonmax) return BandGraph(subgraphs, partitions) end -function build_partition(klab, nsₖ, lgirsₖ::Vector{LGIrrep{D}}, kidx, seen_irs) where D + function build_partition(klab, nsₖ, lgirsₖ :: AbstractVector{LGIrrep{D}}, kidx, seen_irs) where D lgirs = LGIrrep{D}[] multiples = UnitRange{Int}[] max_local_iridx = 0 diff --git a/BandGraphs/src/subsetsum.jl b/BandGraphs/src/subsetsum.jl new file mode 100644 index 00000000..da2c36b3 --- /dev/null +++ b/BandGraphs/src/subsetsum.jl @@ -0,0 +1,253 @@ +using JuMP, HiGHS + +""" + solve_subset_sum_variant(Sʲs, T, R; verbose::Bool = false, optimizer) + +Solves the following variant of the subset sum problem, using binary integer programming. + +## Problem formulation +Given a target list `T` = (T₁, …, T_N), a set of multisets `Sʲs` = {S⁽ʲ⁾}_{j=1}^{J} with +S⁽ʲ⁾ = {s⁽ʲ⁾₁, …, s⁽ʲ⁾_{M⁽ʲ⁾}}, as well as a set `R` = {r₁, …, r_L}, determine if there +exist disjoint subset selections {S⁽ʲ⁾ₙ} and {Rₙ} with S⁽ʲ⁾ₙ ⊂ S⁽ʲ⁾ and Rₙ ⊂ R such that: + + ∑_{s ∈ S⁽ʲ⁾ₙ} s = tₙ + ∑_{r ∈ Rₙ} r + +for all j = 1, …, J and n = 1, …, N and with ⋃ₙ S⁽ʲ⁾ₙ = S⁽ʲ⁾} and ⋃ₙ Rₙ = R (and +S⁽ʲ⁾_n ∩ S⁽ʲ⁾_{n′} = ∅ if n ≠ n′ and similarly for R_n ∩ R_{n′}). + +## Output +- First return value (`::Bool`): whether the problem is feasible (i.e., solvable). +- The model (`::Model`), which contains the solution in binary variables `x` (into `Sʲs`) + and `y` (into `R`). + +## Example +```jl +julia> Sʲs = [[1, 3, 2, 1], [3, 2, 2], [1, 1, 2, 2, 1]]; +julia> T = [3, 3]; +julia> R = [1]; +julia> solve_subset_sum_variant(Sʲs, T, R; verbose=true); +Feasible +T = [3, 3] +R = [1] +Sʲs = [[1, 3, 2, 1], [3, 2, 2], [1, 1, 2, 2, 1]] + +j = 1: + n = 1: S[1] + S[2] = T[1] + (R[1]) + 1 + 3 = 3 + (1) + n = 2: S[3] + S[4] = T[2] + 2 + 1 = 3 +j = 2: + n = 1: S[2] + S[3] = T[1] + (R[1]) + 2 + 2 = 3 + (1) + n = 2: S[1] = T[2] + 3 = 3 +j = 3: + n = 1: S[1] + S[2] + S[3] = T[1] + (R[1]) + 1 + 1 + 2 = 3 + (1) + n = 2: S[4] + S[5] = T[2] + 2 + 1 = 3 +``` +""" +function solve_subset_sum_variant( + Sʲs :: AbstractVector{<:AbstractVector{<:Integer}}, + T :: AbstractVector{<:Integer}, + R :: AbstractVector{<:Integer}; + verbose :: Bool = false, + optimizer = optimizer_with_attributes(HiGHS.Optimizer, "presolve"=>"off") + # (HiGHS' presolve is buggy for this problem; disable it) + ) + + J = length(Sʲs) + N, Mʲs, L = length(T), length.(Sʲs), length(R) + + model = Model(optimizer)::Model + verbose || MOI.set(model, MOI.Silent(), true) + + @variable(model, x[1:sum(Mʲs), 1:N], Bin) # binary "indexing" into `Sʲs` + if L ≠ 0 + @variable(model, y[1:L, 1:N], Bin) # binary "indexing" into `R` + end + + # constraints on `x` and `y` + for n in 1:N + RHS = L ≠ 0 ? T[n] + sum(y[l,n]*R[l] for l in 1:L) : AffExpr(T[n]) + for j in 1:J + mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j]) + @constraint(model, sum(x[mᵍ,n]*Sʲs[j][m] + for (m,mᵍ) in enumerate(mʲ_global_idxs)) == RHS) + end + end + for j in 1:J + mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j]) + @constraint(model, sum(x[mᵍ,n] for mᵍ in mʲ_global_idxs for n in 1:N) == Mʲs[j]) + for mᵍ in mʲ_global_idxs + @constraint(model, sum(x[mᵍ,n] for n in 1:N) == 1) + end + end + + # constraints only on `y` + if L ≠ 0 + @constraint(model, sum(y[l,n] for l in 1:L for n in 1:N) == L) + for l in 1:L + @constraint(model, sum(y[l,n] for n in 1:N) == 1) + end + end + + optimize!(model) + + feasible = is_solved_and_feasible(model) + if !(feasible || termination_status(model) == INFEASIBLE) + error(lazy"unexpected termination status: $(termination_status(m))") + end + + verbose && _print_subset_sum_result(model, Sʲs, T, R) + + return feasible, model +end + +function _print_subset_sum_result(io :: IO, model :: Model, Sʲs, T, R) + J = length(Sʲs) + N, Mʲs, L = length(T), length.(Sʲs), length(R) + + if is_solved_and_feasible(model) + printstyled(io, "Feasible\n"; color=:green) + printstyled(io, "T = ", T, "\nR = ", R, "\nSʲs = ", Sʲs, "\n\n"; + color=:light_black) + x = object_dictionary(model)[:x] # into `Sʲs` + xv = value.(x) + yv = if L ≠ 0 + y = object_dictionary(model)[:y] # into `R` + value.(y) + else + Matrix{Float64}(undef, 0, 0) + end + for j in 1:J + mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j]) + printstyled(io, "j = ", j, ":\n", color=:light_black) + for n in 1:N + printstyled(io, " n = ", n, ": ", color=:light_black) + s = "" + v = "" + for (m, mᵍ) in enumerate(mʲ_global_idxs) + if xv[mᵍ,n] > 0.5 + s = s * (isempty(s) ? "" : " + ") * "S[$m]" + v = v * (isempty(v) ? "" : " + ") * string(Sʲs[j][m]) + end + end + s = *(s, " = T[", string(n), "]") + v = *(v, " = ", string(T[n])) + if L ≠ 0 + for l in 1:L + first = true + if yv[l,n] > 0.5 + if first + s *= " + (" + v *= " + (" + first = false + else + s *= " + " + v *= " + " + end + s *= "R[$l]" + v *= string(R[l]) + end + first || (s *= ")"; v *= ")") + end + end + println(io, s, "\n ", v) + end + end + elseif termination_status(model) == INFEASIBLE + printstyled(io, "Infeasible\n"; color=:red) + end +end +function _print_subset_sum_result(model :: Model, Sʲs, T, R) + _print_subset_sum_result(stdout, model, Sʲs, T, R) +end + +""" + solve_subset_sum_variant_flexibleT(Sʲs, T, R, Nᵀ; verbose::Bool = false, optimizer) + +Akin to [`solve_subset_sum_variant`](@ref), but now allowing a subset selection of `T` also, +in `Nᵀ` sets, rather than N = |T| sets. + +## Problem formulation +Given a target list `T` = (T₁, …, T_N), a set of multisets `Sʲs` = {S⁽ʲ⁾}_{j=1}^{J} with +S⁽ʲ⁾ = {s⁽ʲ⁾₁, …, s⁽ʲ⁾_{M⁽ʲ⁾}}, as well as a set `R` = {r₁, …, r_L}, determine if there +exist disjoint subset selections {S⁽ʲ⁾ₙ}, {Tₙ}, and {Rₙ} with S⁽ʲ⁾ₙ ⊂ S⁽ʲ⁾, Tₙ ⊂ T, and +Rₙ ⊂ R such that: + + ∑_{s ∈ S⁽ʲ⁾ₙ} s = ∑_{t ∈ Tₙ} t + ∑_{r ∈ Rₙ} r + +for n = 1, …, Nᵀ (with Nᵀ≤N) (and otherwise similar to [`solve_subset_sum_variant`](@ref)). +""" +function solve_subset_sum_variant_flexibleT( + Sʲs :: AbstractVector{<:AbstractVector{<:Integer}}, + T :: AbstractVector{<:Integer}, + R :: AbstractVector{<:Integer}, + Nᵀ :: Int = length(T); + verbose :: Bool = false, + optimizer = optimizer_with_attributes(HiGHS.Optimizer, "presolve"=>"off") + # (HiGHS' presolve is buggy for this problem; disable it) + ) + + J = length(Sʲs) + N, Mʲs, L = length(T), length.(Sʲs), length(R) + + model = Model(optimizer)::Model + verbose || MOI.set(model, MOI.Silent(), true) + + @variable(model, x[1:sum(Mʲs), 1:Nᵀ], Bin) # binary "indexing" into `Sʲs` + @variable(model, z[1:N, 1:Nᵀ], Bin) # binary "indexing" into `T` + if L ≠ 0 + @variable(model, y[1:L, 1:Nᵀ], Bin) # binary "indexing" into `R` + end + + # TODO: TEST + # constraints on `x` and `y` and `z` + for n in 1:Nᵀ + RHS_T = sum(z[m,n]*T[m] for m in 1:N) + RHS_R = L ≠ 0 ? sum(y[l,n]*R[l] for l in 1:L) : AffExpr(0) + RHS = RHS_T + RHS_R + for j in 1:J + mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j]) + @constraint(model, sum(x[mᵍ,n]*Sʲs[j][m] + for (m,mᵍ) in enumerate(mʲ_global_idxs)) == RHS) + end + end + for j in 1:J + mʲ_global_idxs = sum(Mʲs[1:j-1])+1:sum(Mʲs[1:j]) + @constraint(model, sum(x[mᵍ,n] for mᵍ in mʲ_global_idxs for n in 1:Nᵀ) == Mʲs[j]) + for mᵍ in mʲ_global_idxs + @constraint(model, sum(x[mᵍ,n] for n in 1:Nᵀ) == 1) + end + end + + # constraints only on `z` + if L ≠ 0 + @constraint(model, sum(z[m,n] for m in 1:N for n in 1:Nᵀ) == N) + for m in 1:N + @constraint(model, sum(z[m,n] for n in 1:Nᵀ) == 1) + end + end + + # constraints only on `y` + if L ≠ 0 + @constraint(model, sum(y[l,n] for l in 1:L for n in 1:Nᵀ) == L) + for l in 1:L + @constraint(model, sum(y[l,n] for n in 1:Nᵀ) == 1) + end + end + + optimize!(model) + + feasible = is_solved_and_feasible(model) + if !(feasible || termination_status(model) == INFEASIBLE) + error(lazy"unexpected termination status: $(termination_status(m))") + end + + # TODO: needs to be modified to allow `Nᵀ` + verbose && _print_subset_sum_result(model, Sʲs, T, R) + + return feasible, model +end \ No newline at end of file diff --git a/BandGraphs/src/types.jl b/BandGraphs/src/types.jl index 554fcddd..79166cdd 100644 --- a/BandGraphs/src/types.jl +++ b/BandGraphs/src/types.jl @@ -18,7 +18,7 @@ irrep labels in `irlabs_nv`, mapping these labels to the full irrep information function SymVector( nv::AbstractVector{<:Integer}, irlabs_nv::AbstractVector{String}, - lgirsd::Dict{String, Vector{LGIrrep{D}}}) where D + lgirsd::Dict{String, <:AbstractVector{LGIrrep{D}}}) where D klabs = unique(klabel.(irlabs_nv)) Nk = length(klabs) mults = [Int[] for _ in 1:Nk] @@ -71,6 +71,19 @@ end Base.collect(n :: SymVector) = reduce(vcat, n.mults) # flatten `n` to a simple vector `nv` +@noinline function _throw_dissimilar_irreps(n₁, n₂) + error("cannot add symmetry vectors with different irreps:\n $(n₁.lgirsv) ≠ $(n₂.lgirsv)") +end +function Base.:+(n₁ :: SymVector{D}, n₂ :: SymVector{D}) where D + length(n₁.lgirsv) == length(n₂.lgirsv) || _throw_dissimilar_irreps(n₁, n₂) + for (lgirs₁, lgirs₂) in zip(n₁.lgirsv, n₂.lgirsv) + length(lgirs₁) == length(lgirs₂) || _throw_dissimilar_irreps(n₁, n₂) + for (lgir₁, lgir₂) in zip(lgirs₁, lgirs₂) + label(lgir₁) == label(lgir₂) || _throw_dissimilar_irreps(n₁, n₂) + end + end + return SymVector(n₁.mults + n₂.mults, n₁.lgirsv, n₁.klabs, n₁.μ + n₂.μ) +end # ---------------------------------------------------------------------------------------- # struct Partition{D} @@ -110,26 +123,14 @@ end function Base.show(io :: IO, s :: SubGraph) print(io, "[") for i in 1:size(s.A, 1) - print(io, label(s.p_max.lgirs[i]), " ↓ ") + print_identified_lgir_label(io, s.p_max.lgirs[i], s.p_max, s.pinned, i) + print(io, " ↓ ") js = findall(!iszero, @view s.A[i,:]) for (nⱼ,j) in enumerate(js) lgirⱼ = s.p_nonmax.lgirs[j] multiplicity = s.A[i,j] // Crystalline.irdim(lgirⱼ) isone(multiplicity) || print(io, multiplicity) - print(io, label(lgirⱼ)) - if !s.pinned - # print a superscript index to identify subduction to nonmaximal irreps for - # which there may be multiple copies; only print an identifier if there - # are multiple "copies", and only if the subgraph is not pinned; the index - # is a local counter within the set of identical/"copied" nonmaximal irreps - for m in s.p_nonmax.multiples - length(m) == 1 && continue # don't print a local index if no copies - local_iridx = findfirst(==(j), m) - isnothing(local_iridx) && continue - printstyled(io, Crystalline.supscriptify(string(local_iridx)); - color=:light_black) - end - end + print_identified_lgir_label(io, lgirⱼ, s.p_nonmax, s.pinned, j) nⱼ == length(js) || print(io, " + ") end i == size(s.A, 1) || print(io, ", ") @@ -137,6 +138,50 @@ function Base.show(io :: IO, s :: SubGraph) print(io, "]") s.pinned && printstyled(io, " (pinned)"; color=:light_black) end + +function print_identified_lgir_label( + io :: IO, lgir :: LGIrrep{D}, p :: Partition{D}, pinned :: Bool, idx :: Int + ) where D + # print a superscript index to identify subduction to nonmaximal irreps for which there + # may be multiple copies; only print an identifier if there are multiple "copies", and + # only if the subgraph is not `pinned`; the index is a local counter within the set of + # identical/"copied" nonmaximal irreps + print(io, label(lgir)) + pinned && return + for m in p.multiples + length(m) == 1 && continue # don't print a local index if no copies + local_iridx = findfirst(==(idx), m) + isnothing(local_iridx) && continue + printstyled(io, Crystalline.supscriptify(string(local_iridx)); color=:light_black) + end +end +function Base.show(io :: IO, ::MIME"text/plain", s :: SubGraph) + summary(io, s) + println(io, ": ", s) + + # max- and nonmax- labels + max_labels = map(1:length(s.p_max)) do i + io′ = IOBuffer() + print_identified_lgir_label(io′, s.p_max.lgirs[i], s.p_max, s.pinned, i) + String(take!(io′)) + end + nonmax_labels = map(1:length(s.p_nonmax)) do i + io′ = IOBuffer() + print_identified_lgir_label(io′, s.p_nonmax.lgirs[i], s.p_nonmax, s.pinned, i) + String(take!(io′)) + end + pretty_table(io, + s.A; + row_labels = max_labels, + header = nonmax_labels, + vlines = [1,], + hlines = [:begin, 1, :end], + row_label_alignment = :l, + alignment = :c, + formatters = (v,i,j) -> iszero(v) ? "·" : string(v), + highlighters = (Highlighter((data,i,j) -> iszero(data[i,j]), crayon"dark_gray"), ) + ) +end # ---------------------------------------------------------------------------------------- # struct BandGraph{D} @@ -150,6 +195,17 @@ function BandGraph( BandGraph{D}(Vector(subgraphs), Vector(partitions)) end +function occupation(bandg::BandGraph) + partitions = bandg.partitions + μ = sum(irdim, first(partitions).lgirs) + for p in @view partitions[2:end] + if sum(irdim, p.lgirs) ≠ μ + error("uneqal number of bands in different partitions!") + end + end + return μ +end + function Base.show(io :: IO, bandg ::BandGraph) print(io, "{") Ns = length(bandg.subgraphs) @@ -183,6 +239,11 @@ mutable struct SubGraphPermutations{D} <: AbstractSubGraph{D} end Base.length(subgraph_ps::SubGraphPermutations) = length(subgraph_ps.As) +function Base.getindex(subgraph_ps::SubGraphPermutations, idx::Integer) + 1 ≤ idx ≤ length(subgraph_ps) || throw(BoundsError(subgraph_ps, idx)) + A = subgraph_ps.As[idx] + return SubGraph(subgraph_ps.p_max, subgraph_ps.p_nonmax, A, subgraph_ps.pinned) +end # ---------------------------------------------------------------------------------------- # struct BandGraphPermutations{D} <: AbstractVector{BandGraph{D}} @@ -192,7 +253,7 @@ end function Base.length(bandgp :: BandGraphPermutations) subgraphs_ps = bandgp.subgraphs_ps - n = 1 + n = BigInt(1) # the number of permutations can be very large; need BigInt to be safe for subgraph_ps in subgraphs_ps n *= length(subgraph_ps) end @@ -207,13 +268,10 @@ end function permutation_info(bandgp :: BandGraphPermutations) subgraphs_ps = bandgp.subgraphs_ps - n = n_no_pins = 1 + n = n_no_pins = BigInt(1) foreach(subgraphs_ps) do subgraph_ps - multiples = subgraph_ps.p_nonmax.multiples - pss = collect.(permutations.(multiples)) - Np = prod(length, pss) # aggregate number of same-irrep permutations, across irreps + Np = length(subgraph_ps) # aggregate number of distinct same-irrep permutations, across irreps if !subgraph_ps.pinned - @assert Np == length(subgraph_ps) n *= Np end n_no_pins *= Np @@ -244,7 +302,7 @@ function Base.getindex( # This is equivalent to the following linear-index-to-subscript-index problem: # CartesianIndices(Tuple(counts))[idx] # or, equivalently, to `collect(Base.product(counts...))[idx]`; both these variants are - # problematic, however, since they are either type-unstable or allocates a lot. + # problematic, however, since they are either type-unstable or allocate a lot. # Instead, we copy the logic from `Base.getindex` on `CartesianIndices`, specifically # from `_ind2sub_recurse`, converting also from a recursive implementation to a loop. @boundscheck 1 ≤ idx ≤ length(bandgp) || throw(BoundsError(bandgp, idx)) diff --git a/BandGraphs/src/unfold.jl b/BandGraphs/src/unfold.jl index 0b04b0ac..2db6def8 100644 --- a/BandGraphs/src/unfold.jl +++ b/BandGraphs/src/unfold.jl @@ -8,7 +8,7 @@ function partition_ids(p :: Partition) end unfold_bandgraph(bandg::BandGraph, kg) = unfold_bandgraph(bandg.subgraphs, bandg.partitions, kg) -function unfold_bandgraph( +@eval BandGraphs function unfold_bandgraph( subgraphs, partitions, kg #= info about connections between k-partitions, from `partition_graph` =# @@ -25,6 +25,13 @@ function unfold_bandgraph( # --- add vertices --- for (trailidx, kidx′) in enumerate(kg_trail.trail) + if kidx′ == -1 + # we've must have more than one disconnected component in the graph; so there + # is not a single disconnected trail we can follow between the vertices; that is + # okay though, there is just no vertex to add to `g_trail` here - and we need to + # keep the absence of a point here in mind when we get to plotting later + continue + end klab = label_for(kg′, kidx′)[1] p = partitionsd[klab] lgir_ids = partition_ids(p) @@ -40,11 +47,20 @@ function unfold_bandgraph( # --- add edges --- subgraphsd = Dict((s.p_max.klab, s.p_nonmax.klab) => s for s in subgraphs) foreach(klab->k_ids[klab]=0, keys(k_ids)) # reset counters + trailidx_begin = firstindex(kg_trail.trail) # peel one iteration - kidx′ = first(kg_trail.trail) + @label CONNECTED_COMPONENT_START + kidx′ = kg_trail.trail[trailidx_begin] p_prev = partitionsd[label_for(kg′, kidx′)[1]] k_id = (k_ids[p_prev.klab] += 1) - for kidx′ in @view kg_trail.trail[2:end] + for trailidx in trailidx_begin+1:lastindex(kg_trail.trail) + kidx′ = kg_trail.trail[trailidx] + if kidx′ == -1 + # nothing to add for trail transitions between disconnected components; restart + # the loop at the next trailidx + trailidx_begin = trailidx + 1 + @goto CONNECTED_COMPONENT_START + end vertex_ids_prev = [(lgir_id..., k_id) for lgir_id in partition_ids(p_prev)] p_curr = partitionsd[label_for(kg′, kidx′)[1]] @@ -68,53 +84,6 @@ function unfold_bandgraph( vertex_ids_curr = vertex_ids_prev end - klabs_trail = [kg′[label_for(kg′, code)].klab for code in kg_trail.trail] + klabs_trail = [code == -1 ? "" : kg′[label_for(kg′, code)].klab for code in kg_trail.trail] return g_trail, klabs_trail -end - - -# WIP below: aim was to remove the nonmaximal vertices and replace them with edge info; but -# stalled due to problem in "NB" below. -# tricky bits include: -# - correctly deal with two max-irreps of high-degen being connected by -# multiple low-degen nonmax-irreps (e.g., `sgnum = 147; brs[end]`); -# - if a max-manifold is not connected via a nonmax manifold, should it be pruned -# (e.g., `sgnum = 6`) or connected via Ω? The latter seems preferable for -# connectivity analysis -function edgify_nonmax_vertices(g_trail) - # FIXME: This does not work in general since the general case would require us to have a - # multigraph structure. To see this, consider e.g. the case that we have an irrep - # Γ₁ which splits into Δ₁+Δ₂ and then recombines to X₁. When we remove the edges - # due to Λ₁ and Λ₂ and seek to make them direct edges between Γ₁ and X₁, we end - # up needing _two_ edges between the Γ₁ and X₁ vertices - - g_decimated = MetaGraph(DiGraph(); - label_type = Tuple{String, Int, Int}, # (irrep label, multiplicity identifier, path-identifier) - vertex_data_type = @NamedTuple{lgir::LGIrrep, maximal::Bool, trailidx::Int}, - edge_data_type = @NamedTuple{lgir::LGIrrep, weight::Int} - ) - original_idxs = Int[] - - for v in vertices(g_trail) - lᵥ = label_for(g_trail, v) - vertex_data = g_trail[lᵥ] - vertex_data.maximal || continue - add_vertex!(g_decimated, lᵥ, vertex_data) - push!(original_idxs, v) - end - - for v in original_idxs - lᵥ = label_for(g_trail, v) - ns = outneighbors(g_trail, v) - for n in ns # over nonmaximal nodes - lₙ = label_for(g_trail, n) - m = only(outneighbors(g_trail, n)) # nonmax nodes can only have one exiting edge - lₘ = label_for(g_trail, m) - edge_data = (; lgir=g_trail[lₙ].lgir, g_trail[lᵥ,lₙ]...) - add_edge!(g_decimated, lᵥ, lₘ, edge_data) || error("edge insertion failure") - end - end - - # original_idxs is useful to e.g., index into the `xy` struct - return g_decimated, original_idxs end \ No newline at end of file diff --git a/BandGraphs/test/Project.toml b/BandGraphs/test/Project.toml index 69328674..c730b2d7 100644 --- a/BandGraphs/test/Project.toml +++ b/BandGraphs/test/Project.toml @@ -4,9 +4,12 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Crystalline = "ae5e2be0-a263-11e9-351e-f94dad1eb351" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" +Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KdotP = "c8045c5a-940e-4537-a0b9-cab4a5eb5ca1" +LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e" MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377" SymmetryBases = "cb68ded9-7e8c-45d3-bc66-19b5867c0467" diff --git a/BandGraphs/test/bandgraph_properties.jl b/BandGraphs/test/bandgraph_properties.jl new file mode 100644 index 00000000..9cb08032 --- /dev/null +++ b/BandGraphs/test/bandgraph_properties.jl @@ -0,0 +1,59 @@ +using Pkg +(dirname(Pkg.project().path) == @__DIR__) || Pkg.activate(@__DIR__) + +using Test +using BandGraphs +using Crystalline +using SymmetryBases + +using Crystalline: irdim + +@testset "Consistency of band graph adjacency matrix blocks" begin +@testset for D in (2,3) + Dᵛ = Val(D) + @testset for sgnum in 1:MAX_SGNUM[D] + # skip 147 and 123 since it takes a long time to do the compatibility basis there + sgnum ∈ (47, 123) && continue + checked_multiplicity = false + for timereversal in (false, true) + subts = subduction_tables(sgnum, Dᵛ; timereversal) + sb, _ = compatibility_basis(sgnum, D; timereversal) + lgirsd = lgirreps(sgnum, Dᵛ) + timereversal && realify!(lgirsd) + for _n in sb + n = SymVector(_n, sb.irlabs, lgirsd) + bandg = build_subgraphs(n, subts, lgirsd) + partitions, subgraphs = bandg.partitions, bandg.subgraphs + + # band occupation & grand-sums of all subgraphs must be constant & equal + occupations = [sum(s.A) for s in subgraphs] + @test all(==(n.μ), occupations) + A = assemble_adjacency(bandg) + @test all(b -> iszero(b) || sum(b) == n.μ, A.blocks) + + # column-wise and row-wise sums of adjacency blocks must give corresponding + # irrep dimensions of columns and rows, respectively + nonmax_irdims = [irdim.(s.p_nonmax.lgirs) for s in subgraphs] + max_irdims = [irdim.(s.p_max.lgirs) for s in subgraphs] + colwise_sums = [vec(sum(s.A; dims=1)) for s in subgraphs] # sums per column (across rows) + rowwise_sums = [vec(sum(s.A; dims=2)) for s in subgraphs] # sums per row (across columns) + + @test colwise_sums == nonmax_irdims + @test rowwise_sums == max_irdims + + # every non-maximal irrep must appear at least twice + if !checked_multiplicity # no need to check every symmetry vector; one is enough + nonmax_multiplicity = [count(s->s.p_nonmax == p, subgraphs) + for p in partitions if !p.maximal] + @test all(≥(2), nonmax_multiplicity) + if !isempty(subgraphs) + @test !isempty(nonmax_multiplicity) + end + check_multiplicity = true + end + end + end + end +end +end # @testset +nothing \ No newline at end of file diff --git a/BandGraphs/test/preliminary_planegroup_work.jl b/BandGraphs/test/preliminary_planegroup_work.jl index 077c42a0..beee833c 100644 --- a/BandGraphs/test/preliminary_planegroup_work.jl +++ b/BandGraphs/test/preliminary_planegroup_work.jl @@ -16,23 +16,27 @@ using GLMakie # Example for p2mm (plane group ⋕6) timereversal = true -sgnum²ᴰ = 17 +sgnum²ᴰ = 12 sgnum³ᴰ = Crystalline.PLANE2SPACE_NUMS[sgnum²ᴰ] subts³ᴰ = subduction_tables(sgnum³ᴰ; timereversal) subts²ᴰ = BandGraphs.SubductionTable{2}[] +cts²ᴰ = BandGraphs.Connection{2}[] for (i, t) in enumerate(subts³ᴰ) if sgnum²ᴰ == 2 xy = [3,1] z = 2 - elseif sgnum²ᴰ == 5 + elseif sgnum²ᴰ in (3, 5) xy = [2,1] z = 3 + elseif sgnum²ᴰ == 4 + xy = [2, 3] + z = 1 else xy = [1,2] z = 3 end - if !iszero(t.c.kᴳ.kv.cnst[z]) || !iszero(t.c.kᴴ.kv.cnst[z])# || !iszero(t.c.kᴴ.kv.free[z,:]) + if !iszero(t.c.kᴳ.kv.cnst[z]) || !iszero(t.c.kᴴ.kv.cnst[z]) #|| !iszero(t.c.kᴴ.kv.free[z,:]) continue # the connection involves k₃≠0 else if !iszero(t.c.kᴴ.kv.free[:,z]) @@ -52,14 +56,16 @@ for (i, t) in enumerate(subts³ᴰ) push!(subts²ᴰ, BandGraphs.SubductionTable{2}(sgnum²ᴰ, c, t.irlabsᴳ, t.irlabsᴴ, table, t.monodromy) ) + push!(cts²ᴰ, c) end +cts²ᴰ ## --------------------------------------------------------------------------------------- # # Visualize "interesting" (i.e., more than 1 band) Hilbert vectors sb, brs = compatibility_basis(sgnum²ᴰ, 2; timereversal) lgirsd = lgirreps(sgnum²ᴰ, Val(2)) -timereversal && (lgirsd = Dict(klab => realify(lgirs) for (klab, lgirs) in lgirsd)) +timereversal && realify!(lgirsd) GLMakie.closeall() for (j, _n) in enumerate(sb) diff --git a/BandGraphs/test/runtests.jl b/BandGraphs/test/runtests.jl index c80af472..bda1fc06 100644 --- a/BandGraphs/test/runtests.jl +++ b/BandGraphs/test/runtests.jl @@ -11,13 +11,15 @@ using GLMakie ## ---------------------------------------------------------------------------------------- ## Testing & visualization of band graph +D = 2 timereversal = true sgnum = 17 -sb, brs = compatibility_basis(sgnum; timereversal) -lgirsd = lgirreps(sgnum) -timereversal && (lgirsd = Dict(klab => realify(lgirs) for (klab, lgirs) in lgirsd)) -subts = subduction_tables(sgnum; timereversal) -_n = sb[1] +sb, brs = compatibility_basis(sgnum, D; timereversal) +lgirsd = lgirreps(sgnum, D) +timereversal && realify!(lgirsd) +subts = subduction_tables(sgnum, D; timereversal) +_n = sb[end] +#_n = brs[end-1] n = SymVector(_n, brs.irlabs, lgirsd) bandg = build_subgraphs(n, subts, lgirsd) @@ -30,7 +32,8 @@ g = assemble_graph(bandg) # structured equiv of `Graph(A)` node_colors = [g[label_for(g, i)].maximal ? :red : :black for i in vertices(g)] f, ax, p = graphplot( g; - nlabels=[(l = label_for(g, i); l[1] * Crystalline.supscriptify(string(l[2]))) for i in vertices(g)], + #nlabels=[(l = label_for(g, i); l[1] * Crystalline.supscriptify(string(l[2]))) for i in vertices(g)], + nlabels=[(l = label_for(g, i); l[1]) for i in vertices(g)], nlabels_distance=4, node_color = node_colors, nlabels_color=node_colors, node_attr = (; strokewidth=4, strokecolor=:white), @@ -65,25 +68,7 @@ graphplot!(ax, dkg, f ## ---------------------------------------------------------------------------------------- - - -sgnum = 19#230 - -sb, brs = compatibility_basis(sgnum; timereversal) -#sb, brs = nontopological_basis(sgnum; timereversal) -lgirsd = lgirreps(sgnum) -timereversal && (lgirsd = Dict(klab => realify(lgirs) for (klab, lgirs) in lgirsd)) -#_n = brs[5] -_n = sb[1] -n = SymVector(_n, brs.irlabs, lgirsd) -subts = subduction_tables(sgnum; timereversal) - -# TODO: Debug segfault(?) of SG 230 (e.g., `sb[end]`) (probably too many permutations for -# some vectors) - -# ---------------------------------------------------------------------------------------- # Create a permuted graph (a single subgraph permutation) and compare -bandg = build_subgraphs(n, subts, lgirsd); subgraphs_ps = permute_subgraphs(bandg.subgraphs); bandgp = BandGraphs.BandGraphPermutations(bandg.partitions, subgraphs_ps); BandGraphs.permutation_info(bandgp) @@ -92,7 +77,7 @@ length(bandgp) > 10 && @info("Be warned: many permutations ($(length(bandgp)))") GLMakie.closeall() xys = nothing maxplot = 4 -fs = Vector{Figure}(undef, maxplot) +fs = Vector{Figure}(undef, min(length(bandgp), maxplot)) for (n, bandg′) in enumerate(bandgp) n > maxplot && break faxp, (; xs, ys) = plot_flattened_bandgraph(bandg′; xys=xys) @@ -106,11 +91,11 @@ end ## ---------------------------------------------------------------------------------------- # Understanding how many permutations remain -let sgnum = 96 # 200 +let sgnum = 96, D = 3 # 200 sb, brs = compatibility_basis(sgnum; timereversal) #sb, brs = nontopological_basis(sgnum; timereversal) lgirsd = Dict(klab=>realify(lgirs) for (klab, lgirs) in lgirreps(sgnum)) - subts = subduction_tables(sgnum; timereversal) + subts = subduction_tables(sgnum, D; timereversal) for nv in sb n = SymVector(nv, brs.irlabs, lgirsd); diff --git a/BandGraphs/test/seperable-irreps-testbed.jl b/BandGraphs/test/seperable-irreps-testbed.jl new file mode 100644 index 00000000..89e1cf64 --- /dev/null +++ b/BandGraphs/test/seperable-irreps-testbed.jl @@ -0,0 +1,136 @@ +using Pkg +(dirname(Pkg.project().path) == @__DIR__) || Pkg.activate(@__DIR__) +using Crystalline +using BandGraphs +using Graphs +using SymmetryBases +using GraphMakie +using GLMakie +using BandGraphs: findall_separable_vertices +using Crystalline: irdim +using KdotP +using Dictionaries + +## --------------------------------------------------------------------------------------- # +#logpath = joinpath((@__DIR__), "log.txt") +#redirect_stdio(; stdout=logpath, stderr=logpath, stdin=devnull) do +## --------------------------------------------------------------------------------------- # + +D = 3 +timereversal = true +criterion = lgir -> isspecial(lgir) && irdim(lgir) == 3 # KdotP.isweyl(lgir; timereversal) +separable_degree = 2 + +sgnums = 1:MAX_SGNUM[D] +contendersd = Dictionary{Int, Vector{SymVector{D}}}() +lgirsdd = Dictionary{Int, Dict{String, Crystalline.IrrepCollection{LGIrrep{D}}}}() +subtsd = Dictionary{Int, Vector{SubductionTable{D}}}() +for sgnum in sgnums + lgirsd = lgirreps(sgnum, D) + timereversal && realify!(lgirsd) + has_relevant_lgirs = any(lgirsd) do (klab, lgirs) + any(criterion, lgirs) + end + if has_relevant_lgirs + insert!(contendersd, sgnum, SymVector{D}[]) + insert!(lgirsdd, sgnum, lgirsd) + insert!(subtsd, sgnum, subduction_tables(sgnum, D; timereversal)) + else + continue + end + printstyled("#", sgnum; bold=true) + + sb, brs = compatibility_basis(sgnum, D; timereversal) + μs = fillings(sb) + for i in eachindex(sb) + μs[i] == 1 && continue # nothing interesting in 1-band cases + + n = SymVector(sb[i], brs.irlabs, lgirsd) + for (lgirs, mults) in zip(n.lgirsv, n.mults) + for (lgir, m) in zip(lgirs, mults) + if m ≠ 0 && criterion(lgir) + push!(contendersd[sgnum], n) + @goto SYMVECTOR_DONE + end + end + end + @label SYMVECTOR_DONE + end + printstyled(" (", length(contendersd[sgnum]), " contenders)\n"; color=:light_black) + flush(stdout) + flush(stderr) +end + + +## --------------------------------------------------------------------------------------- # +separable_summary = Dictionary{Int, Vector{Tuple{Vector{LGIrrep{D}}, SymVector{D}}}}() +separable_details = Dictionary{Int, Vector{Tuple{SymVector{D}, Vector{Tuple{BandGraph{D}, BandGraph{D}, LGIrrep{D}}}}}}() +inseparable = Dictionary{Int, Vector{SymVector{D}}}() + +done = Set{Int}() +too_hard = Set{Int}() +for (sgnum, ns) in pairs(contendersd) + #sgnum == 84 && continue + #sgnum ∈ (178, 181, 179, 180, 214, 199) && continue # untractably many permutations + #sgnum < 211 && continue + #sgnum == 210 && continue # segfault!? + #sgnum ≠ 214 && continue # know there has to be examples here (gyroid)! segfault!? + #sgnum ≠ 230 && continue # know there has to be examples here (double gyroid)! segfault!? + printstyled("#", sgnum, " "; bold=true) + printstyled("(", length(ns), " contenders)\n", color=:light_black) + for (i, n) in enumerate(ns) + #i ≤ 324 && continue + println(" i=", i, ": ", n) + subts = subtsd[sgnum] + lgirsd = lgirsdd[sgnum] + + try + tmp = findall_separable_vertices( + criterion, n, subts, lgirsd; + separable_degree = separable_degree, + max_permutations = 1e5) + if ismissing(tmp) + push!(too_hard, sgnum) + #printstyled(" Skipping any remaining symmetry vectors in this space group!\n"; + # color=:red) + continue + #break # too many permutations + end + has_split, bandg_splits = tmp + if has_split + lgir = getindex.(bandg_splits, 3) + push!(get!(separable_summary, sgnum, Vector{SymVector{D}}()), (lgir, n)) + push!(get!(separable_details, sgnum, + Vector{Tuple{SymVector{D}, typeof(bandg_splits)}}()), + (n, bandg_splits)) + printstyled(" HAS SEPARABLE CONFIGURATIONS!!!\n", color=:green) + else + push!(get!(inseparable, sgnum, Vector{SymVector{D}}()), n) + end + catch e + printstyled(" ", e, "\n", + #" Skipping any remaining symmetry vectors in this space group!\n"; + color=:red) + push!(too_hard, sgnum) + continue + #break # too many permutations + end + end + if sgnum ∉ too_hard + push!(done, sgnum) + end + flush(stdout) + flush(stderr) +end +separable_summary +#end +## --------------------------------------------------------------------------------------- # + +using KdotP +degeneracies = Dictionary{Int, Vector{KdotP.HamiltonianExpansion{D}}}() +for (sgnum, info) in pairs(separable_summary) + lgirs = unique(reduce(vcat, getindex.(info, 1))) + inset!(degeneracies, sgnum, kdotp.(lgirs)) +end + +## --------------------------------------------------------------------------------------- # diff --git a/BandGraphs/test/splittable-2d-irreps-in-planegroups.jl b/BandGraphs/test/splittable-2d-irreps-in-planegroups.jl new file mode 100644 index 00000000..ab825eac --- /dev/null +++ b/BandGraphs/test/splittable-2d-irreps-in-planegroups.jl @@ -0,0 +1,83 @@ +using Pkg +(dirname(Pkg.project().path) == @__DIR__) || Pkg.activate(@__DIR__) +using Crystalline +using BandGraphs +using Graphs +using SymmetryBases +using GraphMakie +using MetaGraphsNext +using GLMakie +using BandGraphs: BandGraphPermutations +using BandGraphs: findall_separable_vertices +using Crystalline: irdim +## --------------------------------------------------------------------------------------- # + +D = 2 +timereversal = false +criterion = (lgir) -> irdim(lgir) == 2 + +contenders = Vector{Pair{Int, SymVector{D}}}() +lgirsdd = Dict{Int, Dict{String, Crystalline.IrrepCollection{LGIrrep{D}}}}() +subtsd = Dict{Int, Vector{SubductionTable{D}}}() +for sgnum in 1:MAX_SGNUM[D] + #sgnum ∈ (4) && continue # TODO: generalize to handle disconnected partitions, but #4 is not separable anyway + lgirsd = lgirreps(sgnum, D) + timereversal && realify!(lgirsd) + + sb, brs = compatibility_basis(sgnum, D; timereversal) + μs = fillings(sb) + had_contender = false + for i in eachindex(sb) + μs[i] == 1 && continue # nothing interesting in 1-band cases + + _n = sb[i] + n = SymVector(_n, brs.irlabs, lgirsd) + for (lgirs, mults) in zip(n.lgirsv, n.mults) + for (lgir, m) in zip(lgirs, mults) + if m ≠ 0 && criterion(lgir) + push!(contenders, sgnum => n) + had_contender = true + @goto SYMVECTOR_DONE + end + end + end + @label SYMVECTOR_DONE + end + if had_contender + lgirsdd[sgnum] = lgirsd + subtsd[sgnum] = subduction_tables(sgnum, D; timereversal) + end +end + + +## --------------------------------------------------------------------------------------- # +separable_summary = Dict{Int, Vector{Tuple{Vector{LGIrrep{D}}, SymVector{D}}}}() +separable_details = Dict{Int, Vector{Tuple{SymVector{D}, Vector{Tuple{BandGraph{D}, BandGraph{D}, LGIrrep{D}}}}}}() +inseparable = Dict{Int, Vector{SymVector{D}}}() + +for (i, (sgnum, n)) in enumerate(contenders) + subts = subtsd[sgnum] + lgirsd = lgirsdd[sgnum] + has_split, bandg_splits = findall_separable_vertices(criterion, n, subts, lgirsd) + if has_split + lgir = getindex.(bandg_splits, 3) + push!(get!(separable_summary, sgnum, Vector{SymVector{D}}()), (lgir, n)) + push!(get!(separable_details, sgnum, + Vector{Tuple{SymVector{D}, typeof(bandg_splits)}}()), + (n, bandg_splits)) + else + push!(get!(inseparable, sgnum, Vector{SymVector{D}}()), n) + end +end +separable_summary + +## --------------------------------------------------------------------------------------- # + +using KdotP +degeneracies = Dict{Int, Vector{KdotP.HamiltonianExpansion{D}}}() +for (sgnum, info) in separable_summary + lgirs = unique(reduce(vcat, getindex.(info, 1))) + degeneracies[sgnum] = kdotp.(lgirs) +end + +## --------------------------------------------------------------------------------------- # diff --git a/BandGraphs/test/subduction-table.jl b/BandGraphs/test/subduction-table.jl index eda85543..d0f543bd 100644 --- a/BandGraphs/test/subduction-table.jl +++ b/BandGraphs/test/subduction-table.jl @@ -6,23 +6,22 @@ using BandGraphs using JLD2 using Test -connectionsd = load_object(joinpath((@__DIR__), "..", "data/connections/3d/connections.jld2")) -subductionsd = load_object(joinpath((@__DIR__), "..", "data/connections/3d/subductions.jld2")) +timereversal = true +trtag = timereversal ? "-tr" : "" +dir = joinpath((@__DIR__), "..", "data/connections/3d") +connectionsd = load_object(joinpath(dir, "connections$(trtag).jld2")) +subductionsd = load_object(joinpath(dir, "subductions$(trtag).jld2")) subductionsd′ = Dict{Int, Vector{SubductionTable{3}}}() failures = Dict{Int, Any}() errors = Dict{Int, Any}() for sgnum in 1:230 subts = subductionsd[sgnum] - # ERROR: nonzero reciprocal vector between nonspecial kv and kv′ requires explicit handling: not yet implemented - sgnum ∈ (26, 27, 28, 29, 30, 31, 32, 33, 34, 39, 40, 41, 100, 101, 102, 103, 104, 106, 116, 117, 118, 226) && continue - # ERROR: DomainError with Provided irreps are not H isspecial(lgir) && irdim(lgir) == 2 + +irs_d = Dict{Int, Vector{LGIrrep{3}}}() +#checked_vectors_d = Dict{Int, Int}() +#n_permutations_d = Dict{Int, BigInt}() +for sgnum in 1:230 + sgnum == 1 && empty!(irs_d) + # --- identify irreps that satisfy `criterion` --- + lgirsd = lgirreps(sgnum) + timereversal && realify!(lgirsd) + for lgirs in values(lgirsd) + isspecial(first(lgirs)) || continue + for lgir in lgirs + if criterion(lgir) + haskey(irs_d, sgnum) || (irs_d[sgnum] = Vector{String}()) + push!(irs_d[sgnum], lgir) + end + end + end + haskey(irs_d, sgnum) || continue # skip further work if `criterion` is not satisfied + printstyled("#", sgnum, "\n"; bold=true) + irlabs = label.(irs_d[sgnum]) + + # --- find Hilbert basis vectors w/ discovered irreps --- + #= + subts = subduction_tables(sgnum; timereversal) + sb, _ = compatibility_basis(sgnum; timereversal) + checked_vectors = 0 + n_permutations = BigInt(0) + for _n in sb + n = SymVector(_n, sb.irlabs, lgirsd) + any(irlab -> contains(string(n), irlab), irlabs) || continue + checked_vectors += 1 + + bandg = build_subgraphs(n, subts, lgirsd) + subgraphs_ps = permute_subgraphs(bandg.subgraphs) + bandgp = BandGraphs.BandGraphPermutations(bandg.partitions, subgraphs_ps) + n_permutations += prod(BigInt.(BandGraphs.permutation_counts(bandgp))) + end + checked_vectors_d[sgnum] = checked_vectors + n_permutations_d[sgnum] = n_permutations + =# +end + +## --------------------------------------------------------------------------------------- # +# --- subgroup check --- +# check whether the `criterion`-valid irreps of supergroup G subduces into one of the +# `criterion`-valid irreps of a subgroup H; if not, the group-relation is not useful +irs_gr = MetaGraph(DiGraph(), Int, Vector{LGIrrep{3}}, Vector{Vector{LGIrrep{3}}}) +for (sgnum, lgirs) in irs_d + add_vertex!(irs_gr, sgnum, lgirs) +end + +for sgnumᴳ in 1:230 + haskey(irs_d, sgnumᴳ) || continue + + # --- basic check of subgroup compatibility --- + gr = maximal_subgroups(sgnumᴳ) + previously_seen_idxs = findall(gr[sgnumᴳ]) do rᴴ + sgnumᴴ = rᴴ.num + sgnumᴴ ∈ keys(irs_d) + end + previously_seen = getfield.(gr[sgnumᴳ][previously_seen_idxs], Ref(:num)) + + # --- print info --- + printstyled("G: ", sgnumᴳ, " (", join(label.(irs_d[sgnumᴳ]), ", "), ")\n"; + bold=true, color=:blue) + #println(" vecs (permutations) : ", + # checked_vectors_d[sgnumᴳ]) + # " (", n_permutations_d[sgnumᴳ], ")") + if !isempty(previously_seen) + print(" supergroup of : ") + join(stdout, previously_seen, ", ") + println() + else + printstyled(" base-case\n"; color=:red) + end + + # --- more careful subgroup/irrep/k-vector check --- + irsᴳ = irs_d[sgnumᴳ] + sgᴳ = reduce_ops(spacegroup(sgnumᴳ, Val(3)), centering(sgnumᴳ, 3)) + for (i, sgnumᴴ) in enumerate(previously_seen) + printstyled(" Gisapprox(kᴳᴴ, kᴴ, nothing, false), kᴳᴴs) + if isnothing(idx_kᴳᴴs_eq_kᴴ) + # kᴳ is not equivalent by either symmetry or reciprocal lattice + # vectors to kᴴ; we cannot compare the little groups + printstyled("÷", color=:red) + printstyled(" (k-points)"; color=:light_black) + continue + end + idx_kᴳᴴs_eq_kᴴ = something(idx_kᴳᴴs_eq_kᴴ) + q_coset = sgᴳ[idx_kᴳᴴs_eq_kᴴ] # relevant coset representative + kᴳ′ᴴ = kᴳᴴs[idx_kᴳᴴs_eq_kᴴ] + + # since kᴳ might not equal kᴴ, we might need to remap the operations in + # irrep according to the coset-representative that connects kᴳ to kᴴ; + # but we need to compare the two k-points in a shared setting; we choose + # to first map kᴴ to the setting of kᴳ, calling this kᴴᴳ + kᴴᴳ = transform(kᴴ, inv(P)) # kᴴ in G basis + tmp = Crystalline.remap_lgirreps_to_point_in_kstar([irᴳ], kᴴᴳ, [q_coset]) + irᴳ′ = tmp[1] + lgᴳ′ = group(irᴳ′) + + # now transform to the basis of H (lgᴳ ops in H setting at kᴴ position) + lgᴳ′ᴴ_ops = transform.(lgᴳ′, Ref(P), Ref(p), false) + + # the operator-ordering of the H and G irreps might differ; find the + # permutation between the two (or, if their little groups are not + # subgroups, bail out) + boolsub, idxsᴳ²ᴴ = Crystalline._findsubgroup(lgᴳ′ᴴ_ops, operations(lgᴴ), cntrᴴ) + if !boolsub # we do not have H < G + printstyled("÷", color=:red) + printstyled(" (little groups)"; color=:light_black) + continue + end + + # even if we now have a subgroup, `lgᴳ′ᴴ` and `lgᴴ` may still differ + # in their operations by primitive lattice vectors; and this may change + # the irreps of nonsymmorphic groups; so we have to correct for that + # by identifying possible translation mismatches between operations + # and multiplying by the appropriate Bloch phase + Δts = translation.(lgᴴ) .- translation.(lgᴳ′ᴴ_ops[idxsᴳ²ᴴ]) + translationsᴳ′ = [copy(τ) for τ in irᴳ′.translations] + #matricesᴳ′ = [copy(m) for m in irᴳ′.matrices] + for (i, Δt) in enumerate(Δts) + norm(Δt) > Crystalline.DEFAULT_ATOL || continue + iᴳ = idxsᴳ²ᴴ[i] + lgᴳ′ᴴ_ops[iᴳ] = SymOperation(Δt) * lgᴳ′ᴴ_ops[iᴳ] + translationsᴳ′[iᴳ] += Δt + #matricesᴳ′[iᴳ] .*= cispi(2*dot(kᴳ′ᴴ.cnst, Δt)) + # TODO: do something similar w/ irᴳ.translations & kᴳ′ᴴ.free + end + lgᴳ′ᴴ = LittleGroup{3}(sgnumᴳ, kᴳ′ᴴ, klabel(lgᴳ), lgᴳ′ᴴ_ops[idxsᴳ²ᴴ]) + irᴳ′ᴴ = LGIrrep{3}(label(irᴳ), lgᴳ′ᴴ, irᴳ.matrices[idxsᴳ²ᴴ], + translationsᴳ′[idxsᴳ²ᴴ], irᴳ.reality, irᴳ.iscorep) + + if subduction_count(irᴳ′ᴴ, irᴴ) ≠ 0 + printstyled("✓"; color=:green) + candidate=true + + if !has_edge(irs_gr, code_for(irs_gr, sgnumᴳ), code_for(irs_gr, sgnumᴴ)) + edge_data = [Vector{LGIrrep{3}}() for _ in 1:length(irsᴳ)] + push!(edge_data[j], irᴴ) + add_edge!(irs_gr, sgnumᴳ, sgnumᴴ, edge_data) + else + push!(irs_gr[sgnumᴳ, sgnumᴴ][j], irᴴ) + end + else + printstyled("÷", color=:red) + printstyled(" (subduction)"; color=:light_black) + end + end + println() + end + end + if !candidate + printstyled(" no relevant irrep subductions\n"; color=:red) + end + end +end + +## ----------------- +function layout_y_position(irs_gr) + sgnums = collect(labels(irs_gr)) + orders = Crystalline.SG_PRIMITIVE_ORDERs[3][sgnums] + layers′ = CrystallineGraphMakieExt.map_numbers_to_oneto(orders; rev = false) + maxlayer = maximum(layers′) + 1:length(orders) .=> (maxlayer+1) .- layers′ +end + +const CrystallineGraphMakieExt = Base.get_extension(Crystalline, :CrystallineGraphMakieExt) +const Zarate = CrystallineGraphMakieExt.Zarate +const solve_positions = CrystallineGraphMakieExt.solve_positions +function zarate_layout(irs_gr′, force_y_layer) + # `solve_positions` implicitly assumes `SimpleGraph`/`SimpleDiGraph` type graphs, so we + # must first convert `gr` to `gr′::SimpleDiGraph` + gr′ = SimpleDiGraphFromIterator(edges(irs_gr′)) + xs, ys, _ = solve_positions(Zarate(), gr′; force_layer=force_y_layer) + # TODO: This is broken in that it does not respect the `force_y_layer` if it is not + # consecutive - so cannot be compared across disconnected components + return GraphMakie.Point2f.(ys, -xs) +end + + +_tmp = induced_subgraph.(Ref(irs_gr), connected_components(irs_gr)) +irs_connected_grs = getindex.(_tmp, 1) +irs_connected_idxs = getindex.(_tmp, 2) +y_layers = layout_y_position(irs_gr) +connected_y_layers = [1:length(idxs) .=> last.(y_layers[idxs]) for idxs in irs_connected_idxs] +f = Figure() +N_connected = length(irs_connected_grs) +for (n, irs_gr′) in enumerate(irs_connected_grs) + nlabels_text = map(vertices(irs_gr′)) do i # label for each vertex + sgnum = label_for(irs_gr′, i) + irlabs = join(label.(irs_d[sgnum]), ", ") + string(sgnum) * "\n(" * irlabs * ")" + end + nlabels_color = map(vertices(irs_gr′)) do i + sgnum = label_for(irs_gr′, i) + isempty(outneighbors(irs_gr′, i)) ? :red : :black + end + elabels_text = map(edges(irs_gr′)) do e + sgnumᴳ, sgnumᴴ = label_for(irs_gr′, src(e)), label_for(irs_gr′, dst(e)) + irlabsᴳ = label.(irs_d[sgnumᴳ]) + irsᴴv = irs_gr′[sgnumᴳ, sgnumᴴ] + irlabsᴴ = map(irsᴴv) do irsᴴ + irlabsᴴ = label.(irsᴴ) + length(irlabsᴴ) == 1 ? irlabsᴴ[1] : "(" * join(irlabsᴴ, ", ") * ")" + end + join([irlabsᴳ[i] * " ↓ " * irlabsᴴ[i] for i in 1:length(irlabsᴳ)], "\n") + end + + ax = Axis(f[1, n]) + p = graphplot!(ax, irs_gr′; + nlabels = nlabels_text, + nlabels_align = (:center, :bottom), + nlabels_offset = GraphMakie.Point(0, .025), + nlabels_color = nlabels_color, + nlabels_fontsize = 13, + #elabels = elabels_text, + #elabels_offset = GraphMakie.Point(0, .025), + #elabels_fontsize = 13, + ) + + # adjust appearance of graph + Makie.hidedecorations!(ax) + Makie.hidespines!(ax) + #xy′ = zarate_layout(irs_gr′, connected_y_layers[n]) + #p[:node_pos][] = xy′ # update layout of vertices + #y_extrema = (-1) .* reverse(extrema(last.(y_layers))) + #y_extrema = y_extrema .+ (0.1 * (-)(y_extrema...)) .* (1, -1) + #Makie.ylims!(y_extrema...) +end +f \ No newline at end of file diff --git a/BandGraphs/test/test-split.jl b/BandGraphs/test/test-split.jl new file mode 100644 index 00000000..5f33c225 --- /dev/null +++ b/BandGraphs/test/test-split.jl @@ -0,0 +1,63 @@ +using Pkg +(dirname(Pkg.project().path) == @__DIR__) || Pkg.activate(@__DIR__) +using Crystalline +using BandGraphs +using Graphs +using SymmetryBases +using GraphMakie +using MetaGraphsNext +using GLMakie + +## --------------------------------------------------------------------------------------- # +include(joinpath((@__DIR__), "..", "src", "complete-split.jl")) + +## --------------------------------------------------------------------------------------- # + +D = 2 +timereversal = true +sgnum = 17 +sb, brs = compatibility_basis(sgnum, D; timereversal) +lgirsd = lgirreps(sgnum, D) +timereversal && realify!(lgirsd) +subts = subduction_tables(sgnum, D; timereversal) +_n = sb[end] +n = SymVector(_n, brs.irlabs, lgirsd) + +bandg = build_subgraphs(n, subts, lgirsd) + +## --------------------------------------------------------------------------------------- # + +v = ("K₃", 1) +bandgs′ = something(complete_split(bandg, v)); + +## --------------------------------------------------------------------------------------- # + +g = assemble_graph(bandgs′[1]) + +node_colors = [g[label_for(g, i)].maximal ? :red : :black for i in vertices(g)] +f, ax, p = graphplot( + g; + nlabels=[(l = label_for(g, i); l[1]) for i in vertices(g)], + nlabels_distance=4, + node_color = node_colors, nlabels_color=node_colors, + node_attr = (; strokewidth=4, strokecolor=:white), + edge_color = :gray, + layout = GraphMakie.NetworkLayout.SFDP(; C=10, K=.2, dim=2) + ) + +make_vertices_dragable!(ax, p) +f + +## ---------------------------------------------------------------------------------------- +# Create a permuted graph (a single subgraph permutation) and compare +subgraphs_ps = permute_subgraphs(bandg.subgraphs); +bandgp = BandGraphs.BandGraphPermutations(bandg.partitions, subgraphs_ps); +BandGraphs.permutation_info(bandgp) +length(bandgp) > 10 && @info("Be warned: many permutations ($(length(bandgp)))") + +v = ("K₃", 1) +bandgs′ = something(complete_split(bandgp[1], v)); +bandg′ = bandgs′[2] + +faxp, (; xs, ys) = plot_flattened_bandgraph(bandg′) +display(GLMakie.Screen(), faxp) diff --git a/BandGraphs/test/weyl-notes.ipynb b/BandGraphs/test/weyl-notes.ipynb index dcd43ee6..e2113773 100644 --- a/BandGraphs/test/weyl-notes.ipynb +++ b/BandGraphs/test/weyl-notes.ipynb @@ -85,7 +85,7 @@ "source": [ "function plot_irrep_filtered_bandgraphs(sgnum, selected_irlabs; timereversal=true, maxplot=Inf)\n", " lgirsd = lgirreps(sgnum)\n", - " timereversal && (lgirsd = Dict(klab=>realify(lgirs) for (klab, lgirs) in lgirsd))\n", + " timereversal && realify!(lgirsd)\n", " subts = subduction_tables(sgnum; timereversal)\n", " sb, brs = compatibility_basis(sgnum; timereversal)\n", "\n", @@ -1724,7 +1724,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.10.2", + "display_name": "Julia 1.10.3", "language": "julia", "name": "julia-1.10" }, @@ -1732,7 +1732,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.10.2" + "version": "1.10.3" } }, "nbformat": 4, diff --git a/BandGraphs/test/weyl_bands.jl b/BandGraphs/test/weyl_bands.jl index d72ab129..391c6c00 100644 --- a/BandGraphs/test/weyl_bands.jl +++ b/BandGraphs/test/weyl_bands.jl @@ -20,7 +20,7 @@ for sgnum in 1:230 sgnum == 1 && empty!(weyl_irs_d) # --- identify space groups with Weyl points --- lgirsd = lgirreps(sgnum) - timereversal && (lgirsd = Dict(klab=>realify(lgirs) for (klab, lgirs) in lgirsd)) + timereversal && realify!(lgirsd) for lgirs in values(lgirsd) isspecial(first(lgirs)) || continue for lgir in lgirs @@ -125,43 +125,50 @@ for sgnumᴳ in 1:230 end idx_kᴳᴴs_eq_kᴴ = something(idx_kᴳᴴs_eq_kᴴ) q_coset = sgᴳ[idx_kᴳᴴs_eq_kᴴ] # relevant coset representative - kᴳ′ᴴ = kᴳᴴs[idx_kᴳᴴs_eq_kᴴ] - # rotate `lgᴳ` to the new representative associated with `q_coset` - lgᴳ′ = transform.(lgᴳ, Ref(rotation(q_coset)), Ref(Crystalline.translation(q_coset)), false) + + # since kᴳ might not equal kᴴ, we might need to remap the operations in + # irrep according to the coset-representative that connects kᴳ to kᴴ; + # but we need to compare the two k-points in a shared setting; we choose + # to first map kᴴ to the setting of kᴳ, calling this kᴴᴳ + kᴴᴳ = transform(kᴴ, inv(P)) # kᴴ in G basis + tmp = Crystalline.remap_lgirreps_to_point_in_kstar([irᴳ], kᴴᴳ, [q_coset]) + irᴳ′ = tmp[1] + lgᴳ′ = group(irᴳ′) + # now transform to the basis of H (lgᴳ ops in H setting at kᴴ position) lgᴳ′ᴴ_ops = transform.(lgᴳ′, Ref(P), Ref(p), false) - - # TODO: this is probably not really right since the kᴳᴴ point is only - # guaranteed to be in the orbit of kᴳ, and so the sorting of - # irreps should probably be remapped under the mapping from kᴳ - # to the representative from the orbit + # the operator-ordering of the H and G irreps might differ; find the + # permutation between the two (or, if their little groups are not + # subgroups, bail out) boolsub, idxsᴳ²ᴴ = Crystalline._findsubgroup(lgᴳ′ᴴ_ops, operations(lgᴴ), cntrᴴ) if !boolsub # we do not have H < G printstyled("÷", color=:red) printstyled(" (little groups)"; color=:light_black) + continue end + # even if we now have a subgroup, `lgᴳ′ᴴ` and `lgᴴ` may still differ # in their operations by primitive lattice vectors; and this may change # the irreps of nonsymmorphic groups; so we have to correct for that # by identifying possible translation mismatches between operations # and multiplying by the appropriate Bloch phase Δts = translation.(lgᴴ) .- translation.(lgᴳ′ᴴ_ops[idxsᴳ²ᴴ]) - matricesᴳ′ = [copy(m) for m in irᴳ.matrices] + translationsᴳ′ = [copy(τ) for τ in irᴳ′.translations] + #matricesᴳ′ = [copy(m) for m in irᴳ′.matrices] for (i, Δt) in enumerate(Δts) norm(Δt) > Crystalline.DEFAULT_ATOL || continue iᴳ = idxsᴳ²ᴴ[i] lgᴳ′ᴴ_ops[iᴳ] = SymOperation(Δt) * lgᴳ′ᴴ_ops[iᴳ] - matricesᴳ′[iᴳ] .*= cispi(2*dot(kᴳ′ᴴ.cnst, Δt)) + translationsᴳ′[iᴳ] += Δt + #matricesᴳ′[iᴳ] .*= cispi(2*dot(kᴳ′ᴴ.cnst, Δt)) # TODO: do something similar w/ irᴳ.translations & kᴳ′ᴴ.free end lgᴳ′ᴴ = LittleGroup{3}(sgnumᴳ, kᴳ′ᴴ, klabel(lgᴳ), lgᴳ′ᴴ_ops[idxsᴳ²ᴴ]) - irᴳ′ᴴ = LGIrrep{3}(label(irᴳ), lgᴳ′ᴴ, matricesᴳ′[idxsᴳ²ᴴ], - irᴳ.translations[idxsᴳ²ᴴ], irᴳ.reality, irᴳ.iscorep) - irᴴ = LGIrrep{3}(label(irᴴ), lgᴴ, irᴴ.matrices, - irᴴ.translations, irᴴ.reality, irᴴ.iscorep) + irᴳ′ᴴ = LGIrrep{3}(label(irᴳ), lgᴳ′ᴴ, irᴳ.matrices[idxsᴳ²ᴴ], + translationsᴳ′[idxsᴳ²ᴴ], irᴳ.reality, irᴳ.iscorep) if subduction_count(irᴳ′ᴴ, irᴴ) ≠ 0 printstyled("✓"; color=:green) @@ -186,7 +193,7 @@ for sgnumᴳ in 1:230 printstyled(" no relevant irrep subductions\n"; color=:red) end end -end +end ## ----------------- function layout_y_position(weyl_gr) diff --git a/build/crawl_and_write_bandpaths.jl b/build/crawl_and_write_bandpaths.jl index f1b149d6..047a364c 100644 --- a/build/crawl_and_write_bandpaths.jl +++ b/build/crawl_and_write_bandpaths.jl @@ -43,7 +43,11 @@ function extract_bandpaths_tables(body::HTMLDocument) matches = eachmatch(select, root) bandpath_htmltable, subduction_htmltable = matches[2], matches[3] - return convert_html_table(bandpath_htmltable), convert_html_table(subduction_htmltable) + bandpath_data = convert_html_table(bandpath_htmltable) + subduction_data = convert_html_table(subduction_htmltable) + canonicalize_subduction_data_structure!(subduction_data) + + return bandpath_data, subduction_data end function convert_html_table(t::HTMLElement{:table}) @@ -91,6 +95,39 @@ function html2unicode_postprocess(s::AbstractString) return replace(replace(s, replacepairlist...), " "=>"") end +function canonicalize_subduction_data_structure!(subduction_data::Matrix) + # The structure of the `subduction_data` table, as returned from BCS, is a bit finnicky. + # If there are no monodromy additions, the table has 5 columns, structured as: + # | kᴳ₁ | kᴳ₁-kᴴ rules | | kᴴ | kᴳ₂-kᴴ rules | kᴳ₂ | + # If there _are_ monodromy additions, the table has 7 columns, structured as: + # | kᴳ₁ | kᴳ₁-kᴴ rules | kᴳ₁′-kᴴ rules | kᴴ | (...)² | kᴳ₂ | + # where by kᴳ₁′ and kᴳ₂′ we denote a monodromy-derived k-vector and where (...)² + # represents represents 2 columns of rules, with a variable structure that depends on + # whether there is a monodromy rule or not (whose presence is only optional): + # A. Has monodromy rule: (...)² = | kᴳ₂-kᴴ rules | kᴳ₂′-kᴴ rules | + # B. No monodromy rule: (...)² = | empty string ("") | kᴳ₂-kᴴ rules | + # This is annoying because the structure is variable and depends on the presence of + # monodromy rules. + # To fix this, we simply reorder the columns to ensure the fixed following structure: + # | kᴳ₁ | kᴳ₁-kᴴ rules | kᴳ₁′-kᴴ rules | kᴴ | kᴳ₂′-kᴴ rules | kᴳ₂-kᴴ rules | kᴳ₂ | + Ncols = size(subduction_data, 2) + if Ncols == 7 + for i in 1:size(subduction_data, 1) + tᵢ₅ = subduction_data[i, 5] + if !isempty(first(tᵢ₅)) + # there is a monodromy rule in this row, so we must flip columns 5 & 6 + tᵢ₆ = subduction_data[i, 6] + subduction_data[i, 5] = tᵢ₆ + subduction_data[i, 6] = tᵢ₅ + end + end + end + + # now with `subduction_data` is in a consistent, simple form and we can parse it in + # `parse_subductions` without worrying about a variable location for the different data + return subduction_data +end + # ---------------------------------------------------------------------------------------- # # Types @@ -128,16 +165,9 @@ end # function to find element of little group with highest-order screw (or, failing any screws, # any glide operation); need this to figure out the k-point of a monodromy-related shifted # k-point +function find_highest_order_nonsymmorphic_operation( + g::Crystalline.AbstractGroup{D}, cntr::Char=centering(g, D)) where D -# function find_highest_order_nonsymmorphic_operation(sgnum::Integer) -# issymmorph(sgnum, 3) && return nothing -# sg = spacegroup(sgnum, Val(3)) -# N = round(Int, inv(det(Bravais.primitivebasismatrix(centering(sgnum, 3))))) -# sg = sg[1:div(length(sg), N)] # only include the first "non-centering copies" operations -# -# return find_highest_order_nonsymmorphic_operation(sg, centering(sgnum, 3)) -# end -function find_highest_order_nonsymmorphic_operation(g::Crystalline.AbstractGroup{3}, cntr) ns = Crystalline.rotation_order.(g) ops = iterated_composition.(g, abs.(ns)) for (j,op) in enumerate(ops) @@ -153,7 +183,7 @@ function find_highest_order_nonsymmorphic_operation(g::Crystalline.AbstractGroup t = ts[i] if norm(t) > Crystalline.DEFAULT_ATOL t′ = float.(rationalize.(t, tol=1e-2)) # hack to get rid of floating point errors - if !all(isinteger, primitivize(RVec(t′), centering(g)).cnst) + if !all(isinteger, primitivize(RVec(t′), cntr).cnst) error("obtained non-integer translation $t for operation $(g[i])") end return (g[i], t′) @@ -191,6 +221,7 @@ function parse_subductions( kᴳ_colidxs = (1, Ncols) rules_colidxs = (2, Ncols-1) + # First, we do the "standard" rules (not monodromy-derived) tables = Vector{SubductionTable{3}}() for row in eachrow(subduction_data) for (kᴳ_colidx, rules_colidx) in zip(kᴳ_colidxs, rules_colidxs) @@ -201,7 +232,8 @@ function parse_subductions( end end - if Ncols == 7 # monodromy additions + # Now we do the monodromy additions + if Ncols == 7 lgs = littlegroups(num, 3) for row in eachrow(subduction_data) for (kᴳ_colidx, rules_colidx) in zip(kᴳ_colidxs, (3, Ncols-2)) @@ -212,29 +244,120 @@ function parse_subductions( # plain rules c, irlabsᴳ, irlabsᴴ, table = tmp - # now, we want to infer the k-point associated with the translated high- - # symmetry k-point `kᴳ′` and the associated connection `c′` - kᴳ, kᴴ = c.kᴳ, c.kᴴ - - kᴴlab = kᴴ.label - lg = lgs[string(kᴴlab)] - # FIXME/TODO: the below approach to get the "translated" k-coordinate - # doesn't seem right: need to fix eventually, but not really - # important per se - t = find_highest_order_nonsymmorphic_operation(lg, centering(num, 3))[2] - kᴳ′ = LabeledKVec(Symbol(kᴳ.label, '′'), kᴳ.kv + t) - c′ = Connection(kᴳ′, kᴴ) - irlabsᴳ′ = map(irlab->replace(irlab, string(kᴳ.label) => string(kᴳ′.label)), irlabsᴳ) + # determine the monodromy-related k-point and the associated connection + c′ = find_monodromy_related_connection(c, lgs) + + irlabsᴳ′ = map(irlabsᴳ) do irlab + replace(irlab, string(c.kᴳ.label) => string(c′.kᴳ.label)) + end push!(tables, SubductionTable(num, c′, irlabsᴳ′, irlabsᴴ, table, true)) end end end end - return unique!(t->t.c, tables) end +# infer the k-point associated with the shifted high-symmetry k-point `kᴳ′` and return the +# associated connection `c′` between `kᴳ′` and `kᴴ` +function find_monodromy_related_connection(c, lgs) + kᴳ, kᴴ = c.kᴳ.kv, c.kᴴ.kv + + kᴴlab = c.kᴴ.label + kᴳlab = c.kᴳ.label + kᴳ′lab = Symbol(kᴳlab, '′') + lgᴴ = lgs[string(kᴴlab)] + + cntr = centering(num(lgᴴ), 3) + kᴳₚ = primitivize(kᴳ, cntr) + kᴴₚ = primitivize(kᴴ, cntr) + # FIXME/TODO: the below approach to get the "translated" k-coordinate + # is still not right: must fix + # (Yes, it doesn't make sense to just add a _real-space_ + # translation to the k-vector. Need to get k⋅t = -2πn with + # smallest possible k (and in primitive setting)) + + t = find_highest_order_nonsymmorphic_operation(lgᴴ, cntr)[2] + tₚ = primitivize(RVec(t), cntr).cnst + + free_parts = vec(transpose(tₚ) * kᴴₚ.free) + idxs = findall(v -> norm(v) > Crystalline.DEFAULT_ATOL, free_parts) + idx = if length(idxs) == 1 + last(idxs) + elseif length(idxs) == 2 + # we cannot uniquely determine which k-point to extend to in this situation; + # for now, we just pick _a_ solution, hoping that it will end up being the one that + # will give us the same subduction table (TO BE VERIFIED) + last(idxs) + else + # hoping that this scenario never occurs + _error_info("unexpected situation when trying to solve k⋅t = 1; too many free parameters", kᴳlab, kᴴlab, kᴳ, kᴳₚ, kᴴ, kᴴₚ, "⋅", "⋅", t, tₚ, cntr) + end + kᴳ′ₚ = search_n2π_solution(free_parts, idx, kᴳₚ, kᴴₚ) + if kᴳ′ₚ === nothing + _error_info("failed to find a kᴳ′ that differs from from kᴳ by a primitive reciprocal lattice vector", kᴳlab, kᴴlab, kᴳ, kᴳₚ, kᴴ, kᴴₚ, kᴳ′, kᴳ′ₚ, t, tₚ, cntr) + end + Δkᴳₚ = (kᴳ′ₚ - kᴳₚ).cnst + if !(dot(Δkᴳₚ, tₚ) ≈ round(dot(Δkᴳₚ, tₚ))) + _error_info("Δkᴳ⋅t ∉ ℤ", kᴳlab, kᴴlab, kᴳ, kᴳₚ, kᴴ, kᴴₚ, kᴳ′, kᴳ′ₚ, t, tₚ, cntr) + end + + kᴳ′ = conventionalize(kᴳ′ₚ, cntr) + c′ = Connection(LabeledKVec(kᴳ′lab, kᴳ′), c.kᴴ) + + return c′ +end + +function search_n2π_solution(free_parts, idx′, kᴳₚ, kᴴₚ, n::Integer=1, nmax::Integer=6) + # solve (kᴳ-kᴳ′)⋅t = -n2π for least positive integer `n`, stopping at `nmax` + n>nmax && return nothing + αβγ_idx′ = -n/free_parts[idx′] + kᴳ′ₚ = KVec(kᴳₚ.cnst + kᴴₚ.free[:,idx′] * αβγ_idx′) + # test that kᴳ′ differs from kᴳ by a primitive reciprocal lattice vector + Δkᴳₚ = (kᴳ′ₚ - kᴳₚ).cnst + if all(v->norm(round(v)-v)i≠idx′ && Crystalline.freeparams(kᴴₚ)[i], 1:3) + if !isnothing(idx′′) + v = kᴴₚ.free[:,idx′′] + for c in (1.0, 0.5, -0.5, -1.0) + kᴳ′ₚ = KVec(kᴳ′ₚ.cnst + c*v) + Δkᴳₚ = (kᴳ′ₚ - kᴳₚ).cnst + if all(v->norm(round(v)-v)