From 6f7986d81ad255876ff17135b444b67519a0748c Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 10 Dec 2020 14:22:40 +0100 Subject: [PATCH 1/2] redone bandstructure + splitbands + indexing yet another bandstructure rewrite runs defects works fixes svd for higher-rank projs svd -> svdvals subbands sketch for the future indexing splitbands docstring fixes orthonormalize subspaces better progmeter tests --- src/KPM.jl | 6 +- src/Quantica.jl | 2 +- src/bandstructure.jl | 1070 ++++++++++++++++++------------------ src/diagonalizer.jl | 20 +- src/hamiltonian.jl | 12 +- src/iterators.jl | 71 ++- src/lattice.jl | 8 +- src/mesh.jl | 38 +- src/model.jl | 14 +- src/plot_makie.jl | 11 +- src/plot_vegalite.jl | 6 +- src/tools.jl | 162 ++---- test/test_bandstructure.jl | 18 +- 13 files changed, 733 insertions(+), 705 deletions(-) diff --git a/src/KPM.jl b/src/KPM.jl index 4a8120d3..405ee3d9 100644 --- a/src/KPM.jl +++ b/src/KPM.jl @@ -232,7 +232,7 @@ resolution`, rounded to the closest integer. Same as above with KPM momenta `μ` as input. Equivalent to `densityKPM(μ; kw...)` except that imaginary parts are dropped. -# See also: +# See also `momentaKPM`, `densityKPM`, `averageKPM` """ function dosKPM(h; resolution = 2, kets = randomkets(1), kw...) @@ -263,7 +263,7 @@ for `dosKPM`, all imaginary parts in `ρ_A` are preserved), where the number of Same as above with the KPM momenta as input (see `momentaKPM`). -# See also: +# See also `dosKPM`, `momentaKPM`, `averageKPM` """ densityKPM(h, A; resolution = 2, kw...) = @@ -295,7 +295,7 @@ energy `E_k`, kBT` is the temperature in energy units and `Ef` the Fermi energy. Same as above with the KPM momenta as input (see `momentaKPM`). -# See also: +# See also `dosKPM`, `momentaKPM`, `averageKPM` """ averageKPM(h, A; kBT = 0, Ef = 0, kw...) = averageKPM(momentaKPM(h, A; kw...); kBT = kBT, Ef = Ef) diff --git a/src/Quantica.jl b/src/Quantica.jl index 19b2f3dc..b3822d44 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -25,7 +25,7 @@ export sublat, bravais, lattice, dims, supercell, unitcell, ket, randomkets, hamiltonian, parametric, bloch, bloch!, similarmatrix, flatten, wrap, transform!, combine, - spectrum, bandstructure, diagonalizer, cuboid, isometric, + spectrum, bandstructure, diagonalizer, cuboid, isometric, splitbands, bands, vertices, energies, states, degeneracy, momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM, diff --git a/src/bandstructure.jl b/src/bandstructure.jl index 2bd5034c..655a02fd 100644 --- a/src/bandstructure.jl +++ b/src/bandstructure.jl @@ -1,17 +1,22 @@ ####################################################################### -# Spectrum +# Subspace ####################################################################### -struct Subspace{C,T,S<:AbstractMatrix{C}} +struct Subspace{D,C,T,S<:AbstractMatrix{C}} energy::T basis::S + basevert::SVector{D,T} end -function Base.show(io::IO, s::Subspace{C,T}) where {C,T} +Subspace(energy::T, basis) where {T} = Subspace(energy, basis, SVector{0,T}()) +Subspace(energy::T, basis, basevert) where {T} = Subspace(energy, basis, SVector(T.(basevert))) + +function Base.show(io::IO, s::Subspace{C,T,D}) where {C,T,D} i = get(io, :indent, "") print(io, -"$(i)Subspace{$C,$T}: eigenenergy subspace +"$(i)Subspace{$C,$T,$D}: eigenenergy subspace on a $(D)D manifold $i Energy : $(s.energy) -$i Degeneracy : $(degeneracy(s))") +$i Degeneracy : $(degeneracy(s)) +$i Bloch/params : $(s.basevert)") end """ @@ -22,15 +27,10 @@ Return the degeneracy of a given energy subspace. It is equal to `size(s.basis, # See also `spectrum`, `bandstructure` """ -degeneracy(s::Subspace) = size(s.basis, 2) - -collect_subspaces(ϵs, ψs, ::Type{T}) where {T} = _collect_subspaces(ϵs, ψs, typeof(Subspace(zero(T), view(ψs, :, 1:1)))) -_collect_subspaces(ϵs, ψs, ::Type{SS}) where {C,T,SS<:Subspace{C,T}} = - SS[Subspace(_convert_energy(ϵs, rng, T), view(ψs, :, rng)) for rng in approxruns(ϵs)] - -_convert_energy(ϵs, rng, ::Type{T}) where {T<:Real} = mean(i -> T(real(ϵs[i])), rng) -_convert_energy(ϵs, rng, ::Type{T}) where {T<:Complex} = mean(i -> T(ϵs[i]), rng) +degeneracy(s::Subspace) = degeneracy(s.basis) +degeneracy(m::AbstractMatrix) = isempty(m) ? 1 : size(m, 2) # To support sentinel empty projs +# destructuring Base.iterate(s::Subspace) = s.energy, Val(:basis) Base.iterate(s::Subspace, ::Val{:basis}) = s.basis, Val(:done) Base.iterate(::Subspace, ::Val{:done}) = nothing @@ -39,22 +39,38 @@ Base.first(s::Subspace) = s.energy Base.last(s::Subspace) = s.basis Base.length(s::Subspace) = 2 -struct Spectrum{C,T,E<:AbstractVector{T},A<:AbstractMatrix{C}} +####################################################################### +# Spectrum +####################################################################### +struct Spectrum{D,C,T,S<:AbstractMatrix{C},E<:AbstractVector{T}} energies::E - states::A + states::S + basevert::SVector{D,T} subs::Vector{UnitRange{Int}} - subs´::Vector{UnitRange{Int}} end -Spectrum(energies, states, subs) = Spectrum(energies, states, subs, copy(subs)) +Spectrum(energies::AbstractVector{T}, states) where {T} = Spectrum(energies, states, zero(SVector{0,real(T)})) + +function Spectrum(energies, states, basevert::SVector{<:Any,T}) where {T} + energies´ = maybereal(energies, T) + subs = collect(approxruns(energies´)) + for rng in subs + orthonormalize!(view(states, :, rng)) + end + return Spectrum(energies´, states, basevert, subs) +end + +maybereal(energies, ::Type{T}) where {T<:Real} = T.(real.(energies)) +maybereal(energies, ::Type{T}) where {T<:Complex} = T.(energies) -function Base.show(io::IO, s::Spectrum{C,T}) where {C,T} +function Base.show(io::IO, s::Spectrum{D,C,T}) where {D,C,T} i = get(io, :indent, "") print(io, -"$(i)Spectrum{$C,$T}: spectrum of a 0D Hamiltonian +"$(i)Spectrum{$D,$C,$T}: spectrum of a $(D)D Hamiltonian $i Orbital type : $C $i Energy type : $T $i Energy range : $(extrema(real, s.energies)) +$i Bloch/params : $(Tuple(s.basevert)) $i Eigenpairs : $(length(s.energies)) $i Subspaces : $(length(s.subs))") end @@ -98,14 +114,11 @@ independent copy. function spectrum(h; method = LinearAlgebraPackage(), transform = missing) diag = diagonalizer(h; method = method) (ϵk, ψk) = diag(()) - subs = collect(approxruns(ϵk)) - s = Spectrum(ϵk, ψk, subs) + s = Spectrum(ϵk, ψk) transform === missing || transform!(transform, s) return s end -Base.Tuple(s::Spectrum) = (s.energies, s.states) - """ transform!(f::Function, s::Spectrum) @@ -113,21 +126,15 @@ Transform the energies of `s` by applying `f` to them in place. """ transform!(f, s::Spectrum) = (map!(f, s.energies, s.energies); s) +# destructuring Base.iterate(s::Spectrum) = s.energies, Val(:states) Base.iterate(s::Spectrum, ::Val{:states}) = s.states, Val(:done) Base.iterate(::Spectrum, ::Val{:done}) = nothing Base.first(s::Spectrum) = s.energies Base.last(s::Spectrum) = s.states +Base.Tuple(s::Spectrum) = (s.energies, s.states) -_subspace(s::Spectrum, rngs) = _subspace.(Ref(s), rngs) - -function _subspace(s::Spectrum, rng::AbstractUnitRange) - ε = mean(j -> s.energies[j], rng) - ψs = view(s.states, :, rng) - Subspace(ε, ψs) -end - -Base.getindex(s::Spectrum, i::Int) = _subspace(s, s.subs[i]) +Base.getindex(s::Spectrum, i::Int) = subspace(s, s.subs[i]) Base.getindex(s::Spectrum, is::Union{AbstractUnitRange,AbstractVector}) = getindex.(Ref(s), is) Base.getindex(s::Spectrum; around) = get_around(s, around) @@ -135,72 +142,76 @@ get_around(s::Spectrum, ε0::Number) = get_around(s, ε0, 1) get_around(s::Spectrum, (ε0, n)::Tuple) = get_around(s, ε0, 1:n) function get_around(s::Spectrum, ε0::Number, which) - copy!(s.subs´, s.subs) - rngs = partialsort!(s.subs´, which, by = rng -> abs(s.energies[first(rng)] - ε0)) - return _subspace(s, rngs) + rngs = partialsort(s.subs, which, by = rng -> abs(s.energies[first(rng)] - ε0)) + return subspace(s, rngs) end +subspace(s::Spectrum, rngs) = subspace.(Ref(s), rngs) + +function subspace(s::Spectrum, rng::AbstractUnitRange) + ε = mean(j -> s.energies[j], rng) + ψs = view(s.states, :, rng) + Subspace(ε, ψs) +end + +nsubspaces(s::Spectrum) = length(s.subs) + +basis_slice_type(s::Spectrum) = typeof(view(s.states, :, 1:0)) +basis_block_type(s::Spectrum) = typeof(view(s.states, 1:0, 1:0)) + ###################################################################### -# BandMesh +# Band ###################################################################### -struct BandMesh{D´,T<:Number} # D´ is dimension of BaseMesh space plus one (energy) - verts::Vector{SVector{D´,T}} # vertices of the band mesh - degs::Vector{Int} # Vertex degeneracies - adjmat::SparseMatrixCSC{Bool,Int} # Undirected adjacency graph: both dest > src and dest < src - sinds::Vector{NTuple{D´,Int}} # indices of verts D´-tuples that form simplices +struct Band{D,C,T,S<:AbstractMatrix{C},S´<:AbstractMatrix,D´} # D´ is dimension of BaseMesh space plus one (energy) + basemesh::CuboidMesh{D,T,D´} + verts::Vector{SVector{D´,T}} + vbases::Vector{S} + vptrs::Array{UnitRange{Int},D} + adjmat::SparseMatrixCSC{Bool,Int} + simps::Vector{NTuple{D´,Int}} + sbases::Vector{NTuple{D´,S´}} + sptrs::Array{UnitRange{Int},D´} end -function Base.show(io::IO, mesh::BandMesh{D}) where {D} +Band(basemesh::CuboidMesh{D,T}, verts::Vector{SVector{D´,T´}}, args...) where {D,D´,T,T´} = + Band(CuboidMesh{D,T´}(basemesh), verts, args...) + +function Base.show(io::IO, b::Band{D}) where {D} i = get(io, :indent, "") print(io, -"$(i)BandMesh{$D}: mesh of a $(D)D band in $(D-1)D parameter space -$i Vertices : $(nvertices(mesh)) -$i Edges : $(nedges(mesh)) -$i Simplices : $(nsimplices(mesh))") +"$(i)Band{$D}: a band over a $(D)D parameter cuboid +$i Vertices : $(nvertices(b)) +$i Edges : $(nedges(b)) +$i Simplices : $(nsimplices(b))") end -nvertices(m::BandMesh) = length(m.verts) +nvertices(m::Band) = length(m.verts) -nedges(m::BandMesh) = div(nnz(m.adjmat), 2) +nedges(m::Band) = div(nnz(m.adjmat), 2) -nsimplices(m::BandMesh) = length(m.sinds) +nsimplices(m::Band) = length(m.simps) -vertices(m::BandMesh) = m.verts +vertices(m::Band) = m.verts edges(adjmat, src) = nzrange(adjmat, src) edgedest(adjmat, edge) = rowvals(adjmat)[edge] -edgevertices(m::BandMesh) = - ((vsrc, m.verts[edgedest(m.adjmat, edge)]) for (i, vsrc) in enumerate(m.verts) for edge in edges(m.adjmat, i)) +edgevertices(b::Band) = + ((vsrc, b.verts[edgedest(b.adjmat, edge)]) for (i, vsrc) in enumerate(b.verts) for edge in edges(b.adjmat, i)) -degeneracy(m::BandMesh, i) = m.degs[i] +degeneracy(m::Band, i) = degeneracy(m.vbases[i]) -transform!(f::Function, m::BandMesh) = (map!(f, vertices(m), vertices(m)); m) +transform!(f::Function, m::Band) = (map!(f, vertices(m), vertices(m)); m) ####################################################################### # Bandstructure ####################################################################### -struct SimplexIndexer{D,T} - basemesh::CuboidMesh{D,T} - sptrs::Array{UnitRange{Int},D} # range of indices of sverts and svecs for each simplex CartesianIndex in base mesh -end - -struct Bandstructure{D,C,T,S<:AbstractMatrix{C},E,D´,B<:BandMesh{D´,T},M<:Diagonalizer} # D is dimension of base mesh, D´ = D+1 +struct Bandstructure{D,C,T,B<:Band{D,C,T},M<:Diagonalizer} # D is dimension of base mesh, D´ = D+1 bands::Vector{B} # band meshes (vertices + adjacencies) - sverts::Vector{NTuple{D´,SVector{D´,T}}} # (base-coords..., energy) of each simplex vertex (groupings of bands.verts) - sbases::Vector{NTuple{D´,S}} # basis on each simplex vertex, possibly degenerate - sprojs::Vector{NTuple{D´,Matrix{E}}} # projection of basis on each simplex vertex to interpolate - indexers::Vector{SimplexIndexer{D,T}} # provides ranges of simplices above corresponding to a given base-mesh minicuboid diag::M # diagonalizer that can be used to add additional base-meshes for refinement end -function Bandstructure(bands, sverts, sbases::Vector{NTuple{D´,S}}, indexers, diag) where {D´,C,S<:AbstractMatrix{C}} - E = eltype(C) - sprojs = Vector{NTuple{D´,Matrix{E}}}(undef, length(sbases)) - return Bandstructure(bands, sverts, sbases, sprojs, indexers, diag) -end - function Base.show(io::IO, bs::Bandstructure) i = get(io, :indent, "") ioindent = IOContext(io, :indent => string(i, " ")) @@ -214,7 +225,116 @@ end Base.summary(::Bandstructure{D}) where {D} = "Bandstructure{$D}: bands of a $(D)D Hamiltonian" -# API # +nvertices(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nvertices, bands(bs)) + +nedges(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nedges, bands(bs)) + +nsimplices(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nsimplices, bands(bs)) + +nbands(bs::Bandstructure) = length(bands(bs)) + +""" + bands(bs::Bandstructure[, i]) + +Return a `bands::Vector{Band}` of all the bands in `bs`, or `bands[i]` if `i` is given. + +# See also + `bandstructure` +""" +bands(bs::Bandstructure) = bs.bands +bands(bs::Bandstructure, i) = bs.bands[i] + +""" + vertices(bs::Bandstructure, i) + +Return the vertices `(k..., ϵ)` of the i-th band in `bs`, in the form of a +`Vector{SVector{L+1}}`, where `L` is the lattice dimension. +""" +vertices(bs::Bandstructure, i) = vertices(bands(bs, i)) + +""" + transform!(f::Function, b::Bandstructure) + +Transform the energies of all bands in `b` by applying `f` to them in place. + + transform!((fk, fε), b::Bandstructure) + +Transform Bloch phases and energies of all bands in `b` by applying `fk` and `fε` to them in +place, respectively. If any of them is `missing`, it will be ignored. + +""" +transform!(f, bs::Bandstructure) = transform!(sanitize_band_transform(f), bs) + +function transform!((fk, fε)::Tuple{Function,Function}, bs::Bandstructure) + for band in bands(bs) + vs = vertices(band) + for (i, v) in enumerate(vs) + vs[i] = SVector((fk(SVector(Base.front(Tuple(v))))..., fε(last(v)))) + end + end + return bs +end + +####################################################################### +# isometric and Brillouin zone points +####################################################################### +function isometric(h::Hamiltonian) + r = qr(bravais(h)).R + r = r * sign(r[1,1]) + ibr = inv(r') + return ϕs -> ibr * ϕs +end + +isometric(h::Hamiltonian{<:Any,L}, nodes) where {L} = _isometric(h, parsenode.(nodes, Val(L))) + +_isometric(h, pts::Tuple) = _isometric(h, [pts...]) + +function _isometric(h, pts::Vector) + br = bravais(h) + pts´ = map(p -> br * p, pts) + pathlength = pushfirst!(cumsum(norm.(diff(pts))), 0.0) + isometric = piecewise_mapping(pathlength) + return isometric +end + +nodeindices(nodes::NTuple{N,Any}) where {N} = ntuple(i -> i-1, Val(N)) + +piecewise_mapping(nodes, ::Val{N}) where {N} = piecewise_mapping(parsenode.(nodes, Val(N))) + +function piecewise_mapping(pts) + N = length(pts) # could be a Tuple or a different container + mapping = x -> begin + x´ = clamp(only(x), 0, N-1) + i = min(floor(Int, x´), N-2) + 1 + p = pts[i] + (x´ - i + 1) * (pts[i+1] - pts[i]) + return p + end + return mapping +end + +parsenode(pt::SVector, ::Val{L}) where {L} = padright(pt, Val(L)) +parsenode(pt::Tuple, val) = parsenode(SVector(float.(pt)), val) + +function parsenode(node::Symbol, val) + pt = get(BZpoints, node, missing) + pt === missing && throw(ArgumentError("Unknown Brillouin zone point $pt, use one of $(keys(BZpoints))")) + pt´ = parsenode(pt, val) + return pt´ +end + +const BZpoints = + ( Γ = (0,) + , X = (pi,) + , Y = (0, pi) + , Z = (0, 0, pi) + , K = (2pi/3, -2pi/3) + , K´ = (4pi/3, 2pi/3) + , M = (pi, 0) + ) + +####################################################################### +# bandstructure API +####################################################################### """ bandstructure(h::Hamiltonian; subticks = 13, kw...) @@ -257,7 +377,7 @@ Curried form of the above equivalent to `bandstructure(h[, mesh]; kw...)`. The default options are - (mapping = missing, minoverlap = 0.3, method = LinearAlgebraPackage(), transform = missing, showprogress = true) + (mapping = missing, method = LinearAlgebraPackage(), transform = missing, splitbands = true, showprogress = true) `mapping`: when not `missing`, `mapping = v -> p` is a function that map base mesh vertices `v` to Bloch phases and/or parameters `p`. The structure of `p` is whatever is accepted by @@ -266,9 +386,6 @@ Bloch phases. For `h::ParametricHamiltonian`, `p = (ϕs..., (; ps))` or `p = (ϕ combine Bloch phases `ϕs` and keyword parameters `ps` of `ph`. This allows to compute a bandstructure along a cut in the Brillouin zone/parameter space of `ph`, see examples below. -The option `minoverlap` determines the minimum overlap between eigenstates to connect -them into a common subband. - `method`: it is chosen automatically if unspecified, and can be one of the following method diagonalization function @@ -287,6 +404,8 @@ also do `transform -> (fφ, fε)` to transform also mesh vertices with fφ. Addi momenta, assuming they represent Bloch phases. This works both in full bandstructures and linecuts. +`splitbands`: split all bands into disconnected subbands. See also `splitbands!` + `showprogress`: indicate whether progress bars are displayed during the calculation # Indexing @@ -312,7 +431,7 @@ julia> bandstructure(h; subticks = 25, method = LinearAlgebraPackage()) Bandstructure{2}: collection of 2D bands Bands : 8 Element type : scalar (Complex{Float64}) - BandMesh{2}: mesh of a 2-dimensional manifold + Band{2}: mesh of a 2-dimensional manifold Vertices : 625 Edges : 1776 @@ -320,7 +439,7 @@ julia> bandstructure(h, :Γ, :X, :Y, :Γ; subticks = (10,15,10)) Bandstructure{2}: collection of 1D bands Bands : 18 Element type : scalar (Complex{Float64}) - BandMesh{1}: mesh of a 1-dimensional manifold + Band{1}: mesh of a 1-dimensional manifold Vertices : 33 Edges : 32 @@ -329,7 +448,7 @@ julia> bandstructure(h, mesh((0, 2π); subticks = 13); mapping = φ -> (φ, 0)) Bandstructure{2}: collection of 1D bands Bands : 18 Element type : scalar (Complex{Float64}) - BandMesh{1}: mesh of a 1-dimensional manifold + Band{1}: mesh of a 1-dimensional manifold Vertices : 11 Edges : 10 @@ -339,127 +458,14 @@ julia> bandstructure(ph, mesh((0, 2π); subticks = 13); mapping = φ -> (φ, 0, Bandstructure{2}: collection of 1D bands Bands : 18 Element type : scalar (Complex{Float64}) - BandMesh{1}: mesh of a 1-dimensional manifold + Band{1}: mesh of a 1-dimensional manifold Vertices : 11 Edges : 10 ``` # See also - `cuboid`, `diagonalizer`, `bloch`, `parametric` -""" -bandstructure - -nvertices(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nvertices, bands(bs)) - -nedges(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nedges, bands(bs)) - -nsimplices(bs::Bandstructure) = isempty(bands(bs)) ? 0 : sum(nsimplices, bands(bs)) - -nbands(bs::Bandstructure) = length(bands(bs)) - -""" - bands(bs::Bandstructure[, i]) - -Return a `bands::Vector{BandMesh}` of all the bands in `bs`, or `bands[i]` if `i` is given. - -# See also - `bandstructure` -""" -bands(bs::Bandstructure) = bs.bands -bands(bs::Bandstructure, i) = bs.bands[i] - -""" - vertices(bs::Bandstructure, i) - -Return the vertices `(k..., ϵ)` of the i-th band in `bs`, in the form of a -`Vector{SVector{L+1}}`, where `L` is the lattice dimension. -""" -vertices(bs::Bandstructure, i) = vertices(bands(bs, i)) - -""" - transform!(f::Function, b::Bandstructure) - -Transform the energies of all bands in `b` by applying `f` to them in place. - - transform!((fk, fε), b::Bandstructure) - -Transform Bloch phases and energies of all bands in `b` by applying `fk` and `fε` to them in -place, respectively. If any of them is `missing`, it will be ignored. - + `cuboid`, `diagonalizer`, `bloch`, `parametric`, `splitbands!` """ -transform!(f, bs::Bandstructure) = transform!(sanitize_transform(f), bs) - -function transform!((fk, fε)::Tuple{Function,Function}, bs::Bandstructure) - for band in bands(bs) - vs = vertices(band) - for (i, v) in enumerate(vs) - vs[i] = SVector((fk(SVector(Base.front(Tuple(v))))..., fε(last(v)))) - end - alignnormals!(band.sinds, vs) - end - return bs -end - -####################################################################### -# isometric and Brillouin zone points -####################################################################### -function isometric(h::Hamiltonian) - r = qr(bravais(h)).R - r = r * sign(r[1,1]) - ibr = inv(r') - return ϕs -> ibr * ϕs -end - -isometric(h::Hamiltonian{<:Any,L}, nodes) where {L} = _isometric(h, parsenode.(nodes, Val(L))) - -_isometric(h, pts::Tuple) = _isometric(h, [pts...]) - -function _isometric(h, pts::Vector) - br = bravais(h) - pts´ = map(p -> br * p, pts) - pathlength = pushfirst!(cumsum(norm.(diff(pts))), 0.0) - isometric = piecewise_mapping(pathlength) - return isometric -end - -nodeindices(nodes::NTuple{N,Any}) where {N} = ntuple(i -> i-1, Val(N)) - -piecewise_mapping(nodes, ::Val{N}) where {N} = piecewise_mapping(parsenode.(nodes, Val(N))) - -function piecewise_mapping(pts) - N = length(pts) # could be a Tuple or a different container - mapping = x -> begin - x´ = clamp(only(x), 0, N-1) - i = min(floor(Int, x´), N-2) + 1 - p = pts[i] + (x´ - i + 1) * (pts[i+1] - pts[i]) - return p - end - return mapping -end - -parsenode(pt::SVector, ::Val{L}) where {L} = padright(pt, Val(L)) -parsenode(pt::Tuple, val) = parsenode(SVector(float.(pt)), val) - -function parsenode(node::Symbol, val) - pt = get(BZpoints, node, missing) - pt === missing && throw(ArgumentError("Unknown Brillouin zone point $pt, use one of $(keys(BZpoints))")) - pt´ = parsenode(pt, val) - return pt´ -end - -const BZpoints = - ( Γ = (0,) - , X = (pi,) - , Y = (0, pi) - , Z = (0, 0, pi) - , K = (2pi/3, -2pi/3) - , K´ = (4pi/3, 2pi/3) - , M = (pi, 0) - ) - -####################################################################### -# bandstructure building -####################################################################### function bandstructure(h::Hamiltonian{<:Any, L}; subticks = 13, kw...) where {L} base = cuboid(filltuple((-π, π), Val(L))...; subticks = subticks) return bandstructure(h, base; kw...) @@ -469,246 +475,221 @@ function bandstructure(h::Hamiltonian{<:Any,L}, node1, node2, nodes...; subticks allnodes = (node1, node2, nodes...) mapping´ = piecewise_mapping(allnodes, Val(L)) base = cuboid(nodeindices(allnodes); subticks = subticks) - transform´ = sanitize_transform(transform, h, allnodes) + transform´ = sanitize_band_transform(transform, h, allnodes) return bandstructure(h, base; mapping = mapping´, transform = transform´, kw...) end function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, basemesh::CuboidMesh; - transform = missing, showprogress = true, kw...) + transform = missing, splitbands = true, showprogress = true, kw...) diag = diagonalizer(h; kw...) - b = bandstructure(diag, basemesh, showprogress) - if transform !== missing - transform´ = sanitize_transform(transform, h) - transform!(transform´, b) - end + b = bandstructure(diag, basemesh, splitbands, showprogress) + transform === missing || transform!(sanitize_band_transform(transform, h), b) return b end function bandstructure(matrixf::Function, basemesh::CuboidMesh; - mapping = missing, transform = missing, showprogress = true) + mapping = missing, transform = missing, splitbands = true, showprogress = true) matrixf´ = wrapmapping(mapping, matrixf) dimh = size(samplematrix(matrixf´, basemesh), 1) diag = diagonalizer(matrixf´, dimh) - b = bandstructure(diag, basemesh, showprogress) + b = bandstructure(diag, basemesh, splitbands, showprogress) transform === missing || transform!(transform, b) return b end + @inline map_phiparams(mapping::Missing, basevertex) = sanitize_phiparams(basevertex) @inline map_phiparams(mapping::Function, basevertex) = sanitize_phiparams(mapping(basevertex...)) wrapmapping(mapping::Missing, matrixf::Function) = matrixf wrapmapping(mapping::Function, matrixf::Function) = basevertex -> matrixf(toSVector(mapping(basevertex...))) -sanitize_transform(::Missing, args...) = (identity, identity) -sanitize_transform(f::Function, args...) = (identity, f) -sanitize_transform(f::typeof(isometric), args...) = (isometric(args...), identity) -sanitize_transform((_,f)::Tuple{typeof(isometric),Function}, args...) = (isometric(args...), f) -sanitize_transform(fs::Tuple{Function,Function}, args...) = fs -sanitize_transform((_,f)::Tuple{Missing,Function}, args...) = (identity, f) -sanitize_transform((f,_)::Tuple{Function,Missing}, args...) = (f, identity) +sanitize_band_transform(::Missing, args...) = (identity, identity) +sanitize_band_transform(f::Function, args...) = (identity, f) +sanitize_band_transform(f::typeof(isometric), args...) = (isometric(args...), identity) +sanitize_band_transform((_,f)::Tuple{typeof(isometric),Function}, args...) = (isometric(args...), f) +sanitize_band_transform(fs::Tuple{Function,Function}, args...) = fs +sanitize_band_transform((_,f)::Tuple{Missing,Function}, args...) = (identity, f) +sanitize_band_transform((f,_)::Tuple{Function,Missing}, args...) = (f, identity) samplematrix(matrixf, basemesh) = matrixf(Tuple(first(vertices(basemesh)))) -function bandstructure(diag::Diagonalizer, basemesh::CuboidMesh, showprogress) - # Step 1/2 - Diagonalising: - subspaces = bandstructure_diagonalize(diag, basemesh, showprogress) - # Step 2/2 - Knitting bands: - bands, cuboidinds, linearinds = bandstructure_knit(diag, basemesh, subspaces, showprogress) - - sverts, sbases, sptrs = bandstructure_collect(subspaces, bands, cuboidinds) - - indexers = [SimplexIndexer(basemesh, sptrs)] - - return Bandstructure(bands, sverts, sbases, indexers, diag) -end - ####################################################################### -# bandstructure_diagonalize +# bandstructure building ####################################################################### -function bandstructure_diagonalize(diag, basemesh::CuboidMesh, showprogress = false) +function bandstructure(diag::Diagonalizer, basemesh::CuboidMesh, split, showprogress) + # Step 1/2 - Diagonalise: + spectra = bandstructure_diagonalize(diag, basemesh, showprogress) + # Step 2/2 - Knit bands: + bands = bandstructure_knit(basemesh, spectra, showprogress) + bs = Bandstructure(bands, diag) + split && (bs = splitbands(bs)) + return bs +end + +## Diagonalize bands step + +function bandstructure_diagonalize(diag, basemesh::CuboidMesh, showprogress) prog = Progress(length(basemesh), "Step 1/2 - Diagonalising: ") - subspaces = [build_subspaces(diag, vertex, showprogress, prog) for vertex in vertices(basemesh)] - return subspaces + spectra = [build_spectrum(diag, vertex, showprogress, prog) for vertex in vertices(basemesh)] + return spectra end -function build_subspaces(diag::Diagonalizer, vertex::SVector{E,T}, showprog, prog) where {E,T} +function build_spectrum(diag, vertex, showprog, prog) (ϵs, ψs) = diag(Tuple(vertex)) - subspaces = collect_subspaces(ϵs, ψs, T) - showprog && ProgressMeter.next!(prog; showvalues = ()) - return subspaces -end + spectrum = Spectrum(ϵs, ψs, vertex) + showprog && ProgressMeter.next!(prog) + return spectrum +end + +## Knit bands step +function bandstructure_knit(basemesh::CuboidMesh{D}, spectra::AbstractArray{SP}, showprog) where {D,C,T,SP<:Spectrum{D,C,T}} + simpitr = marchingsimplices(basemesh) # D+1 - dimensional iterator over simplices + + S = basis_slice_type(first(spectra)) + verts = SVector{D+1,T}[] + vbases = S[] + vptrs = Array{UnitRange{Int}}(undef, size(vertices(basemesh))) + simps = NTuple{D+1,Int}[] + sbases = NTuple{D+1,Matrix{C}}[] + sptrs = Array{UnitRange{Int}}(undef, size(simpitr)) + + cbase = eachindex(basemesh) + lbase = LinearIndices(cbase) + + # Collect vertices + for csrc in cbase + len = length(verts) + push_verts!((verts, vbases), spectra[csrc]) + vptrs[csrc] = len+1:length(verts) + end -####################################################################### -# bandstructure_knit -####################################################################### -struct BandLinearIndex - bandidx::Int - vertidx::Int -end - -Base.zero(::BandLinearIndex) = zero(BandLinearIndex) -Base.zero(::Type{BandLinearIndex}) = BandLinearIndex(0, 0) - -Base.iszero(b::BandLinearIndex) = iszero(b.bandidx) - -struct BandCuboidIndex{D} - baseidx::CartesianIndex{D} - colidx::Int -end - -Base.Tuple(ci::BandCuboidIndex) = (Tuple(ci.baseidx)..., ci.colidx) - -function bandstructure_knit(diag, basemesh::CuboidMesh{D,T}, subspaces::Array{Vector{S},D}, showprog = false) where {D,T,C,S<:Subspace{C}} - nverts = sum(length, subspaces) - prog = Progress(nverts, "Step 2/2 - Knitting bands: ") - - bands = BandMesh{D+1,T}[] - pending = Tuple{BandCuboidIndex{D},BandCuboidIndex{D}}[] # pairs of neighboring vertex indices src::IT, dst::IT - linearinds = [zeros(BandLinearIndex, length(ss)) for ss in subspaces] # 0 == unclassified, > 0 vertex index - cuboidinds = BandCuboidIndex{D}[] # cuboid indices of processed vertices - I = Int[]; J = Int[] # To build adjacency matrices - P = real(eltype(C)) # type of projections between states - maxsubs = maximum(length, subspaces) - projinds = Vector{Tuple{P,Int}}(undef, maxsubs) # Reusable list of projections for sorting - - bandidx = 0 - while true - bandidx += 1 - seedidx = next_unprocessed(linearinds, subspaces) - seedidx === nothing && break - resize!(pending, 1) - resize!(I, 0) - resize!(J, 0) - pending[1] = (seedidx, seedidx) # source CartesianIndex for band search, with no originating vertex - bandmesh = knit_band(bandidx, basemesh, subspaces, diag.minoverlap, pending, cuboidinds, linearinds, I, J, projinds, showprog, prog) - iszero(nsimplices(bandmesh)) || push!(bands, bandmesh) + # Store subspace projections in vertex adjacency matrix + prog = Progress(length(cbase), "Step 2/2 - Knitting bands: ") + S´ = basis_block_type(first(spectra)) + I, J, V = Int[], Int[], S´[] + for csrc in cbase + for cdst in neighbors_forward(basemesh, csrc) + push_adjs!((I, J, V), spectra[csrc], spectra[cdst], vptrs[csrc], vptrs[cdst]) + end + showprog && ProgressMeter.next!(prog) + end + n = length(verts) + adjprojs = sparse(I, J, V, n, n) + adjmat = sparsealiasbool(adjprojs) + + # Build simplices from stable cycles around base simplices + buffers = NTuple{D+1,Int}[], NTuple{D+1,Int}[] + emptybases = filltuple(fill(zero(C), 0, 0), Val(D+1)) # sentinel bases for deg == 1 simplices + for (csimp, vs) in zip(CartesianIndices(simpitr), simpitr) # vs isa NTuple{D+1,CartesianIndex{D}} + len = length(simps) + ranges = getindex.(Ref(vptrs), vs) + push_simps!((simps, sbases, verts), buffers, ranges, adjprojs, emptybases) + sptrs[csimp] = len+1:length(simps) end - return bands, cuboidinds, linearinds + bands = [Band(basemesh, verts, vbases, vptrs, adjmat, simps, sbases, sptrs)] + return bands end -function next_unprocessed(linearinds, subspaces) - ci = CartesianIndices(linearinds) - @inbounds for (n, vs) in enumerate(linearinds), i in eachindex(subspaces[n]) - iszero(vs[i]) && return BandCuboidIndex(ci[n], i) +function push_verts!((verts, vbases), spectrum) + for rng in spectrum.subs + vert = SVector(Tuple(spectrum.basevert)..., spectrum.energies[first(rng)]) + basis = view(spectrum.states, :, rng) + push!(verts, vert) + push!(vbases, basis) end return nothing end -function knit_band(bandidx, basemesh::CuboidMesh{D,T}, subspaces, minoverlap, pending, cuboidinds, linearinds, I, J, projinds, showprog, prog) where {D,T} - verts = SVector{D+1,T}[] - degs = Int[] - vertcounter = 0 - while !isempty(pending) - src, dst = pop!(pending) - n, i = dst.baseidx, dst.colidx - n0, i0 = src.baseidx, src.colidx - # process dst only if unclassified (otherwise simply link) - if !iszero(linearinds[n][i]) - append_adjacent!(I, J, linearinds[n0][i0], linearinds[n][i]) - continue - end - - vert = vcat(vertex(basemesh, n), SVector(subspaces[n][i].energy)) - push!(verts, vert) - push!(degs, degeneracy(subspaces[n][i])) - push!(cuboidinds, dst) - vertcounter += 1 - linearinds[n][i] = BandLinearIndex(bandidx, vertcounter) - src == dst || append_adjacent!(I, J, linearinds[n0][i0], linearinds[n][i]) - showprog && ProgressMeter.next!(prog; showvalues = ()) - - subdst = subspaces[n][i] - deg = degeneracy(subdst) - for n´ in neighbors(basemesh, n) - deg == 1 && n´ == n0 && continue # Only if deg == 1 is no-backstep justified (think deg at BZ boundary) - sorted_valid_projections!(projinds, subspaces[n´], subdst, minoverlap, bandidx, linearinds[n´]) - cumdeg´ = 0 - for (p, i´) in projinds - i´ == 0 && break - push!(pending, (dst, BandCuboidIndex(n´, i´))) - cumdeg´ += degeneracy(subspaces[n´][i´]) - cumdeg´ >= deg && break # links on each column n´ = cumulated deg at most equal to deg links +function push_adjs!((I, J, V), ssrc, sdst, rsrc, rdst) + ψdst = sdst.states + ψsrc = ssrc.states + proj = ψdst' * ψsrc + proj´ = copy(proj') + for (is, rs) in enumerate(ssrc.subs) + srcdim = length(rs) + for (id, rd) in enumerate(sdst.subs) + crange = CartesianIndices((rd, rs)) + crange´ = CartesianIndices((rs, rd)) + rank = rankproj(proj, crange) + if !iszero(rank) + srcdim -= rank + srcdim < 0 && @warn("Unexpected band connectivity between $(ssrc.basevert) and $(sdst.basevert). Rank $rank in $(size(crange)) projector.") + append!(I, (rdst[id], rsrc[is])) + append!(J, (rsrc[is], rdst[id])) + append!(V, (view(proj, crange), view(proj´, crange´))) end end end - - adjmat = sparse(I, J, true) - - sinds = band_simplices(verts, adjmat) - - return BandMesh(verts, degs, adjmat, sinds) -end - -function append_adjacent!(I, J, msrc, mdst) - append!(I, (mdst.vertidx, msrc.vertidx)) - append!(J, (msrc.vertidx, mdst.vertidx)) return nothing end -function sorted_valid_projections!(projinds, subs::Vector{<:Subspace}, sub0::Subspace{C}, minoverlap, bandidx, linearindscol) where {C} - nsubs = length(subs) - realzero = zero(real(eltype(C))) - complexzero = zero(eltype(C)) - fill!(projinds, (realzero, 0)) - for (j, sub) in enumerate(subs) - bandidx´ = linearindscol[j].bandidx - bandidx´ == 0 || bandidx´ == bandidx || continue - p = proj_squared(sub.basis, sub0.basis, realzero, complexzero) - p > minoverlap && (projinds[j] = (p, j)) - end - sort!(projinds, rev = true, alg = Base.DEFAULT_UNSTABLE) - return projinds -end - -# non-allocating version of `sum(abs2, ψ' * ψ0)` -function proj_squared(ψ, ψ0, realzero, complexzero) - size(ψ, 1) == size(ψ0, 1) || throw(error("Internal error: eigenstates of different sizes")) - p = realzero - # project the smaller-dim subspace onto the larger-dim subspace - for j0 in axes(ψ0, 2), j in axes(ψ, 2) - p0 = complexzero - @simd for i0 in axes(ψ0, 1) - @inbounds p0 += dot(ψ[i0,j], ψ0[i0,j0]) - end - p += abs2(p0) - end - return p +# equivalent to r = round(Int, tr(m'm)) over crange, but if r > 0, must compute and count singular values +function rankproj(m, crange) + r = round(Int, sum(c -> abs2(m[c]), crange)) + (iszero(r) || minimum(size(crange)) == 1) && return r + s = svdvals(view(m, crange)) + r = count(s -> abs2(s) >= 0.5, s) + return r end -###################################################################### -# Simplices -###################################################################### -function band_simplices(vertices::Vector{SVector{D´,T}}, adjmat) where {D´,T} - D´ > 0 || throw(ArgumentError("Need a positive number of simplex vertices")) - nverts = length(vertices) - D´ == 1 && return Tuple.(1:nverts) - sinds = NTuple{D´,Int}[] - if nverts >= D´ - buffer = (NTuple{D´,Int}[], NTuple{D´,Int}[]) - for srcind in eachindex(vertices) - newsinds = vertex_simplices!(buffer, adjmat, srcind) - D´ > 2 && alignnormals!(newsinds, vertices) - append!(sinds, newsinds) +sparsealiasbool(s) = SparseMatrixCSC(s.m, s.n, s.colptr, s.rowval, fill(true, length(s.nzval))) + +function push_simps!((simps, sbases, verts), buffers, ranges, adjprojs, emptybases) + cycles = build_cycles!(buffers, ranges, adjprojs) + for cycle in cycles + done = false + projs = getindex.(Ref(adjprojs), shiftleft(cycle), cycle) + degs = degeneracy.(projs) + if all(isequal(1), degs) # optimization with sentinel emptybases + push!(simps, cycle) + push!(sbases, emptybases) + continue + end + d, start = findmin(degs) # minimum degeneracy of source + projs = shiftleft(projs, start - 1) + bases = shiftright(accumulate(orthomul, projs; init = I(d))) + while !done + if rankbasis(first(bases)) < d + basis_start = discard_zero_cols(first(bases)) + d = rankbasis(basis_start) + d == 0 && break + bases = shiftright(accumulate(orthomul, projs; init = basis_start)) + else + done = true + end + end + if done + push!(simps, cycle) + bases = shiftright(bases, start - 1) + push!(sbases, bases) end end - return sinds + return nothing end -# Add (greater) neighbors to last vertex of partials that are also neighbors of all members of partial, till N -function vertex_simplices!(buffer::Tuple{P,P}, adjmat, srcind) where {D´,P<:AbstractArray{<:NTuple{D´}}} - partials, partials´ = buffer - resize!(partials, 0) - push!(partials, padright((srcind,), Val(D´))) - for pass in 2:D´ - resize!(partials´, 0) +# build cycles of elements in ranges that are connected by adjprojs +function build_cycles!((partials, partials´), ranges::NTuple{D´}, adjprojs) where {D´} + resize!(partials, length(first(ranges))) + for (n, i) in enumerate(first(ranges)) + partials[n] = padright((i,), Val(D´)) + end + rows = rowvals(adjprojs) + for n in 1:D´ + empty!(partials´) for partial in partials - nextsrc = partial[pass - 1] - for edge in edges(adjmat, nextsrc), neigh in edgedest(adjmat, edge) - valid = neigh > nextsrc && isconnected(neigh, partial, adjmat) - valid || continue - newinds = tuplesplice(partial, pass, neigh) - push!(partials´, newinds) + if n == D´ + seed = first(partial) + nextrng = seed:seed + else + nextrng = ranges[n+1] + end + src = partial[n] + for ptr in nzrange(adjprojs, src) + dst = rows[ptr] + dst in nextrng || continue + partial´ = n == D´ ? partial : tuplesplice(partial, n+1, dst) + push!(partials´, partial´) end end partials, partials´ = partials´, partials @@ -716,111 +697,199 @@ function vertex_simplices!(buffer::Tuple{P,P}, adjmat, srcind) where {D´,P<:Abs return partials end -# equivalent to all(n -> n in neighbors(adjmat, neigh), partial) -function isconnected(neigh, partial, adjmat) - connected = all(partial) do ind - ind == 0 && return true - for edge in edges(adjmat, neigh), neigh´ in edgedest(adjmat, edge) - ind == neigh´ && return true +function discard_zero_cols(mat) + ncols = rankbasis(mat) + mat´ = similar(mat, size(mat, 1), ncols) + j = 0 + for col in eachcol(mat) + if !iszero(col) + j += 1 + mat´[:, j] .= col end - return false end - return connected + return mat´ end -function alignnormals!(simplices, vertices) - for (i, s) in enumerate(simplices) - volume = elementvolume(vertices, s) - volume < 0 && (simplices[i] = switchlast(s)) - end - return simplices +# cheap-rank for number of non-zero columns +rankbasis(m) = count(!iszero, eachcol(m)) + +# The minimum projected norm2 per src state (m column) is 0.5/dim(src) by default. +# 0.5 due to adjacency, and 1/dim(src) because the projection eigenstate has a minimum +# 1/√din(src) component over some src basis vector +function orthomul(m1, m2) + m = m2 * m1 + threshold = 0.5 / size(m, 2) + return orthonormalize!(m, threshold) end -# Project N-1 edges onto (N-1)-dimensional vectors to have a deterministic volume -elementvolume(verts, s::NTuple{N,Int}) where {N} = - elementvolume(hcat(ntuple(i -> padright(SVector(verts[s[i+1]] - verts[s[1]]), Val(N-1)), Val(N-1))...)) -elementvolume(mat::SMatrix{N,N}) where {N} = det(mat) +####################################################################### +# splitbands +####################################################################### +""" + splitbands(bs::Bandstructure) -switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N)) +Splits the bands in `bs` into disconnected subbands that share no vertices. See also `splitbands` option in `bandstructure`. -###################################################################### -# bandstructure_collect -###################################################################### -function bandstructure_collect(subspaces::Array{Vector{Subspace{C,T,S}},D}, bands, cuboidinds) where {C,T,S,D} - nsimps = isempty(bands) ? 0 : sum(nsimplices, bands) - sverts = Vector{NTuple{D+1,SVector{D+1,T}}}(undef, nsimps) - sbases = Vector{NTuple{D+1,S}}(undef, nsimps) - sptrs = fill(1:0, size(subspaces) .- 1) # assuming non-periodic basemesh - s0inds = Vector{CartesianIndex{D}}(undef, nsimps) # base cuboid index for reference vertex in simplex, for sorting - - scounter = 0 - ioffset = 0 - for band in bands - for s in band.sinds - scounter += 1 - let ioffset = ioffset # circumvent boxing, JuliaLang/#15276 - baseinds = (i -> cuboidinds[ioffset + i].baseidx).(s) - pbase = sortperm(SVector(baseinds)) - s0inds[scounter] = minimum(baseinds) # or equivalently, baseinds[first(pbase)] - sverts[scounter] = ntuple(i -> band.verts[s[pbase[i]]], Val(D+1)) - sbases[scounter] = ntuple(Val(D+1)) do i - c = cuboidinds[ioffset + s[pbase[i]]] - subspaces[c.baseidx][c.colidx].basis - end - end +# See also + `bandstructure` +""" +splitbands(b::Bandstructure) = Bandstructure(splitbands(b.bands), b.diag) + +function splitbands(bands::Vector{<:Band}) + bands´ = similar(bands, 0) + splitbands!.(Ref(bands´), bands) + return bands´ +end + +function splitbands!(bands, band) + bs, is, nbands = subgraphs(band.adjmat) + nbands == 1 && return push!(bands, band) + adjmats = split_adjmats(band.adjmat, bs, is, nbands) + verts, vbases, vptrs = split_list(band.verts, band.vbases, band.vptrs, bs, is, nbands) + simps, sbases, sptrs = split_list(band.simps, band.sbases, band.sptrs, bs, is, nbands) + for n in 1:nbands + newband = Band(band.basemesh, verts[n], vbases[n], vptrs[n], adjmats[n], simps[n], sbases[n], sptrs[n]) + nsimplices(newband) > 0 && push!(bands, newband) + end + return bands +end + +struct BandIndex + n::Int +end + +Base.:*(_, b::BandIndex) = b +Base.:*(b::BandIndex, _) = b +Base.:+(b1::BandIndex, b2::BandIndex) = BandIndex(min(b1.n, b2.n)) +Base.zero(b::Type{BandIndex}) = BandIndex(typemax(Int)) + +# for each row or column i in a, compute its subgraph bs[i] and the column/row index within that subgraph, is[i] +function subgraphs(a) + bs = BandIndex.(1:size(a, 2)) + bs´ = copy(bs) + total = sum(i -> i.n, bs) + done = false + while !done + mul!(bs´, a, bs, true, true) + bs, bs´ = bs´, bs + total´ = sum(i -> i.n, bs) + done = total == total´ + total = total´ + end + dict = Dict(1 => 1) + m = 1 + for b in bs + if !haskey(dict, b.n) + m += 1 + dict[b.n] = m end - ioffset += nvertices(band) end - - psimps = sortperm(s0inds; alg = Base.DEFAULT_UNSTABLE) - permute!(s0inds, psimps) - permute!(sverts, psimps) - permute!(sbases, psimps) - - for rng in equalruns(s0inds) - baseind = s0inds[first(rng)] - baseind in CartesianIndices(sptrs) || continue # TODO: should not be necessary - sptrs[baseind] = rng + is = similar(bs, Int) + ngraphs = m + counters = fill(0, ngraphs) + for (i, b) in enumerate(bs) + n = dict[b.n] + bs[i] = BandIndex(n) + counters[n] += 1 + is[i] = counters[n] end - - return sverts, sbases, sptrs + return bs, is, ngraphs +end + +function split_adjmats(adjmat, bs, is, nbands) + I0, J0, _ = findnz(adjmat) + Is = [Int[] for _ in 1:nbands] + Js = copy.(Is) + for (i, j) in zip(I0, J0) + n = bs[i].n + n == bs[j].n || throw(ArgumentError("Unexpected error in subgraphs function")) + push!(Is[n], is[i]) + push!(Js[n], is[j]) + end + nverts = [count(isequal(BandIndex(n)), bs) for n in 1:nbands] + adjmats = [sparse(Is[n], Js[n], true, nverts[n], nverts[n]) for n in 1:nbands] + return adjmats +end + +function split_list(vs, vbs, vptrs, bs, is, nbands) + vs´ = [similar(vs, 0) for _ in 1:nbands] + vbs´ = [similar(vbs, 0) for _ in 1:nbands] + vptrs´ = [similar(vptrs) for _ in 1:nbands] + counters = fill(0, nbands) + counters´ = copy(counters) + for (iptr, rng) in enumerate(vptrs) + for i in rng + n, v = get_vert_or_simp(vs[i], i, bs, is) + push!(vs´[n], v) + push!(vbs´[n], vbs[i]) + end + counters´ = length.(vs´) + for n in 1:nbands + vptrs´[n][iptr] = counters[n]+1:counters´[n] + end + counters, counters´ = counters´, counters + end + return vs´, vbs´, vptrs´ end +# return band index and vertex/simplex +get_vert_or_simp(v::SVector, i, bs, is) = bs[i].n, v +get_vert_or_simp(s::NTuple, i, bs, is) = bs[first(s)].n, getindex.(Ref(is), s) + ####################################################################### # Bandstructure indexing ####################################################################### Base.getindex(bs::Bandstructure, ϕs::Tuple; around = missing) = interpolate_bandstructure(bs, ϕs, around) -function interpolate_bandstructure(bs::Bandstructure{D,C,T}, ϕs, around) where {D,C,T} +function interpolate_bandstructure(bs::Bandstructure{D,C,T}, ϕs, around) where {D,C<:Complex,T} + subs = Subspace{D,C,T,Matrix{promote_type(C,T)}}[] + foreach(band -> interpolate_band!(subs, band, ϕs), bs.bands) + filter_around!(subs, around) + return subs +end + +function interpolate_band!(subs, band::Band{D}, ϕs) where {D} D == length(ϕs) || throw(ArgumentError("Bandstructure needs a NTuple{$D} of base coordinates for interpolation")) - found = false - inds = filltuple(0, Val(D)) - indexer = first(bs.indexers) - for outer indexer in bs.indexers - inds = find_basemesh_interval.(ϕs, indexer.basemesh.ticks) - found = !any(iszero, inds) - found && break - end - found || throw(ArgumentError("Cannot interpolate $ϕs within any of the bandstructure's base meshes")) - rng = indexer.sptrs[CartesianIndex(inds)] - subs = Subspace{C,T,Matrix{C}}[] - # Since we have sorted simplices to canonical base vertex order in bandstructure_collect - # we can avoid unncecessary simplices (and double matches with ϕs at a simplex boundary) - # by demanding equal vertices - simplexbase = find_basemesh_simplex(ϕs, bs, rng) - if simplexbase !== nothing - basevertices, dϕinds = simplexbase - for i in rng - sverts = bs.sverts[i] - basevertices´ = frontSVector.(sverts) - basevertices´ == basevertices || continue - sbases = bs.sbases[i] - sprojs = get_or_add_projection_basis!(bs, i) - push_interpolated_subspace!(subs, dϕinds, sverts, sbases, sprojs) + s = find_basemesh_simplex(promote(ϕs...), band.basemesh) + s === nothing && return subs + (ϕs, inds, dϕinds) = s + for i in band.sptrs[inds...] + si = band.simps[i] + ε0, εs = firsttail(last.(getindex.(Ref(band.verts), si))) + b0, bs = firsttail(getindex.(Ref(band.vbases), si)) + p0, ps = firsttail(band.sbases[i]) + energy = ε0 * (1 - sum(dϕinds)) + sum(εs .* dϕinds) + basis = isempty(p0) ? copy(b0) : b0 * p0 + lmul!(1 - sum(dϕinds), basis) + if isempty(p0) # deg == 1 sentinel optimization + ps´ = ntuple(i -> cis(angle(dot(bs[i], b0))), Val(D)) + for i in 1:D + basis .+= bs[i] .* (cis(angle(dot(bs[i], b0))) * dϕinds[i]) + end + else + foreach(i -> mul!(basis, bs[i], ps[i], dϕinds[i], true), 1:D) end + push!(subs, Subspace(energy, basis, ϕs)) end - subs´ = filter_around!(subs, around) - return subs´ + return subs +end + +# The simplex sought has unitperm such that reverse(dϕs[unitperm]) is sorted +function find_basemesh_simplex(ϕs::NTuple{D}, basemesh) where {D} + baseinds = find_basemesh_interval.(ϕs, basemesh.ticks) + any(iszero, baseinds) && return nothing + vertex0 = getindex.(basemesh.ticks, baseinds) + vertex1 = getindex.(basemesh.ticks, baseinds .+ 1) + dϕs = (ϕs .- vertex0) ./ (vertex1 .- vertex0) + dϕperm = Tuple(sortperm(SVector(dϕs), rev = true)) + baseperms = unitperms(basemesh) + i = findfirst(isequal(dϕperm), baseperms) + i === nothing && return nothing + sverts = SVector.(Tuple.(Base.tail(unitsimplex(basemesh, i)))) + smat = hcat(sverts...) + dϕinds = inv(smat) * SVector(dϕs) + inds = (baseinds..., i) + return ϕs, inds, dϕinds end function find_basemesh_interval(ϕ, ticks) @@ -831,77 +900,6 @@ function find_basemesh_interval(ϕ, ticks) return 0 end -function find_basemesh_simplex(ϕs, bs, rng) - for i in rng - basevertices = frontSVector.(bs.sverts[i]) - dϕinds = normalized_base_inds(ϕs, basevertices) - insimplex(dϕinds) && return basevertices, dϕinds - end - return nothing -end - -function normalized_base_inds(ϕs, baseverts) - dbase = tuple_minus_first(baseverts) - smat = hcat(dbase...) - dϕvec = SVector(ϕs) - first(baseverts) - dϕinds = inv(smat) * dϕvec - return dϕinds -end - -insimplex(dϕinds) = sum(dϕinds) <= 1 && all(i -> 0 <= i <= 1, dϕinds) - -add_projection_bases!(bs) = foreach(i -> get_or_add_projection_basis!(bs, i), eachindex(bs.sprojs)) - -function get_or_add_projection_basis!(bs, i) - isassigned(bs.sprojs, i) && return bs.sprojs[i] - bases = bs.sbases[i] - firstbasis = first(bases) - projbasis = firstbasis - for basis in Base.tail(bases) - size(basis, 2) < size(projbasis, 2) && (projbasis = basis) - end - sprojs = projection_basis.(bases, Ref(projbasis)) - bs.sprojs[i] = sprojs - return sprojs -end - -# Projects src onto dst, and finds an orthonormal basis of the projection subspace -function projection_basis(dst, src) - t = dst' * src # n´ x n matrix, so that dst * t spans the projection subspace - src === dst && return t # assumes orthonormal src == dst - n = size(src, 2) - n´ = size(dst, 2) - @inbounds for j in 1:n - col = view(t, :, j) - for j´ in 1:j-1 - col´ = view(t, :, j´) - r = dot(col´, col)/dot(col´, col´) - col .-= r .* col´ - end - normalize!(col) - end - return t -end - -function push_interpolated_subspace!(subs, dϕinds, verts, bases, projs) - ϵs = last.(verts) - dϵs = SVector(tuple_minus_first(ϵs)) - ϵ0 = first(ϵs) - energy = ϵ0 + dot(dϕinds, dϵs) - basis = interpolate_subspace_basis(dϕinds, bases, projs) - push!(subs, Subspace(energy, basis)) - return nothing -end - -function interpolate_subspace_basis(dϕinds, bases, projs) - ibasis = first(bases) * first(projs) - ibasis .*= (1 - sum(dϕinds)) - for i in 2:length(bases) - mul!(ibasis, bases[i], projs[i], dϕinds[i-1], 1) - end - return ibasis -end - filter_around!(ss, ::Missing) = sort!(ss; by = s -> s.energy, alg = Base.DEFAULT_UNSTABLE) filter_around!(ss, ε0::Number) = filter_around!(ss, ε0, 1) diff --git a/src/diagonalizer.jl b/src/diagonalizer.jl index 2759697e..f6b81c92 100644 --- a/src/diagonalizer.jl +++ b/src/diagonalizer.jl @@ -12,9 +12,8 @@ end (f::HamiltonianBlochFunctor)(vertex) = bloch!(f.matrix, f.h, map_phiparams(f.mapping, vertex)) -struct Diagonalizer{M<:AbstractDiagonalizeMethod,T<:Real,F} +struct Diagonalizer{M<:AbstractDiagonalizeMethod,F} method::M - minoverlap::T perm::Vector{Int} # reusable permutation vector matrixf::F # functor or function matrixf(φs) that produces matrices to be diagonalized end @@ -50,21 +49,22 @@ true # See also `bandstructure`, `spectrum` """ -function diagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}; method = LinearAlgebraPackage(), mapping = missing, minoverlap = 0.3) +function diagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}; method = LinearAlgebraPackage(), mapping = missing) matrix = similarmatrix(h, method_matrixtype(method, h)) matrixf = HamiltonianBlochFunctor(h, matrix, mapping) perm = Vector{Int}(undef, size(matrix, 2)) - return Diagonalizer(method, float(minoverlap), perm, matrixf) + return Diagonalizer(method, perm, matrixf) end -function diagonalizer(matrixf::Function, dimh; method = LinearAlgebraPackage(), minoverlap = 0.3) +function diagonalizer(matrixf::Function, dimh; method = LinearAlgebraPackage()) perm = Vector{Int}(undef, dimh) - return Diagonalizer(method, float(minoverlap), perm, matrixf) + return Diagonalizer(method, perm, matrixf) end @inline function (d::Diagonalizer)(φs) ϵ, ψ = diagonalize(d.matrixf(φs), d.method) issorted(ϵ, by = real) || sorteigs!(d.perm, ϵ, ψ) + fixphase!(ψ) return ϵ, ψ end @@ -78,6 +78,14 @@ function sorteigs!(perm, ϵ::AbstractVector, ψ::AbstractMatrix) return ϵ, ψ end +fixphase!(ψ::AbstractArray{<:Real}) = ψ +function fixphase!(ψ::AbstractArray{<:Complex}) + for col in eachcol(ψ) + col ./= cis(angle(mean(col))) + end + return ψ +end + function Base.show(io::IO, d::Diagonalizer) print(io, "Diagonalizer with method : $(summary(d.method))") end diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index bfd35350..fe1b1768 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -517,7 +517,7 @@ julia> h[(3,3)][[1,2],[1,2]] .= Ref(@SMatrix[1 2; 2 1]) [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im] [1.0+0.0im 2.0+0.0im; 2.0+0.0im 1.0+0.0im] ``` -# See also: +# See also `onsite`, `hopping`, `bloch`, `bloch!` """ hamiltonian(lat::AbstractLattice, ts...; orbitals = missing, kw...) = @@ -829,7 +829,7 @@ supercell that contains cells that are both in the supercell of `h1` and `h2` Equivalent of the above for `Superlattice`s -# See also: +# See also `|`, `xor` """ (Base.:&)(s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice}) = @@ -845,7 +845,7 @@ supercell that contains cells that are either in the supercell of `h1` or `h2` Equivalent of the above for `Superlattice`s -# See also: +# See also `&`, `xor` """ (Base.:|)(s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice}) = @@ -862,7 +862,7 @@ both Equivalent of the above for `Superlattice`s -# See also: +# See also `&`, `|` """ (Base.xor)(s1::Hamiltonian{<:Superlattice}, s2::Hamiltonian{<:Superlattice}) = @@ -1544,7 +1544,7 @@ julia> h = LatticePresets.honeycomb() |> hamiltonian(onsite(1) + hopping(2)) |> [2, 2] = 1.0+0.0im ``` -# See also: +# See also `bloch!`, `similarmatrix` """ bloch(ϕs, axis = 0) = h -> bloch(h, ϕs, axis) @@ -1602,7 +1602,7 @@ julia> bloch!(similarmatrix(ph, flatten), ph, (0, 0, (; α = 2))) [3, 3] = 0.0+0.0im ``` -# See also: +# See also `bloch`, `similarmatrix` """ bloch!(matrix, h::Hamiltonian, ϕs, axis = 0) = _bloch!(matrix, h, toSVector(ϕs), axis) diff --git a/src/iterators.jl b/src/iterators.jl index e6204ea7..0b3ca559 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -364,6 +364,7 @@ Base.length(s::SiteSublats) = nsites(s.u) ####################################################################### # Runs ####################################################################### +# iteration yields ranges of subsequent xs elements such that istogether on them gives true struct Runs{T,F} xs::Vector{T} istogether::F @@ -402,6 +403,7 @@ struct MarchingNeighbors{D,G} len::Int end +# c is a CartesianRange over vertices, and n is a CartesianIndex from which to find neighbors function marchingneighbors(c, n) rp = max(n, first(c)):min(n + oneunit(n), last(c)) rm = max(n - oneunit(n), first(c)):min(n, last(c)) @@ -410,6 +412,13 @@ function marchingneighbors(c, n) return MarchingNeighbors(n, gen, len) end +function marchingneighbors_forward(c, n) + rp = max(n, first(c)):min(n + oneunit(n), last(c)) + gen = Iterators.filter(!isequal(n), rp) + len = length(rp) - 2 + return MarchingNeighbors(n, gen, len) +end + Base.iterate(m::MarchingNeighbors, s...) = iterate(m.gen, s...) Base.IteratorSize(::MarchingNeighbors) = Base.HasLength() @@ -418,4 +427,64 @@ Base.IteratorEltype(::MarchingNeighbors) = Base.HasEltype() Base.eltype(::MarchingNeighbors{D}) where {D} = CartesianIndex{D} -Base.length(m::MarchingNeighbors) = m.len \ No newline at end of file +Base.length(m::MarchingNeighbors) = m.len + +####################################################################### +# MarchingSimplices +####################################################################### +struct MarchingSimplices{D,D´} + cinds::CartesianIndices{D´,NTuple{D´,Base.OneTo{Int}}} + simps::Vector{NTuple{D´,CartesianIndex{D}}} # simplices verts in unit cuboid + unitperms::Vector{NTuple{D,Int}} # unit vector order of simplices +end + +# c is a CartesianIndices over all vertices +function marchingsimplices(c::CartesianIndices{D}) where {D} + unitperms = marchingsimplices_unitperms(Val(D)) + edgeperms = map(edge -> unitvector.(CartesianIndex{D}, edge), unitperms) + simps = marchingsimplices_verts.(edgeperms) + srange = eachindex(simps) + refverts = ntuple(i -> Base.OneTo(c.indices[i][1:end-1]), Val(D)) + cinds = CartesianIndices((refverts..., srange)) + return MarchingSimplices(cinds, simps, unitperms) +end + +# indices of unit vectors in a unit cuboid, taken in any order +marchingsimplices_unitperms(::Val{D}) where {D} = permutations(ntuple(identity, Val(D))) + +# simplex vertices as the sum of unit vectors in any permutation +marchingsimplices_verts(edges::NTuple{D}) where {D} = + (zero(CartesianIndex{D}), edges...) |> cumsum |> orientsimp + +function orientsimp(sverts) + volume = det(hcat(SVector.(Tuple.(Base.tail(sverts)))...)) + sverts´ = ifelse(volume >= 0, sverts, switchlast(sverts)) + return sverts´ +end + +switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N)) + +function Base.iterate(m::MarchingSimplices, s...) + it = iterate(m.cinds, s...) + it === nothing && return nothing + c, s´ = it + t = Tuple(c) + nsimp = last(t) + refvert = CartesianIndex(Base.front(t)) + verts = Ref(refvert) .+ m.simps[nsimp] + return verts, s´ +end + +Base.IteratorSize(::MarchingSimplices{D}) where {D} = Base.HasShape{D}() + +Base.axes(m::MarchingSimplices, i...) = axes(m.cinds, i...) + +Base.size(m::MarchingSimplices, i...) = size(m.cinds, i...) + +Base.CartesianIndices(m::MarchingSimplices) = m.cinds + +Base.IteratorEltype(::MarchingSimplices) = Base.HasEltype() + +Base.eltype(::MarchingSimplices{D,D´}) where {D,D´} = NTuple{D´,CartesianIndex{D}} + +Base.length(m::MarchingSimplices) = prod(length.(axes(m))) \ No newline at end of file diff --git a/src/lattice.jl b/src/lattice.jl index b55b1fce..eb1d3113 100644 --- a/src/lattice.jl +++ b/src/lattice.jl @@ -171,7 +171,7 @@ Bravais{2,2,Float64} : set of 2 Bravais vectors in 2D space. Matrix : [1.0 3.0; 2.0 4.0] ``` -# See also: +# See also `lattice` """ bravais(lat::AbstractLattice) = lat.bravais.matrix @@ -261,7 +261,7 @@ Lattice{3,2,Float64} : 2D lattice in 3D space Sites : (1) --> 1 total per unit cell ``` -# See also: +# See also `LatticePresets`, `bravais`, `sublat`, `supercell`, `intracell` """ function lattice(ss::Sublat...; bravais = (), kw...) @@ -607,7 +607,7 @@ Superlattice{2,2,Float64,2} : 2D lattice in 2D space, filling a 2D supercell Supersites : 9 ``` -# See also: +# See also `unitcell`, `siteselector` """ supercell(v...; kw...) = lat -> supercell(lat, v...; kw...) @@ -841,7 +841,7 @@ Lattice{2,2,Float64} : 2D lattice in 2D space Sites : (9) --> 9 total per unit cell ``` -# See also: +# See also `supercell`, `siteselector` """ unitcell(v::Union{SMatrix,Tuple,SVector,Integer}...; kw...) = lat -> unitcell(lat, v...; kw...) diff --git a/src/mesh.jl b/src/mesh.jl index 0bb008c0..eca63f61 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -1,8 +1,18 @@ ###################################################################### # CuboidMesh ###################################################################### -struct CuboidMesh{D,T} +struct CuboidMesh{D,T,D´} ticks::NTuple{D,Vector{T}} + simpitr::MarchingSimplices{D,D´} +end + +function Base.show(io::IO, mesh::CuboidMesh{D}) where {D} + i = get(io, :indent, "") + print(io, +"$(i)CuboidMesh{$D}: a mesh of a $(D)D parameter cuboid +$i Ranges : $(extrema.(mesh.ticks)) +$i Axes ticks : $(length.(mesh.ticks)) +$i Simplices : $(size(mesh.simpitr)) -> $(length(mesh.simpitr))") end """ @@ -46,25 +56,35 @@ function _cuboid(subticks, axesticks::Vararg{Tuple,L}) where {L} end allaxisticks end - return CuboidMesh(allticks) + simpitr = marchingsimplices(CartesianIndices(eachindex.(allticks))) + return CuboidMesh(allticks, simpitr) end Base.eachindex(mesh::CuboidMesh) = CartesianIndices(eachindex.(mesh.ticks)) +Base.getindex(mesh::CuboidMesh, n::CartesianIndex) = SVector(getindex.(mesh.ticks, Tuple(n))) +Base.getindex(mesh::CuboidMesh, i...) = mesh[eachindex(mesh)[i...]] + Base.size(mesh::CuboidMesh, i...) = size(eachindex(mesh), i...) Base.length(mesh::CuboidMesh) = prod(length.(mesh.ticks)) -vertex(mesh::CuboidMesh, n::CartesianIndex) = SVector(getindex.(mesh.ticks, Tuple(n))) - vertices(mesh::CuboidMesh) = (SVector(v) for v in Iterators.product(mesh.ticks...)) nvertices(mesh::CuboidMesh) = length(vertices(mesh)) neighbors(mesh::CuboidMesh, n::CartesianIndex) = marchingneighbors(eachindex(mesh), n) -function neighbors(mesh::CuboidMesh, j::Int) - c = eachindex(mesh) - l = LinearIndices(c) - return (l[i] for i in neighbors(mesh, c[j])) -end \ No newline at end of file +neighbors_forward(mesh::CuboidMesh, n::CartesianIndex) = marchingneighbors_forward(eachindex(mesh), n) + +# function neighbors(mesh::CuboidMesh, j::Int) +# c = eachindex(mesh) +# l = LinearIndices(c) +# return (l[i] for i in neighbors(mesh, c[j])) +# end + +marchingsimplices(m::CuboidMesh) = m.simpitr + +unitperms(m::CuboidMesh) = m.simpitr.unitperms + +unitsimplex(m::CuboidMesh, i) = m.simpitr.simps[i] \ No newline at end of file diff --git a/src/model.jl b/src/model.jl index 53a9297b..7ee124f3 100644 --- a/src/model.jl +++ b/src/model.jl @@ -18,7 +18,7 @@ closest pairs of sites in a lattice, irrespective of their sublattice. Obtain the actual nth-nearest-neighbot distance between sites in lattice `lat`. -# See also: +# See also `hopping` """ nrange(n::Int) = NeighborRange(n) @@ -613,7 +613,7 @@ Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space Coordination : 0.0 ``` -# See also: +# See also `hopping` """ onsite(o; kw...) = onsite(o, siteselector(; kw...)) @@ -726,7 +726,7 @@ Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space Coordination : 9.0 ``` -# See also: +# See also `onsite`, `nrange` """ function hopping(t; plusadjoint = false, range = nrange(1), kw...) @@ -837,7 +837,7 @@ do `@onsite!(...) + @hopping!(...)`). `ElementModifier`s are not model terms but transformations of an existing Hamiltonian that are meant to be applied sequentially (the order of application usually matters). -# See also: +# See also `@hopping!`, `parametric` """ macro onsite!(kw, f) @@ -864,7 +864,7 @@ do `@onsite!(...) + @hopping!(...)`). `ElementModifier`s are not model terms but transformations of an existing Hamiltonian that are meant to be applied sequentially (the order of application usually matters). -# See also: +# See also `@onsite!`, `parametric` """ macro hopping!(kw, f) @@ -998,7 +998,7 @@ KetModel{2}: model with 2 terms Coefficient : -1 ``` -# See also: +# See also `onsite`, `Vector`, `Matrix` """ ket(f; normalized = true, maporbitals::Bool = false, kw...) = KetModel(onsite(f; kw...), normalized, Val(maporbitals)) @@ -1053,7 +1053,7 @@ independently to each orbital. The remaining keywords `kw` are passed to `ket` a used to constrain the random amplitude to a subset of sites. `normalized`, however, is always `false`. -# See also: +# See also `ket` """ randomkets(n::Int, f = r -> cis(2pi*rand()); kw...) = diff --git a/src/plot_makie.jl b/src/plot_makie.jl index a92e23b0..cb93f33d 100644 --- a/src/plot_makie.jl +++ b/src/plot_makie.jl @@ -1,6 +1,9 @@ -using .Makie using GeometryBasics -import .Makie.AbstractPlotting: plot!, plot, to_value +using .Makie: AbstractPlotting +using .Makie.AbstractPlotting: to_value, RGBAf0, Vec3f0, FRect, @recipe, LineSegments, Theme, + lift, campixel, SceneSpace, Node, Axis, text!, on, mouse_selection, poly!, scale!, + translate!, linesegments!, mesh!, scatter!, meshscatter! +import .Makie.AbstractPlotting: plot!, plot """ plot(h::Hamiltonian) @@ -279,7 +282,7 @@ function plot!(plot::BandPlot2D) for (nb, color) in zip(bands, colors) band = bs.bands[nb] vertices = band.verts - simplices = band.sinds + simplices = band.simps linesegments!(plot, (t -> vertices[first(t)] => vertices[last(t)]).(simplices), linewidth = plot[:linethickness][], color = color) end @@ -306,7 +309,7 @@ function plot!(plot::BandPlot3D) for (nb, color) in zip(bandinds, colors) band = bs.bands[nb] vertices = band.verts - connectivity = [s[j] for s in band.sinds, j in 1:3] + connectivity = [s[j] for s in band.simps, j in 1:3] if isempty(connectivity) scatter!(plot, vertices, color = color) else diff --git a/src/plot_vegalite.jl b/src/plot_vegalite.jl index 39ef8a25..be26039f 100644 --- a/src/plot_vegalite.jl +++ b/src/plot_vegalite.jl @@ -180,10 +180,10 @@ function bandtable(b::Bandstructure{1,C,T}, (scalingx, scalingy), bandsiter) whe table = NT[] for (nb, band) in enumerate(bands(b)) verts = vertices(band) - sinds = band.sinds - isempty(sinds) && continue + simps = band.simps + isempty(simps) && continue s0 = (0, 0) - for s in sinds + for s in simps if first(s) == last(s0) || s0 == (0, 0) push!(table, (; x = verts[first(s)][1] * scalingx, y = verts[first(s)][2] * scalingy, band = nb, tooltip = degeneracy(band, first(s)))) else diff --git a/src/tools.jl b/src/tools.jl index 47e5ecc5..a191e575 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -18,8 +18,9 @@ toSVector(::Type{T}, ::Tuple{}) where {T} = SVector{0,T}() # Dynamic dispatch toSVector(v::AbstractVector) = SVector(Tuple(v)) -unitvector(::Type{SVector{L,T}}, i) where {L,T} = - SVector{L,T}(ntuple(j -> j == i ? one(T) : zero(T), Val(L))) +unitvector(::Type{SVector{L,T}}, i) where {L,T} = SVector{L,T}(unitvector(NTuple{L,T}, i)) +unitvector(::Type{CartesianIndex{L}}, i) where {L} = CartesianIndex(unitvector(NTuple{L,Int}, i)) +unitvector(::Type{NTuple{L,T}}, i) where {L,T} =ntuple(j -> j == i ? one(T) : zero(T), Val(L)) ensuretuple(s::Tuple) = s ensuretuple(s) = (s,) @@ -40,6 +41,9 @@ filltuple(x, ::NTuple{N,Any}) where {N} = ntuple(_ -> x, Val(N)) tuplesplice(s::NTuple{N,T}, ind, el) where {N,T} = ntuple(i -> i === ind ? T(el) : s[i], Val(N)) +shiftleft(s::NTuple{N,Any}, n = 1) where {N} = ntuple(i -> s[mod1(i+n, N)], Val(N)) +shiftright(s, n = 1) = shiftleft(s, -n) + tupleproduct(p1, p2) = tupleproduct(ensuretuple(p1), ensuretuple(p2)) tupleproduct(p1::NTuple{M,Any}, p2::NTuple{N,Any}) where {M,N} = ntuple(i -> (p1[1+fld(i-1, N)], p2[1+mod(i-1, N)]), Val(M * N)) @@ -63,9 +67,12 @@ end tuple_minus_first(t::Tuple{T,Vararg{T,D}}) where {D,T} = ntuple(i -> ifelse(t[i+1] ≈ t[1], zero(T), t[i+1] - t[1]), Val(D)) +firsttail(t::Tuple) = first(t), Base.tail(t) + frontSVector(s::SVector) = SVector(Base.front(Tuple(s))) mergetuples(ts...) = keys(merge(tonamedtuple.(ts)...)) + tonamedtuple(ts::Val{T}) where {T} = NamedTuple{T}(filltuple(0,T)) function deletemultiple_nocheck(dn::SVector{N}, axes::NTuple{M,Int}) where {N,M} @@ -205,35 +212,6 @@ function normalize_columns!(kmat::AbstractMatrix, cols) return kmat end -############################################################################################ - -# function pushapproxruns!(runs::AbstractVector{<:UnitRange}, list::AbstractVector{T}, -# offset = 0, degtol = sqrt(eps(real(T)))) where {T} -# len = length(list) -# len < 2 && return runs -# rmin = rmax = 1 -# prev = list[1] -# @inbounds for i in 2:len -# next = list[i] -# if abs(next - prev) < degtol -# rmax = i -# else -# rmin < rmax && push!(runs, (offset + rmin):(offset + rmax)) -# rmin = rmax = i -# end -# prev = next -# end -# rmin < rmax && push!(runs, (offset + rmin):(offset + rmax)) -# return runs -# end - -# function hasapproxruns(list::AbstractVector{T}, degtol = sqrt(eps(real(T)))) where {T} -# for i in 2:length(list) -# abs(list[i] - list[i-1]) < degtol && return true -# end -# return false -# end - eltypevec(::AbstractMatrix{T}) where {T<:Number} = T eltypevec(::AbstractMatrix{S}) where {N,T<:Number,S<:SMatrix{N,N,T}} = SVector{N,T} @@ -241,6 +219,24 @@ tuplesort((a,b)::Tuple{<:Number,<:Number}) = a > b ? (b, a) : (a, b) tuplesort(t::Tuple) = t tuplesort(::Missing) = missing +# Gram-Schmidt but with column normalization only when norm^2 >= threshold (otherwise zero!) +function orthonormalize!(m::AbstractMatrix, threshold = 0) + @inbounds for j in axes(m, 2) + col = view(m, :, j) + for j´ in 1:j-1 + col´ = view(m, :, j´) + norm2´ = dot(col´, col´) + iszero(norm2´) && continue + r = dot(col´, col)/norm2´ + col .-= r .* col´ + end + norm2 = real(dot(col, col)) + factor = ifelse(norm2 < threshold, zero(norm2), 1/sqrt(norm2)) + col .*= factor + end + return m +end + # Like copyto! but with potentially different tensor orders (adapted from Base.copyto!, see #33588) function copyslice!(dest::AbstractArray{T1,N1}, Rdest::CartesianIndices{N1}, src::AbstractArray{T2,N2}, Rsrc::CartesianIndices{N2}, by = identity) where {T1,T2,N1,N2} @@ -265,12 +261,25 @@ function append_slice!(dest::AbstractArray, src::AbstractArray{T,N}, Rsrc::Carte Rdest = (length(dest) + 1):(length(dest) + length(Rsrc)) resize!(dest, last(Rdest)) src′ = Base.unalias(dest, src) - @inbounds for (Is, Id) in zip(Rsrc, Rdest) + for (Is, Id) in zip(Rsrc, Rdest) @inbounds dest[Id] = src′[Is] end return dest end +permutations(ss::NTuple) = permutations!(typeof(ss)[], ss, ()) + +function permutations!(p, s1, s2) + for (i, s) in enumerate(s1) + permutations!(p, delete(s1, i), (s2..., s)) + end + return p +end + +permutations!(p, ::Tuple{}, s2) = push!(p, s2) + +delete(t::NTuple{N,Any}, i) where {N} = ntuple(j -> j < i ? t[j] : t[j+1], Val(N-1)) + ###################################################################### # convert a matrix/number block to a matrix/inlinematrix string ###################################################################### @@ -300,95 +309,6 @@ end numberstring(x) = _isreal(x) ? string(" ", real(x)) : _isimag(x) ? string(" ", imag(x), "im") : string(" ", x) -###################################################################### -# Permutations (taken from Combinatorics.jl) -###################################################################### - -struct Permutations{T} - a::T - t::Int -end - -Base.eltype(::Type{Permutations{T}}) where {T} = Vector{eltype(T)} - -Base.length(p::Permutations) = (0 <= p.t <= length(p.a)) ? factorial(length(p.a), length(p.a)-p.t) : 0 - -""" - permutations(a) -Generate all permutations of an indexable object `a` in lexicographic order. Because the number of permutations -can be very large, this function returns an iterator object. -Use `collect(permutations(a))` to get an array of all permutations. -""" -permutations(a) = Permutations(a, length(a)) - -""" - permutations(a, t) -Generate all size `t` permutations of an indexable object `a`. -""" -function permutations(a, t::Integer) - if t < 0 - t = length(a) + 1 - end - Permutations(a, t) -end - -function Base.iterate(p::Permutations, s = collect(1:length(p.a))) - (!isempty(s) && max(s[1], p.t) > length(p.a) || (isempty(s) && p.t > 0)) && return - nextpermutation(p.a, p.t ,s) -end - -function nextpermutation(m, t, state) - perm = [m[state[i]] for i in 1:t] - n = length(state) - if t <= 0 - return(perm, [n+1]) - end - s = copy(state) - if t < n - j = t + 1 - while j <= n && s[t] >= s[j]; j+=1; end - end - if t < n && j <= n - s[t], s[j] = s[j], s[t] - else - if t < n - reverse!(s, t+1) - end - i = t - 1 - while i>=1 && s[i] >= s[i+1]; i -= 1; end - if i > 0 - j = n - while j>i && s[i] >= s[j]; j -= 1; end - s[i], s[j] = s[j], s[i] - reverse!(s, i+1) - else - s[1] = n+1 - end - end - return (perm, s) -end - -# Taken from Combinatorics.jl -# TODO: This should really live in Base, otherwise it's type piracy -""" - factorial(n, k) - -Compute ``n!/k!``. -""" -function Base.factorial(n::T, k::T) where T<:Integer - if k < 0 || n < 0 || k > n - throw(DomainError((n, k), "n and k must be nonnegative with k ≤ n")) - end - f = one(T) - while n > k - f = Base.checked_mul(f, n) - n -= 1 - end - return f -end - -Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...) - ############################################################################################ ######## fast sparse copy # Revise after #33589 is merged ################################# ############################################################################################ diff --git a/test/test_bandstructure.jl b/test/test_bandstructure.jl index dcc03f5b..64df1fd0 100644 --- a/test/test_bandstructure.jl +++ b/test/test_bandstructure.jl @@ -40,6 +40,14 @@ using Arpack h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(10) b = bandstructure(h, :K, :K´, subticks = 7, method = ArpackPackage(sigma = 0.1im, nev = 12), showprogress = false) @test nbands(b) == 1 + + # number of simplices contant across BZ (computing their degeneracy with the projs at each vertex j) + h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(3) + nev = size(h, 1) + b = bandstructure(h, splitbands = false) + @test nbands(b) == 1 + band = bands(b, 1) + @test all(j -> all(rng -> sum(i -> degeneracy(band.sbases[i][j]), rng) == nev, band.sptrs), 1:3) end @testset "functional bandstructures" begin @@ -48,7 +56,8 @@ end hf((x,)) = bloch!(matrix, hc, (x, -x)) m = cuboid((0, 1)) b = bandstructure(hf, m, showprogress = false) - @test nbands(b) == 2 + @test nbands(b) == 2 + @test nsimplices(b) == 24 hc2 = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) hp2 = parametric(hc2, @hopping!((t; s) -> s*t)) @@ -56,7 +65,8 @@ end hf2((s, x)) = bloch!(matrix2, hp2(s = s), (x, x)) m2 = cuboid((0, 1), (0, 1)) b = bandstructure(hf2, m2, showprogress = false) - @test nbands(b) == 1 + @test nbands(b) == 1 + @test nsimplices(b) == 576 end @testset "bandstructures lifts & transforms" begin @@ -141,7 +151,7 @@ end h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(2) bs = bandstructure(h, subticks = 13) @test sum(degeneracy, bs[(1,2)]) == size(h,1) - @test_broken sum(degeneracy, bs[(0.2,0.3)]) == size(h,1) + @test sum(degeneracy, bs[(0.2,0.3)]) == size(h,1) end @testset "diagonalizer" begin @@ -152,5 +162,5 @@ end d = diagonalizer(h) @test first(d()) ≈ [-3,3] @test first(d(())) ≈ [-3,3] - @test d() == Tuple(spectrum(h)) + @test all(d() .≈ Tuple(spectrum(h))) end \ No newline at end of file From b447792c2dac8cd43923aa5c0c3bcbab1d81cb94 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 10 Dec 2020 14:53:12 +0100 Subject: [PATCH 2/2] up min julia version --- .github/workflows/ci.yml | 2 +- Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6054bb75..c6292c5a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -10,8 +10,8 @@ jobs: strategy: matrix: version: - - '1.4' - '1.5' + - 'nightly' os: - ubuntu-latest - macOS-latest diff --git a/Project.toml b/Project.toml index baebac1f..38df31ad 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ ProgressMeter = "1.2" Requires = "1.0" SpecialFunctions = "0.10, 1.0" StaticArrays = "0.12.3, 0.13" -julia = "1.4" +julia = "1.5" [extras] ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"