diff --git a/Project.toml b/Project.toml index c95b857a..e9d758d6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ authors = ["Pablo San-Jose"] version = "0.4.0" [deps] -DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -20,7 +19,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] -DualNumbers = "0.6" ExprTools = "^0.1" GeometryBasics = "^0.3" LinearMaps = "2.6" diff --git a/src/Quantica.jl b/src/Quantica.jl index 9a536b6f..5fe80d95 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -11,7 +11,7 @@ function __init__() end using StaticArrays, NearestNeighbors, SparseArrays, LinearAlgebra, OffsetArrays, - ProgressMeter, LinearMaps, Random, SpecialFunctions, DualNumbers + ProgressMeter, LinearMaps, Random, SpecialFunctions using ExprTools @@ -25,7 +25,7 @@ export sublat, bravais, lattice, dims, supercell, unitcell, ket, randomkets, hamiltonian, parametric, bloch, bloch!, similarmatrix, flatten, wrap, transform!, combine, - spectrum, bandstructure, mesh, isometric, + spectrum, bandstructure, cuboid, isometric, bands, vertices, energies, states, momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM, diff --git a/src/bandstructure.jl b/src/bandstructure.jl index 28e54de3..5d334ea0 100644 --- a/src/bandstructure.jl +++ b/src/bandstructure.jl @@ -63,37 +63,82 @@ Transform the energies of `s` by applying `f` to them in place. """ transform!(f, s::Spectrum) = (map!(f, s.energies, s.energies); s) +###################################################################### +# BandMesh +###################################################################### +struct BandMesh{D´,T<:Number} # D´ is dimension of BaseMesh space plus one (energy) + verts::Vector{SVector{D´,T}} + adjmat::SparseMatrixCSC{Bool,Int} # Undirected graph: both dest > src and dest < src + simpinds::Vector{NTuple{D´,Int}} +end + +function Base.show(io::IO, mesh::BandMesh{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))") +end + +nvertices(m::BandMesh) = length(m.verts) + +nedges(m::BandMesh) = div(nnz(m.adjmat), 2) + +nsimplices(m::BandMesh) = length(m.simpinds) + +vertices(m::BandMesh) = m.verts + +edges(adjmat, src) = nzrange(adjmat, src) + +# neighbors(adjmat::SparseMatrixCSC, src::Int) = view(rowvals(adjmat), 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)) + +transform!(f::Function, m::BandMesh) = (map!(f, vertices(m), vertices(m)); m) + ####################################################################### # Bandstructure ####################################################################### -struct Band{D,M,A<:AbstractVector{M},MD<:Mesh{D},S<:AbstractArray} - mesh::MD # Mesh with missing vertices removed - simplices::S # Tuples of indices of mesh vertices that define mesh simplices - states::A # Must be resizeable container to build & refine band - dimstates::Int # Needed to extract the state at a given vertex from vector `states` +struct Simplices{D´,T,S<:SubArray,D} + sverts::Vector{NTuple{D´,SVector{D´,T}}} + sstates::Vector{NTuple{D´,S}} + sptrs::Array{UnitRange{Int},D} # range of indices of sverts and svecs for each simplex CartesianIndex in base mesh end -function Band(mesh::Mesh{D}, states::AbstractVector{M}, dimstates::Int) where {M,D} - simps = simplices(mesh, Val(D)) - return Band(mesh, simps, states, dimstates) -end - -struct Bandstructure{D,M,B<:Band{D,M}} # D is dimension of ε + parameter space +struct Bandstructure{D,T,M<:CuboidMesh{D},D´,B<:BandMesh{D´,T},S<:Simplices{D´,T}} # D is dimension of base mesh, D´ = D+1 + base::M bands::Vector{B} + simplices::S end -function Base.show(io::IO, b::Bandstructure{D,M}) where {D,M} +function Base.show(io::IO, bs::Bandstructure{D,M}) where {D,M} i = get(io, :indent, "") ioindent = IOContext(io, :indent => string(i, " ")) - print(io, i, summary(b), "\n", -"$i Bands : $(length(b.bands)) -$i Element type : $(displayelements(M))") + print(io, i, summary(bs), "\n", +"$i Bands : $(length(bs.bands)) +$i Element type : $(displayelements(M)) +$i Vertices : $(nvertices(bs)) +$i Edges : $(nedges(bs)) +$i Simplices : $(nsimplices(bs))") end -Base.summary(b::Bandstructure{D,M}) where {D,M} = - "Bandstructure{$D}: collection of $(D-1)D bands in a $(D)D space" +Base.summary(::Bandstructure{D,M}) where {D,M} = + "Bandstructure{$D}: bands of a $(D)D Hamiltonian" # API # + +nvertices(bs::Bandstructure) = sum(nvertices, bands(bs)) + +nedges(bs::Bandstructure) = sum(nedges, bands(bs)) + +nsimplices(bs::Bandstructure) = sum(nsimplices, bands(bs)) + +nbands(bs::Bandstructure) = length(bands(bs)) + """ bands(bs::Bandstructure) @@ -109,8 +154,6 @@ Return the vertices `(k..., ϵ)` of the i-th band in `bs`, in the form of a """ vertices(bs::Bandstructure, i) = vertices(bands(bs)[i]) -vertices(b::Band) = vertices(b.mesh) - """ energies(b::Bandstructure) @@ -129,7 +172,7 @@ Return the states of each vertex of the i-th band in `bs`, in the form of a `Mat """ states(bs::Bandstructure, i) = states(bands(bs)[i]) -states(b::Band) = reshape(b.states, b.dimstates, :) +# states(b::Band) = reshape(b.statess, b.dimstates, :) """ transform!(f::Function, b::Bandstructure) @@ -150,29 +193,86 @@ function transform!((fk, fε)::Tuple{Function,Function}, bs::Bandstructure) for (i, v) in enumerate(vs) vs[i] = SVector((fk(SVector(Base.front(Tuple(v))))..., fε(last(v)))) end - alignnormals!(band.simplices, vs) + alignnormals!(band.simpinds, 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 ####################################################################### """ - bandstructure(h::Hamiltonian; points = 13, kw...) + bandstructure(h::Hamiltonian; subticks = 13, kw...) -Compute `bandstructure(h, mesh((-π,π)...; points = points); kw...)` using a mesh over `h`'s -full Brillouin zone with the specified `points` along each [-π,π] reciprocal axis. +Compute `bandstructure(h, mesh((-π,π)...; subticks = subticks); kw...)` using a mesh over `h`'s +full Brillouin zone with the specified `subticks` along each [-π,π] reciprocal axis. - bandstructure(h::Hamiltonian, nodes...; points = 13, kw...) + bandstructure(h::Hamiltonian, nodes...; subticks = 13, kw...) Create a linecut of a bandstructure of `h` along a polygonal line connecting two or more `nodes`. Each node is either a `Tuple` or `SVector` of Bloch phases, or a symbolic name for a Brillouin zone point (`:Γ`,`:K`, `:K´`, `:M`, `:X`, `:Y` or `:Z`). Each segment in the -polygon has the specified number of `points`. Different `points` per segments can be -specified with `points = (p1, p2...)`. +polygon has the specified number of `subticks`. Different `subticks` per segments can be +specified with `subticks = (p1, p2...)`. - bandstructure(h::Hamiltonian, mesh::Mesh; mapping = missing, kw...) + bandstructure(h::Hamiltonian, mesh::BandMesh; mapping = missing, kw...) Compute the bandstructure `bandstructure(h, mesh; kw...)` of Bloch Hamiltonian `bloch(h, ϕ)`, with `ϕ = v` taken on each vertex `v` of `mesh` (or `ϕ = mapping(v...)` if a `mapping` @@ -184,7 +284,7 @@ Compute the bandstructure of a `ph`. Unless all parameters have default values, is required between mesh vertices and Bloch/parameters for `ph`, see details on `mapping` below. - bandstructure(matrixf::Function, mesh::Mesh; kw...) + bandstructure(matrixf::Function, mesh::BandMesh; kw...) Compute the bandstructure of the Hamiltonian matrix `m = matrixf(ϕ)`, with `ϕ` evaluated on the vertices `v` of the `mesh`. Note that `ϕ` in `matrixf(ϕ)` is an unsplatted container. @@ -199,7 +299,7 @@ Curried form of the above equivalent to `bandstructure(h, [mesh]; kw...)`. The default options are - (mapping = missing, minoverlap = 0, method = LinearAlgebraPackage(), transform = missing) + (mapping = missing, minoverlap = 0.3, method = LinearAlgebraPackage(), transform = missing, showprogress = true) `mapping`: when not `missing`, `mapping = v -> p` is a function that map mesh vertices `v` to Bloch phases and/or parameters `p`. The structure of `p` is whatever is accepted by @@ -229,42 +329,44 @@ 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. +`showprogress`: indicate whether progress bars are displayed during the calculation + # Examples ```jldoctest julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(3); -julia> bandstructure(h; points = 25, method = LinearAlgebraPackage()) +julia> bandstructure(h; subticks = 25, method = LinearAlgebraPackage()) Bandstructure{2}: collection of 2D bands Bands : 8 Element type : scalar (Complex{Float64}) - Mesh{2}: mesh of a 2-dimensional manifold + BandMesh{2}: mesh of a 2-dimensional manifold Vertices : 625 Edges : 1776 -julia> bandstructure(h, :Γ, :X, :Y, :Γ; points = (10,15,10)) +julia> bandstructure(h, :Γ, :X, :Y, :Γ; subticks = (10,15,10)) Bandstructure{2}: collection of 1D bands Bands : 18 Element type : scalar (Complex{Float64}) - Mesh{1}: mesh of a 1-dimensional manifold + BandMesh{1}: mesh of a 1-dimensional manifold Vertices : 33 Edges : 32 -julia> bandstructure(h, mesh((0, 2π); points = 11); mapping = φ -> (φ, 0)) - # Equivalent to bandstructure(h, :Γ, :X; points = 11) +julia> bandstructure(h, mesh((0, 2π); subticks = 13); mapping = φ -> (φ, 0)) + # Equivalent to bandstructure(h, :Γ, :X; subticks = 13) Bandstructure{2}: collection of 1D bands Bands : 18 Element type : scalar (Complex{Float64}) - Mesh{1}: mesh of a 1-dimensional manifold + BandMesh{1}: mesh of a 1-dimensional manifold Vertices : 11 Edges : 10 julia> ph = parametric(h, @hopping!((t; α) -> t * α)); -julia> bandstructure(ph, mesh((0, 2π); points = 11); mapping = φ -> (φ, 0, (; α = 2φ))) +julia> bandstructure(ph, mesh((0, 2π); subticks = 13); mapping = φ -> (φ, 0, (; α = 2φ))) Bandstructure{2}: collection of 1D bands Bands : 18 Element type : scalar (Complex{Float64}) - Mesh{1}: mesh of a 1-dimensional manifold + BandMesh{1}: mesh of a 1-dimensional manifold Vertices : 11 Edges : 10 ``` @@ -272,35 +374,26 @@ Bandstructure{2}: collection of 1D bands # See also `mesh`, `bloch`, `parametric` """ -function bandstructure(h::Hamiltonian{<:Any, L}; points = 13, kw...) where {L} - m = mesh(filltuple((-π, π), Val(L))...; points = points) - return bandstructure(h, m; kw...) +function bandstructure(h::Hamiltonian{<:Any, L}; subticks = 13, kw...) where {L} + base = cuboid(filltuple((-π, π), Val(L))...; subticks = subticks) + return bandstructure(h, base; kw...) end -function bandstructure(h::Hamiltonian{<:Any,L}, node1, node2, nodes...; points = 13, transform = missing, kw...) where {L} +function bandstructure(h::Hamiltonian{<:Any,L}, node1, node2, nodes...; subticks = 13, transform = missing, kw...) where {L} allnodes = (node1, node2, nodes...) mapping´ = piecewise_mapping(allnodes, Val(L)) + base = cuboid(nodeindices(allnodes); subticks = subticks) transform´ = sanitize_transform(transform, h, allnodes) - return bandstructure(h, piecewise_mesh(allnodes, points); mapping = mapping´, transform = transform´, kw...) + return bandstructure(h, base; mapping = mapping´, transform = transform´, kw...) 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) - ) - -function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh; - method = LinearAlgebraPackage(), minoverlap = 0, mapping = missing, transform = missing) +function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, basemesh::CuboidMesh; + method = LinearAlgebraPackage(), minoverlap = 0.3, mapping = missing, transform = missing, showprogress = true) # ishermitian(h) || throw(ArgumentError("Hamiltonian must be hermitian")) matrix = similarmatrix(h, method_matrixtype(method, h)) - diag = diagonalizer(h, matrix, mesh, method, minoverlap, mapping) - matrixf(ϕsmesh) = bloch!(matrix, h, map_phiparams(mapping, ϕsmesh)) - b = _bandstructure(matrixf, matrix, mesh, diag) + diag = diagonalizer(matrix, method, minoverlap) + matrixf(vertex) = bloch!(matrix, h, map_phiparams(mapping, vertex)) + b = bandstructure(matrixf, basemesh, diag, showprogress) if transform !== missing transform´ = sanitize_transform(transform, h) transform!(transform´, b) @@ -308,20 +401,20 @@ function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh; return b end -function bandstructure(matrixf::Function, mesh::Mesh; - method = LinearAlgebraPackage(), minoverlap = 0, mapping = missing, transform = missing) +function bandstructure(matrixf::Function, basemesh::CuboidMesh; + method = LinearAlgebraPackage(), minoverlap = 0.3, mapping = missing, transform = missing, showprogress = true) matrixf´ = wrapmapping(mapping, matrixf) - matrix = _samplematrix(matrixf´, mesh) - diag = diagonalizer(matrixf´, matrix, mesh, method, minoverlap, missing) - b = _bandstructure(matrixf´, matrix, mesh, diag) + matrix = samplematrix(matrixf´, basemesh) + diag = diagonalizer(matrix, method, minoverlap) + b = bandstructure(matrixf´, basemesh, diag, showprogress) transform === missing || transform!(transform, b) return b end -@inline map_phiparams(mapping::Missing, ϕsmesh) = sanitize_phiparams(ϕsmesh) -@inline map_phiparams(mapping::Function, ϕsmesh) = sanitize_phiparams(mapping(ϕsmesh...)) +@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) = ϕsmesh -> matrixf(toSVector(mapping(ϕsmesh...))) +wrapmapping(mapping::Function, matrixf::Function) = basevertex -> matrixf(toSVector(mapping(basevertex...))) sanitize_transform(::Missing, args...) = (identity, identity) sanitize_transform(f::Function, args...) = (identity, f) @@ -331,173 +424,286 @@ 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) -_samplematrix(matrixf, mesh) = matrixf(Tuple(first(vertices(mesh)))) +samplematrix(matrixf, basemesh) = matrixf(Tuple(first(vertices(basemesh)))) + +function bandstructure(matrixf::Function, basemesh::CuboidMesh, diago::Diagonalizer, showprogress) + # Step 1/3 - Diagonalising: + subspaces, nverts = bandstructure_diagonalize(matrixf, basemesh, diago, showprogress) + # Step 2/3 - Knitting bands: + bands, cuboidinds, linearinds = bandstructure_knit(basemesh, diago, subspaces, nverts, showprogress) + # Step 3/3 - Collecting simplices: + simplices = bandstructure_collect(subspaces, bands, cuboidinds, showprogress) -function _bandstructure(matrixf::Function, matrix´::AbstractMatrix{M}, mesh::MD, d::Diagonalizer) where {M,D,T,MD<:Mesh{D,T}} - nϵ = 0 # Temporary, to be reassigned - ϵks = Matrix{T}(undef, 0, 0) # Temporary, to be reassigned - ψks = Array{M,3}(undef, 0, 0, 0) # Temporary, to be reassigned + return Bandstructure(basemesh, bands, simplices) +end + +####################################################################### +# bandstructure_diagonalize +####################################################################### +struct Subspace{C,T,S<:SubArray{C}} + energy::T + states::S +end - lenψ = size(matrix´, 1) - nk = nvertices(mesh) - # function to apply to eigenvalues when building bands (depends on momenta type) - by = _maybereal(T) +degeneracy(s::Subspace) = size(s.states, 2) - p = Progress(nk, "Step 1/2 - Diagonalising: ") - for (n, ϕs) in enumerate(vertices(mesh)) - matrix = matrixf(Tuple(ϕs)) +function bandstructure_diagonalize(matrixf::Function, basemesh::CuboidMesh{D,T}, diago::Diagonalizer{M,S}, showprogress = false) where {D,T,M,C,S<:SubArray{C}} + prog = Progress(length(basemesh), "Step 1/3 - Diagonalising: ") + subspaces = Array{Vector{Subspace{C,T,S}},D}(undef, size(basemesh)...) + nverts = 0 + for n in eachindex(basemesh) + matrix = matrixf(Tuple(vertex(basemesh, n))) # (ϵk, ψk) = diagonalize(Hermitian(matrix), d) ## This is faster (!) - (ϵk, ψk) = diagonalize(matrix, d.method) - resolve_degeneracies!(ϵk, ψk, ϕs, d.codiag) - if n == 1 # With first vertex can now know the number of eigenvalues... Reassign - nϵ = size(ϵk, 1) - ϵks = Matrix{T}(undef, nϵ, nk) - ψks = Array{M,3}(undef, lenψ, nϵ, nk) - end - copyslice!(ϵks, CartesianIndices((1:nϵ, n:n)), - ϵk, CartesianIndices((1:nϵ,)), by) - copyslice!(ψks, CartesianIndices((1:lenψ, 1:nϵ, n:n)), - ψk, CartesianIndices((1:lenψ, 1:nϵ))) - ProgressMeter.next!(p; showvalues = ()) + (ϵs, ψs) = diagonalize(matrix, diago.method) + subspaces[n] = collect_subspaces(ϵs, ψs, Subspace{C,T,S}) + nverts += length(subspaces[n]) + showprogress && ProgressMeter.next!(prog; showvalues = ()) end + return subspaces, nverts +end + +collect_subspaces(ϵs, ψs, ::Type{SS}) where {SS} = + SS[Subspace(ϵs[first(rng)], view(ψs, :, rng)) for rng in approxruns(ϵs)] + +####################################################################### +# 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 + +function bandstructure_knit(basemesh::CuboidMesh{D,T}, diago::Diagonalizer{M,S}, subspaces, nverts, showprog = false) where {D,T,M,S} + prog = Progress(nverts, "Step 2/3 - Knitting bands: ") - p = Progress(nϵ * nk, "Step 2/2 - Connecting bands: ") - pcounter = 0 - bands = Band{D+1,M,Vector{M},Mesh{D+1,T,Vector{SVector{D+1,T}}},Vector{NTuple{D+1,Int}}}[] - vertindices = zeros(Int, nϵ, nk) # 0 == unclassified, -1 == different band, > 0 vertex index - pending = Tuple{Int,CartesianIndex{2}}[] # (originating vertex index, (ϵ, k)) - dests = Int[]; srcs = Int[] # To build adjacency matrices - sizehint!(pending, nk) + 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(eltype(S))) # 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 - src = findfirst(iszero, vertindices) - src === nothing && break + bandidx += 1 + seedidx = next_unprocessed(linearinds, subspaces) + seedidx === nothing && break resize!(pending, 1) - resize!(dests, 0) - resize!(srcs, 0) - pending[1] = (0, src) # source CartesianIndex for band search, with no originating vertex - band = extract_band(mesh, ϵks, ψks, vertindices, d.minoverlap, pending, dests, srcs) - nverts = nvertices(band.mesh) - nverts > D && push!(bands, band) # avoid bands with no simplices - pcounter += nverts - ProgressMeter.update!(p, pcounter; showvalues = ()) + 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, diago.minoverlap, pending, cuboidinds, linearinds, I, J, projinds, showprog, prog) + push!(bands, bandmesh) end - return Bandstructure(bands) -end -_maybereal(::Type{<:Complex}) = identity -_maybereal(::Type{<:Real}) = real + return bands, cuboidinds, linearinds +end -####################################################################### -# extract_band -####################################################################### +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) + end + return nothing +end -function extract_band(kmesh::Mesh{D,T}, ϵks::AbstractArray{T}, ψks::AbstractArray{M}, vertindices, minoverlap, pending, dests, srcs) where {D,T,M} - lenψ, nϵ, nk = size(ψks) - kverts = vertices(kmesh) - states = eltype(ψks)[] - sizehint!(states, nk * lenψ) +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}[] - lenverts = 0 - sizehint!(verts, nk) - adjmat = SparseMatrixBuilder{Bool}() - srcidx = 0 # represents the index of the last added vertex (used to search for the nexts) + vertcounter = 0 while !isempty(pending) - origin, src = pop!(pending) # origin is the vertex index that originated this src, 0 if none (first) - srcidx = vertindices[src] - if srcidx != 0 - append_adjacent!(dests, srcs, origin, srcidx) + 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 - ϵ, k = Tuple(src) # src == CartesianIndex(ϵ::Int, k::Int) - vertex = vcat(kverts[k], SVector(ϵks[src])) - push!(verts, vertex) - srcidx = length(verts) - vertindices[ϵ, k] = srcidx - append_slice!(states, ψks, CartesianIndices((1:lenψ, ϵ:ϵ, k:k))) - append_adjacent!(dests, srcs, origin, srcidx) - added_vertices = 0 - for edgek in edges(kmesh, k) - k´ = edgedest(kmesh, edgek) - proj, ϵ´ = findmostparallel(ψks, k´, ϵ, k) - # if unclassified and sufficiently parallel add it to pending list - if proj >= minoverlap && !iszero(ϵ´) && iszero(vertindices[ϵ´, k´]) - push!(pending, (srcidx, CartesianIndex(ϵ´, k´))) - added_vertices += 1 + + vert = vcat(vertex(basemesh, n), SVector(subspaces[n][i].energy)) + push!(verts, vert) + 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) + found = false + for n´ in neighbors(basemesh, n) + deg == 1 && n´ == n0 && continue # Only if deg == 1 is this 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 + found = true end end - # In 1D we avoid backsteps, to keep nicely continuous bands - D == 1 && added_vertices == 0 && break end - for (i, vi) in enumerate(vertindices) - @inbounds vi > 0 && (vertindices[i] = -1) # mark as classified in a different band - end - adjmat = sparse(dests, srcs, true) - mesh = Mesh(verts, adjmat) - return Band(mesh, states, lenψ) + + adjmat = sparse(I, J, true) + + simpinds = band_simplices(verts, adjmat) + + return BandMesh(verts, adjmat, simpinds) end -function append_adjacent!(dests, srcs, origin, srcidx) - if origin != 0 && srcidx != 0 - append!(dests, (origin, srcidx)) - append!(srcs, (srcidx, origin)) - end +function append_adjacent!(I, J, msrc, mdst) + append!(I, (mdst.vertidx, msrc.vertidx)) + append!(J, (msrc.vertidx, mdst.vertidx)) return nothing end -function findmostparallel(ψks::Array{M,3}, destk, srcb, srck) where {M} - T = real(eltype(M)) - dimh, nϵ, nk = size(ψks) - maxproj = zero(T) - destb = 0 - @inbounds for nb in 1:nϵ - proj = zero(M) - for i in 1:dimh - proj += ψks[i, nb, destk]' * ψks[i, srcb, srck] +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(sub.states, sub0.states, 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(ψ, ψ0, realzero, complexzero) + size(ψ, 1) == size(ψ0, 1) || throw(error("Internal error: eigenstates of different sizes")) + p = realzero + for j0 in axes(ψ0, 2), j in axes(ψ, 2) + p0 = complexzero + @simd for i0 in axes(ψ0, 1) + @inbounds p0 += ψ[i0,j] ⋅ ψ0[i0,j0] end - absproj = T(abs(tr(proj))) - if maxproj <= absproj # must happen at least once - destb = nb - maxproj = absproj + p += abs2(p0) + end + return p +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) + simpinds = NTuple{D´,Int}[] + if nverts >= D´ + buffer = (NTuple{D´,Int}[], NTuple{D´,Int}[]) + for srcind in eachindex(vertices) + newsimps = vertex_simplices!(buffer, adjmat, srcind) + D´ > 2 && alignnormals!(newsimps, vertices) + append!(simpinds, newsimps) end end - return maxproj, destb + return simpinds end -####################################################################### -# resolve_degeneracies -####################################################################### -# Tries to make states continuous at crossings. Here, ϵ needs to be sorted -function resolve_degeneracies!(ϵ, ψ, ϕs, codiag) - issorted(ϵ, by = real) || sorteigs!(codiag.perm, ϵ, ψ) - hasapproxruns(ϵ, codiag.degtol) || return ϵ, ψ - ranges, ranges´ = codiag.rangesA, codiag.rangesB - resize!(ranges, 0) - pushapproxruns!(ranges, ϵ, 0, codiag.degtol) # 0 is an offset - for n in codiag.matrixindices - v = codiag.comatrix(ϕs, n) - resize!(ranges´, 0) - for (i, r) in enumerate(ranges) - subspace = view(ψ, :, r) - vsubspace = subspace' * v * subspace - veigen = eigen!(Hermitian(vsubspace)) - if hasapproxruns(veigen.values, codiag.degtol) - roffset = minimum(r) - 1 # Range offset within the ϵ vector - pushapproxruns!(ranges´, veigen.values, roffset, codiag.degtol) +# 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) + 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) end - subspace .= subspace * veigen.vectors end - ranges, ranges´ = ranges´, ranges - isempty(ranges) && break + partials, partials´ = partials´, partials end - return ψ + return partials end -# Could perhaps be better/faster using a generalized CoSort -function sorteigs!(perm, ϵ::Vector, ψ::Matrix) - resize!(perm, length(ϵ)) - p = sortperm!(perm, ϵ, by = real) - # permute!(ϵ, p) - sort!(ϵ, by = real) - Base.permutecols!!(ψ, p) - return ϵ, ψ +# 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 + end + return false + end + return connected end + +function alignnormals!(simplices, vertices) + for (i, s) in enumerate(simplices) + volume = elementvolume(vertices, s) + volume < 0 && (simplices[i] = switchlast(s)) + end + return simplices +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) + +switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N)) + +###################################################################### +# bandstructure_collect +###################################################################### +function bandstructure_collect(subspaces::Array{Vector{Subspace{C,T,S}},D}, bands, cuboidinds, showprog) where {C,T,S,D} + nsimplices = sum(band -> length(band.simpinds), bands) + prog = Progress(nsimplices, "Step 3/3 - Collecting simplices: ") + + sverts = Vector{NTuple{D+1,SVector{D+1,T}}}(undef, nsimplices) + sstates = Vector{NTuple{D+1,S}}(undef, nsimplices) + sptrs = fill(1:0, size(subspaces) .- 1) # assuming non-periodic basemesh + s0inds = Vector{CartesianIndex{D}}(undef, nsimplices) # base cuboid index for reference vertex in simplex, for sorting + + scounter = 0 + ioffset = 0 + for band in bands + for s in band.simpinds + scounter += 1 + let ioffset = ioffset # circumvent boxing, JuliaLang/#15276 + s0inds[scounter] = minimum(i -> cuboidinds[ioffset + i].baseidx, s) + sverts[scounter] = ntuple(i -> band.verts[s[i]], Val(D+1)) + sstates[scounter] = ntuple(Val(D+1)) do i + c = cuboidinds[ioffset + s[i]] + subspaces[c.baseidx][c.colidx].states + end + end + showprog && ProgressMeter.next!(prog; showvalues = ()) + end + ioffset += nvertices(band) + end + + p = sortperm(s0inds; alg = Base.DEFAULT_UNSTABLE) + permute!(s0inds, p) + permute!(sverts, p) + permute!(sstates, p) + + for rng in equalruns(s0inds) + sptrs[s0inds[first(rng)]] = rng + end + + return Simplices(sverts, sstates, sptrs) +end \ No newline at end of file diff --git a/src/diagonalizer.jl b/src/diagonalizer.jl index 153194f1..29592ec7 100644 --- a/src/diagonalizer.jl +++ b/src/diagonalizer.jl @@ -4,14 +4,13 @@ ####################################################################### abstract type AbstractDiagonalizeMethod end -struct Diagonalizer{S<:AbstractDiagonalizeMethod,C} - method::S - codiag::C - minoverlap::Float64 +struct Diagonalizer{M<:AbstractDiagonalizeMethod,S<:SubArray,T<:Real} + method::M + minoverlap::T + matviewtype::Type{S} end -diagonalizer(h, matrix, mesh, method, minoverlap, mapping) = - Diagonalizer(method, codiagonalizer(h, matrix, mesh, mapping), Float64(minoverlap)) +diagonalizer(matrix, method, minoverlap) = Diagonalizer(method, float(minoverlap), method_matviewtype(method, matrix)) ## Diagonalize methods ## @@ -85,159 +84,14 @@ function diagonalize(matrix::AbstractMatrix{M}, method::KrylovKitPackage) where return ϵ´, ψ´ end -### matrix types +### matrix/vector types similarmatrix(h, method::AbstractDiagonalizeMethod) = similarmatrix(h, method_matrixtype(method, h)) method_matrixtype(::LinearAlgebraPackage, h) = Matrix{blockeltype(h)} method_matrixtype(::AbstractDiagonalizeMethod, h) = flatten -####################################################################### -# shift and invert methods -####################################################################### - -function shiftandinvert(matrix::AbstractMatrix{Tv}, origin) where {Tv} - cols, rows = size(matrix) - # Shift away from real axis to avoid pivot point error in factorize - matrix´ = diagshift!(parent(matrix), origin + im) - F = factorize(matrix´) - lmap = LinearMap{Tv}((x, y) -> ldiv!(x, F, y), cols, rows, - ismutating = true, ishermitian = false) - return lmap -end - -function diagshift!(matrix::AbstractMatrix, origin) - matrix´ = parent(matrix) - vals = nonzeros(matrix´) - rowval = rowvals(matrix´) - for col in 1:size(matrix, 2) - found_diagonal = false - for ptr in nzrange(matrix´, col) - if col == rowval[ptr] - found_diagonal = true - vals[ptr] -= origin * I # To respect non-scalar eltypes - break - end - end - found_diagonal || throw(error("Sparse work matrix must include the diagonal. Possible bug in `similarmatrix`.")) - end - return matrix -end - -function invertandshift(ϵ::Vector{T}, origin) where {T} - ϵ´ = similar(ϵ, real(T)) - ϵ´ .= real(inv.(ϵ) .+ (origin + im)) # Caution: we assume a real spectrum - return ϵ´ -end - -####################################################################### -# Codiagonalizer -####################################################################### - -## Codiagonalizer -## Uses velocity operators along different directions. If not enough, use finite differences -## along mesh directions -struct Codiagonalizer{T,F<:Function} - comatrix::F - matrixindices::UnitRange{Int} - degtol::T - rangesA::Vector{UnitRange{Int}} # Prealloc buffer for degeneray ranges - rangesB::Vector{UnitRange{Int}} # Prealloc buffer for degeneray ranges - perm::Vector{Int} # Prealloc for sortperm! -end - -# mapping = missing is assumed when h is a Function that generates matrices, instead of a Hamiltonian or ParametricHamiltonian -function codiagonalizer(h, matrix::AbstractMatrix{T}, mesh, mapping) where {T} - sample_phiparams = map_phiparams(mapping, first(vertices(mesh))) - dirs = codiag_directions(sample_phiparams) - degtol = sqrt(eps(real(eltype(T)))) - delta = meshdelta(mesh) - delta = iszero(delta) ? degtol : delta - comatrix, matrixindices = codiag_function(h, matrix, mapping, dirs, delta) - return Codiagonalizer(comatrix, matrixindices, degtol, UnitRange{Int}[], UnitRange{Int}[], Int[]) -end - -function codiag_function(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, mapping, dirs, delta) - hdual = dual_if_parametric(h) - matrixdual = dualarray(matrix) - anyold = anyoldmatrix(matrix) - ndirs = length(dirs) - matrixindices = 1:(ndirs + ndirs + 1) - comatrix(meshϕs, n) = - if n <= ndirs # automatic differentiation using dual numbers - ϕsparams = dual_phisparams(map_phiparams(mapping, meshϕs), dirs[n]) - dualpart.(bloch!(matrixdual, hdual, ϕsparams)) - elseif n - ndirs <= ndirs # resort to finite differences - ϕsparams = delta_phisparams(map_phiparams(mapping, meshϕs), delta * dirs[n - ndirs]) - bloch!(matrix, h, ϕsparams) - else # use a fixed arbitrary matrix - anyold - end - return comatrix, matrixindices -end - -dual_if_parametric(ph::ParametricHamiltonian) = Dual(ph) -dual_if_parametric(h::Hamiltonian) = h - -# In the Function case we cannot know what directions to scan (arguments of matrixf). Also, -# we cannot be sure that dual numbers propagate. We thus restrict to finite differences in the mesh -# Note that mapping is already wrapped into matrixf in the calling bandstructure(::Function,...) -function codiag_function(matrixf::Function, matrix, mapping::Missing, meshdirs, delta) - anyold = anyoldmatrix(matrix) - ndirs = length(meshdirs) - matrixindices = 1:(ndirs + 1) - comatrix(meshϕs, n) = - if n <= ndirs # finite differences - matrixf(meshϕs + delta * meshdirs[n]) - else # use a fixed arbitrary matrix - anyold - end - return comatrix, matrixindices -end - -val_length(::SVector{N}, nt::NamedTuple) where {N} = Val(N+length(nt)) - -codiag_directions(phiparams) = codiag_directions(val_length(phiparams...), phiparams) - -function codiag_directions(::Val{L}, phiparams, direlements = 0:1) where {L} - directions = vec(SVector{L,Int}.(Iterators.product(ntuple(_ -> direlements, Val(L))...))) - mask_dirs!(directions, phiparams) - filter!(ispositive, directions) - unique!(directions) - sort!(directions, by = norm) - return directions -end - -# Zeros out any direction that cannot modify a param because it is not a number -function mask_dirs!(dirs::Vector{S}, pp) where {L,S<:SVector{L}} - valparams = values(last(pp)) - valids = valparams .!= maybe_sum.(valparams, 1) - n = length(first(pp)) - mask = SVector(ntuple(i -> i <= n || valids[i - n] , Val(L))) - map!(dir -> mask .* dir, dirs, dirs) - return dirs -end - -dual_phisparams(ϕs::SVector{N}, dir) where {N} = Dual.(ϕs, frontsvec(dir, Val(N))) -dual_phisparams(params::NamedTuple, dir) = NamedTuple{keys(params)}(maybe_dual.(values(params), tailtuple(dir, Val(length(params))))) -dual_phisparams((ϕs, params)::Tuple, dir) = (dual_phisparams(ϕs, dir), dual_phisparams(params, dir)) - -maybe_dual(param::Number, ε) = Dual(param, ε) -maybe_dual(param, ε) = param - -delta_phisparams(ϕs::SVector{N}, dir) where {N} = ϕs + frontsvec(dir, Val(N)) -delta_phisparams(params::NamedTuple, dir) = NamedTuple{keys(params)}(maybe_sum.(values(params), tailtuple(dir, Val(length(params))))) -delta_phisparams((ϕs, params)::Tuple, dir) = (delta_phisparams(ϕs, dir), delta_phisparams(params, dir)) - -maybe_sum(param::Number, ε) = param + ε -maybe_sum(param, ε) = param - -meshdelta(mesh::Mesh{<:Any,T}) where {T} = T(0.1) * norm(first(minmax_edge(mesh))) - -function anyoldmatrix(matrix::SparseMatrixCSC, rng = MersenneTwister(1)) - s = copy(matrix) - rand!(rng, nonzeros(s)) - return s -end +# Type of states in bandstructure. Should be Matrix view, because eigvecs are dense and we need to support degeneracies +method_matviewtype(::AbstractDiagonalizeMethod, ::AbstractMatrix{T}) where {T} = subarray_matrix_type(orbitaltype(T)) -anyoldmatrix(m::DenseArray, rng = MersenneTwister(1)) = rand!(rng, copy(m)) \ No newline at end of file +subarray_matrix_type(::Type{M}) where {M} = typeof(view(Matrix{eltype(M)}(undef, 2, 2), :, 1:0)) \ No newline at end of file diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index a9e0a5d3..9b101d21 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -314,18 +314,18 @@ check_unflatten_eltypes(v::AbstractVector{T}, h) where {T} = T === orbitaltype(h) || throw(ArgumentError("Eltype of desination array is inconsistent with Hamiltonian")) -## maybe_unflatten: call unflatten but only if we cannot do it without copying -maybe_unflatten(vflat, h) = _maybe_unflatten(vflat, h, orbitaltype(h), h.orbitals) +## unflatten_or_reinterpret: call unflatten but only if we cannot do it without copying +unflatten_or_reinterpret(vflat, h) = _unflatten_or_reinterpret(vflat, h, orbitaltype(h), h.orbitals) # source is already of the correct orbitaltype(h) -function _maybe_unflatten(v::AbstractVector{T}, h, ::Type{T}, orbs) where {T} +function _unflatten_or_reinterpret(v::AbstractVector{T}, h, ::Type{T}, orbs) where {T} check_unflatten_dst_dims(v, h) return v end # source can be reinterpreted, because the number of orbitals is the same M for all N sublattices -_maybe_unflatten(v::AbstractVector{T}, h, ::Type{S}, ::NTuple{N,NTuple{M}}) where {N,M,T<:Number,S<:SVector{M}} = +_unflatten_or_reinterpret(v::AbstractVector{T}, h, ::Type{S}, ::NTuple{N,NTuple{M}}) where {N,M,T<:Number,S<:SVector{M}} = reinterpret(SVector{M,T}, v) # otherwise call unflatten -_maybe_unflatten(v, h, S, orbs) = unflatten(v, h) +_unflatten_or_reinterpret(v, h, S, orbs) = unflatten(v, h) ####################################################################### # similarmatrix @@ -625,6 +625,9 @@ _orbitaltype(t::Type{SVector{1,Tv}}) where {Tv} = Tv orbitaltype(h::Hamiltonian{LA,L,M}) where {N,T,LA,L,M<:SMatrix{N,N,T}} = SVector{N,T} orbitaltype(h::Hamiltonian{LA,L,M}) where {LA,L,M<:Number} = M +orbitaltype(::Type{M}) where {M<:Number} = M +orbitaltype(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = SVector{N,T} + coordination(ham) = nhoppings(ham) / nsites(ham) function nhoppings(ham) @@ -721,12 +724,6 @@ SparseArrays.issparse(h::Hamiltonian{LA,L,M,A}) where {LA,L,M,A} = false Base.parent(h::Hamiltonian) = h -# Dual numbers # - -DualNumbers.Dual(h::Hamiltonian) = Hamiltonian(h.lattice, Dual.(h.harmonics), h.orbitals) - -DualNumbers.Dual(h::HamiltonianHarmonic) = HamiltonianHarmonic(h.dn, dualarray(h.h)) - # Iterators # function nonzero_indices(h::Hamiltonian, rowrange = 1:size(h, 1), colrange = 1:size(h, 2)) diff --git a/src/iterators.jl b/src/iterators.jl index b6239792..efad2c07 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -360,3 +360,62 @@ Base.IteratorEltype(::SiteSublats) = Base.HasEltype() Base.eltype(::SiteSublats) = Tuple{Int,Int} Base.length(s::SiteSublats) = nsites(s.u) + +####################################################################### +# Runs +####################################################################### +struct Runs{T,F} + xs::Vector{T} + istogether::F +end + +approxruns(xs::Vector{T}) where {T} = Runs(xs, (x, y) -> isapprox(x, y; atol = sqrt(eps(T)))) +equalruns(xs) = Runs(xs, ==) + +function last_in_run(xs::Vector{T}, i, istogether) where {T} + xi = xs[i] + for j in i:length(xs) + (j == length(xs) || !istogether(xs[j+1], xi)) && return j + end + return i +end + +# Base.iterate(s::Runs) = iterate(s, 1) +function Base.iterate(s::Runs, j = 1) + j > length(s.xs) && return nothing + lst = last_in_run(s.xs, j, s.istogether) + return j:lst, lst + 1 +end + +Base.IteratorSize(::Runs) = Base.SizeUnknown() + +Base.IteratorEltype(::Runs) = Base.HasEltype() + +Base.eltype(::Runs) = UnitRange{Int} + +####################################################################### +# MarchingNeighbors +####################################################################### +struct MarchingNeighbors{D,G} + n::CartesianIndex{D} + gen::G + len::Int +end + +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)) + gen = Iterators.filter(!isequal(n), Iterators.flatten((rp, rm))) + len = length(rp) + length(rm) - 2 + return MarchingNeighbors(n, gen, len) +end + +Base.iterate(m::MarchingNeighbors, s...) = iterate(m.gen, s...) + +Base.IteratorSize(::MarchingNeighbors) = Base.HasLength() + +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 diff --git a/src/mesh.jl b/src/mesh.jl index db414a6b..0bb008c0 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -1,363 +1,70 @@ ###################################################################### -# Mesh -####################################################################### - -abstract type AbstractMesh{D} end - -struct Mesh{D,T<:Number,V<:AbstractArray{SVector{D,T}}} <: AbstractMesh{D} # D is dimension of parameter space - vertices::V # Iterable vertex container with SVector{D,T} eltype - adjmat::SparseMatrixCSC{Bool,Int} # Undirected graph: both dest > src and dest < src -end - -# const Mesh{D,T} = Mesh{D,T,Vector{SVector{D,T}},Vector{Tuple{Int,Vararg{Int,D}}}} - -# Mesh{D,T}() where {D,T} = Mesh(SVector{D,T}[], sparse(Int[], Int[], Bool[]), NTuple{D+1,Int}[]) - -function Base.show(io::IO, mesh::Mesh{D}) where {D} - i = get(io, :indent, "") - print(io, -"$(i)Mesh{$D}: mesh of a $D-dimensional manifold -$i Vertices : $(nvertices(mesh)) -$i Edges : $(nedges(mesh))") +# CuboidMesh +###################################################################### +struct CuboidMesh{D,T} + ticks::NTuple{D,Vector{T}} end """ - mesh(minmaxaxes...; axes = 1.0 * I, points = 13) + cuboid(ticks...; subticks = 13) -Create a `Mesh` of L-dimensional marching-tetrahedra over a parallelepiped with axes given -by the columns of `axes`. The dimension `L` is given by the number of `minmaxaxes`, each of -the form `(min, max)`. The points along each axis are distributed between the corresponding -`min` and `max`. The number of points on each axis is given by `points`, or `points[i]` if -several are given. +Create a `CuboidMesh` of L-dimensional marching-tetrahedra over a cuboid aligned with the +Cartesian axes. The dimension `L` is given by the number of `ticks`, each of the form `(x₁, +x₂,...)`. The interval between `xⱼ` and `xⱼ₊₁` ticks in axis `i` are further subdivided to +have a number of subticks including endpoints. The number is `subticks` if `subticks` is an +`Integer`, `subticks[i]` if `subticks = (s₁, s₂,...)` or `subticks[i][j]` if `subticks = +((s₁₁, s₁₂,...), (s₂₁, s₂₂,...), ...)`. # Examples ```jldoctest -julia> buildmesh(mesh((-π, π), (0,2π); points = 25)) -Mesh{2}: mesh of a 2-dimensional manifold - Vertices : 625 - Edges : 1776 +julia> cuboid((-π, π), (0, 2π); subticks = 25) -julia> buildmesh(mesh((-π, π), (0,2π); points = (10,10))) -Mesh{2}: mesh of a 2-dimensional manifold - Vertices : 100 - Edges : 261 +julia> cuboid((-π, π), (0, 2π); subticks = (10, 10)) ``` # External links - Marching tetrahedra (https://en.wikipedia.org/wiki/Marching_tetrahedra) in Wikipedia """ -function mesh(minmaxaxes::Vararg{Tuple{Number,Number},D}; axes = 1.0 * I, points = 13) where {D} - ranges = ((b, r)->range(b...; length = r)).(minmaxaxes, points) - npoints = length.(ranges) - cs = CartesianIndices(ntuple(n -> 1:npoints[n], Val(D))) - ls = LinearIndices(cs) - csinner = CartesianIndices(ntuple(n -> 1:npoints[n]-1, Val(D))) - - # edge vectors for marching tetrahedra in D-dimensions (skip zero vector [first]) - uedges = [c for c in CartesianIndices(ntuple(_ -> 0:1, Val(D)))][2:end] - # tetrahedra built from the D unit-length uvecs added in any permutation - perms = permutations( - ntuple(i -> CartesianIndex(ntuple(j -> i == j ? 1 : 0, Val(D))), Val(D))) - utets = [cumsum(pushfirst!(perm, zero(CartesianIndex{D}))) for perm in perms] - - # We don't use generators because their non-inferreble eltype causes problems elsewhere - verts = [axes * SVector(getindex.(ranges, Tuple(c))) for c in cs] - - s = SparseMatrixBuilder{Bool}(length(cs), length(cs)) - for c in cs - for u in uedges - dest = c + u # dest > src - dest in cs && pushtocolumn!(s, ls[dest], true) - dest = c - u # dest < src - dest in cs && pushtocolumn!(s, ls[dest], true) +cuboid(ticks::Vararg{Tuple,L}; subticks = 13) where {L} = _cuboid(sanitize_subticks(subticks, ticks), ticks...) + +sanitize_subticks(st::NTuple{L,Any}, t::NTuple{L,Any}) where {L} = _sanitize_subticks.(st, t) +sanitize_subticks(st, t) = _sanitize_subticks.(Ref(st), t) +_sanitize_subticks(st::Number, t::Tuple{Vararg{Number}}) = Base.tail(_sanitize_subticks.(st, t)) +_sanitize_subticks(st::Number, ::Number) = Int(st) +_sanitize_subticks(st::Tuple{Vararg{Number,L´}}, ::Tuple{Vararg{Number,L}}) where {L,L´} = L´ == L - 1 ? st : + throw(ArgumentError("Malformed `subticks`. The number of subticks for each axis should be one less than the number of ticks for that axis")) + +function _cuboid(subticks, axesticks::Vararg{Tuple,L}) where {L} + allticks = ntuple(Val(L)) do i + allaxisticks = typeof(1.0)[] # We want the machine's float type, without committing to Float64 + axisticks = axesticks[i] + nticks = length(axisticks) + foreach(1:nticks-1) do j + append!(allaxisticks, range(axisticks[j], axisticks[j+1], length = subticks[i][j])) + j == nticks-1 || pop!(allaxisticks) end - finalizecolumn!(s) + allaxisticks end - adjmat = sparse(s) - - return Mesh(verts, adjmat) + return CuboidMesh(allticks) end -mesh(; kw...) = throw(ArgumentError("Need a finite number of ranges to define a marching mesh")) - -nvertices(m::Mesh) = length(m.vertices) - -nedges(m::Mesh) = div(nnz(m.adjmat), 2) - -nsimplices(m::Mesh) = length(simplices(m)) - -vertices(m::Mesh) = m.vertices - -edges(m::Mesh, src) = nzrange(m.adjmat, src) - -edgedest(m::Mesh, edge) = rowvals(m.adjmat)[edge] - -edgevertices(m::Mesh) = - ((vsrc, m.vertices[edgedest(m, edge)]) for (i, vsrc) in enumerate(m.vertices) for edge in edges(m, i)) - -function minmax_edge(m::Mesh{D,T}) where {D,T<:Real} - minlen2 = typemax(T) - maxlen2 = zero(T) - verts = vertices(m) - minedge = zero(first(verts)) - maxedge = zero(first(verts)) - for src in eachindex(verts), edge in edges(m, src) - dest = edgedest(m, edge) - dest > src || continue # Need only directed graph - vec = verts[dest] - verts[src] - norm2 = vec' * vec - norm2 < minlen2 && (minlen2 = norm2; minedge = vec) - norm2 > maxlen2 && (maxlen2 = norm2; maxedge = vec) - end - return minedge, maxedge -end - -transform!(f::Function, m::Mesh) = (map!(f, vertices(m), vertices(m)); m) - -###################################################################### -# Compute N-simplices (N = number of vertices) -###################################################################### -function simplices(mesh::Mesh{D}, ::Val{N} = Val(D+1)) where {D,N} - N > 0 || throw(ArgumentError("Need a positive number of simplex vertices")) - N == 1 && return Tuple.(1:nvertices(mesh)) - simps = NTuple{N,Int}[] - if nvertices(mesh) >= N - buffer = (NTuple{N,Int}[], NTuple{N,Int}[], Int[]) - for src in eachindex(vertices(mesh)) - append!(simps, _simplices(buffer, mesh, src)) - end - N > 2 && alignnormals!(simps, vertices(mesh)) - end - return simps -end - -# Add (greater) neighbors to last vertex of partials that are also neighbors of scr, till N -function _simplices(buffer::Tuple{P,P,V}, mesh, src) where {N,P<:AbstractArray{<:NTuple{N}},V} - partials, partials´, srcneighs = buffer - resize!(srcneighs, 0) - resize!(partials, 0) - for edge in edges(mesh, src) - srcneigh = edgedest(mesh, edge) - srcneigh > src || continue # Directed graph, to avoid simplex duplicates - push!(srcneighs, srcneigh) - push!(partials, padright((src, srcneigh), 0, Val(N))) - end - for pass in 3:N - resize!(partials´, 0) - for partial in partials - nextsrc = partial[pass - 1] - for edge in edges(mesh, nextsrc) - dest = edgedest(mesh, edge) - dest > nextsrc || continue # If not directed, no need to check - dest in srcneighs && push!(partials´, tuplesplice(partial, pass, dest)) - end - end - partials, partials´ = partials´, partials - end - return partials -end - -function alignnormals!(simplices, vertices) - for (i, s) in enumerate(simplices) - volume = elementvolume(vertices, s) - volume < 0 && (simplices[i] = switchlast(s)) - end - return simplices -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) - -switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N)) - -####################################################################### -# piecewise and isometric -####################################################################### -# a 1D mesh from 0 to N-1 with `points[i+1]` points in each segment -piecewise_mesh(nodes::NTuple{N,Any}, points::Int) where {N} = piecewise_mesh(nodes, filltuple(points, Val(N-1))) - -function piecewise_mesh(nodes::NTuple{N,Any}, points) where {N} - vsegments = ntuple(Val(N-1)) do i - v = SVector.(range(i-1, i, length = points[i])) - i == N-1 || pop!(v) - return v - end - verts = vcat(vsegments...) - nv = length(verts) - adjmat = sparse(vcat(1:nv-1, 2:nv), vcat(2:nv, 1:nv-1), true, nv, nv) - return Mesh(verts, adjmat) -end - -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 - -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 - -parsenode(pt::SVector, ::Val{N}) where {N} = padright(pt, Val(N)) -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 - -# ####################################################################### -# # LinearMeshSpec -# ####################################################################### -# struct LinearMeshSpec{N,L,T,R} <: MeshSpec{1} -# vertices::SVector{N,SVector{L,T}} -# samelength::Bool -# closed::Bool -# points::R -# end - -# """ -# linearmesh(nodes...; points = 13, samelength = false, closed = false) - -# Create a `MeshSpec` for a one-dimensional `Mesh` connecting the `nodes` with straight -# segments, each containing a number `points` of points (endpoints included). If a different -# number of points for each of the `N` segments is required, use `points::NTuple{N,Int}`. -# If `samelength` each segment has equal length in mesh coordinates. If `closed` the last node -# is connected to the first node (must be equal) - -# # Examples - -# ```jldoctest -# julia> buildmesh(linearmesh(:Γ, :K, :M, :Γ; points = (101, 30, 30)), HamiltonianPresets.graphene()) -# Mesh{1}: mesh of a 1-dimensional manifold -# Vertices : 159 -# Edges : 158 -# ``` - -# # See also -# `marchingmesh`, `buildmesh` -# """ -# linearmesh(nodes...; points = 13, samelength::Bool = false, closed::Bool = false) = -# LinearMeshSpec(sanitize_BZpts(nodes, closed), samelength, closed, points) +Base.eachindex(mesh::CuboidMesh) = CartesianIndices(eachindex.(mesh.ticks)) -# function sanitize_BZpts(pts, closed) -# pts´ = parse_BZpoint.(pts) -# if closed -# all(isapprox.(first(pts´), last(pts´))) || -# throw(ArgumentError("Closed linear meshes should have equal first and last nodes.")) -# end -# dim = maximum(length.(pts´)) -# pts´´ = SVector(padright.(pts´, Val(dim))) -# return pts´´ -# end +Base.size(mesh::CuboidMesh, i...) = size(eachindex(mesh), i...) -# parse_BZpoint(p::SVector) = float.(p) -# parse_BZpoint(p::Tuple) = SVector(float.(p)) -# function parse_BZpoint(p::Symbol) -# pt = get(BZpoints, p, missing) -# pt === missing && throw(ArgumentError("Unknown Brillouin zone point $p, use one of $(keys(BZpoints))")) -# return SVector(float.(pt)) -# end -# const BZpoints = -# ( Γ = (0,) -# , X = (pi,) -# , Y = (0, pi) -# , Z = (0, 0, pi) -# , K = (2pi/3, -2pi/3) -# , Kp = (4pi/3, 2pi/3) -# , M = (pi, 0) -# ) +Base.length(mesh::CuboidMesh) = prod(length.(mesh.ticks)) -# linearmesh_nodes(l, br) = cumsum(SVector(0, segment_lengths(l, br)...)) +vertex(mesh::CuboidMesh, n::CartesianIndex) = SVector(getindex.(mesh.ticks, Tuple(n))) -# function segment_lengths(s::LinearMeshSpec{N,LS,TS}, br::SMatrix{E,LB,TB}) where {TS,TB,N,E,LS,LB} -# T = promote_type(TS, TB) -# verts = padright.(s.vertices, Val(LB)) -# dϕs = ntuple(i -> verts[i + 1] - verts[i], Val(N-1)) -# if s.samelength -# ls = filltuple(T(1/(N-1)), Val(N-1)) -# else -# ibr = pinverse(br)' -# ls = (dϕ -> norm(ibr * dϕ)).(dϕs) -# ls = ls ./ sum(ls) -# end -# return ls -# end +vertices(mesh::CuboidMesh) = (SVector(v) for v in Iterators.product(mesh.ticks...)) -# function idx_to_node(s, br) -# nodes = SVector.(linearmesh_nodes(s, br)) -# nmax = length(nodes) -# nodefunc = nvec -> begin -# n = only(nvec) -# node = if n >= nmax -# nodes[nmax] -# else -# nc = max(n, 1) -# i = Int(floor(nc)) -# nodes[i] + rem(nc, 1) * (nodes[i+1] - nodes[i]) -# end -# return node -# end -# return nodefunc -# end +nvertices(mesh::CuboidMesh) = length(vertices(mesh)) -# function buildmesh(s::LinearMeshSpec{N}, br::SMatrix) where {N} -# ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), s.points) -# verts = SVector.(first(ranges)) -# for r in Base.tail(ranges) -# pop!(verts) -# append!(verts, SVector.(r)) -# end -# s.closed && pop!(verts) -# nv = length(verts) -# nodefunc = idx_to_node(s, br) -# verts .= nodefunc.(verts) -# adjmat = sparse(vcat(1:nv-1, 2:nv), vcat(2:nv, 1:nv-1), true, nv, nv) -# s.closed && (adjmat[end, 1] = adjmat[1, end] = true) -# return Mesh(verts, adjmat) -# end +neighbors(mesh::CuboidMesh, n::CartesianIndex) = marchingneighbors(eachindex(mesh), n) -# function buildlift(s::LinearMeshSpec, br::SMatrix{E,L}) where {E,L} -# ls = segment_lengths(s, br) -# nodes = linearmesh_nodes(s, br) -# verts = padright.(s.vertices, Val(L)) -# l = sum(ls) -# liftfunc = x -> begin -# xc = clamp(only(x), 0, l) -# for (i, node) in enumerate(nodes) -# if node > xc -# p = verts[i-1] + (xc - nodes[i-1])/ls[i-1] * (verts[i]-verts[i-1]) -# return p -# end -# end -# return last(verts) -# end -# return liftfunc -# end \ No newline at end of file +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 diff --git a/src/parametric.jl b/src/parametric.jl index a72df590..23fee886 100644 --- a/src/parametric.jl +++ b/src/parametric.jl @@ -201,9 +201,6 @@ Base.size(ph::ParametricHamiltonian, n...) = size(ph.h, n...) Base.eltype(ph::ParametricHamiltonian) = eltype(ph.h) -DualNumbers.Dual(p::ParametricHamiltonian) = - ParametricHamiltonian(Dual(p.baseh), Dual(p.h), p.modifiers, p.ptrdata, p.allptrs, p.parameters) - ####################################################################### # bloch! for parametric ####################################################################### diff --git a/src/plot_makie.jl b/src/plot_makie.jl index 83d27248..fd041007 100644 --- a/src/plot_makie.jl +++ b/src/plot_makie.jl @@ -247,14 +247,14 @@ end # plot(::Bandstructure) ####################################################################### -function plot(bs::Bandstructure{2}; kw...) +function plot(bs::Bandstructure{1}; kw...) scene = bandplot2d(bs; kw...) axis = scene[Axis] axis[:names, :axisnames] = ("φ", "ε") return scene end -function plot(bs::Bandstructure{3}; kw...) +function plot(bs::Bandstructure{2}; kw...) scene = bandplot3d(bs; kw...) axis = scene[Axis] to_value(scene[:show_axis]) && (axis[:names, :axisnames] = ("φ₁", "φ₂", "ε")) @@ -263,7 +263,7 @@ end @recipe(BandPlot2D, bandstructure) do scene Theme( - linewidth = 3, + linethickness = 3.0, wireframe = true, colors = map(t -> RGBAf0((0.8 .* t)...), ((0.973, 0.565, 0.576), (0.682, 0.838, 0.922), (0.742, 0.91, 0.734), @@ -278,18 +278,19 @@ function plot!(plot::BandPlot2D) colors = Iterators.cycle(plot[:colors][]) for (nb, color) in zip(bands, colors) band = bs.bands[nb] - vertices = band.mesh.vertices - simplices = band.simplices + vertices = band.verts + simplices = band.simpinds linesegments!(plot, (t -> vertices[first(t)] => vertices[last(t)]).(simplices), - linewidth = plot[:linewidth], color = color) + linewidth = plot[:linethickness][], color = color) end return plot end @recipe(BandPlot3D, bandstructure) do scene Theme( - linewidth = 1, - wireframe = false, + linethickness = 1.0, + wireframe = true, + linedarken = 0.5, ssao = true, ambient = Vec3f0(0.55), diffuse = Vec3f0(0.4), colors = map(t -> RGBAf0(t...), ((0.973, 0.565, 0.576), (0.682, 0.838, 0.922), (0.742, 0.91, 0.734), @@ -300,20 +301,20 @@ end function plot!(plot::BandPlot3D) bs = to_value(plot[1]) - bands = haskey(plot, :bands) ? to_value(plot[:bands]) : eachindex(bs.bands) + bandinds = haskey(plot, :bands) ? to_value(plot[:bands]) : eachindex(bs.bands) colors = Iterators.cycle(plot[:colors][]) - for (nb, color) in zip(bands, colors) + for (nb, color) in zip(bandinds, colors) band = bs.bands[nb] - vertices = band.mesh.vertices - connectivity = [s[j] for s in band.simplices, j in 1:3] + vertices = band.verts + connectivity = [s[j] for s in band.simpinds, j in 1:3] if isempty(connectivity) scatter!(plot, vertices, color = color) else mesh!(plot, vertices, connectivity, color = color, transparency = false, - ssao = plot[:ssao][], ambient = plot[:ambient][], diffuse = plot[:diffuse][]) + ssao = plot[:ssao][], ambient = plot[:ambient][], diffuse = plot[:diffuse][]) if plot[:wireframe][] - edgevertices = collect(Quantica.edgevertices(band.mesh)) - linesegments!(plot, edgevertices, linewidth = plot[:linewidth]) + edgevertices = collect(Quantica.edgevertices(band)) + linesegments!(plot, edgevertices, color = darken(color, plot[:linedarken][]), linewidth = plot[:linethickness][]) end end end diff --git a/src/plot_vegalite.jl b/src/plot_vegalite.jl index 3bfe38bb..fe3be7b2 100644 --- a/src/plot_vegalite.jl +++ b/src/plot_vegalite.jl @@ -94,7 +94,7 @@ link_shader(shader::Missing, psi, h) = (row, col) -> 0.0 # vlplot ####################################################################### """ - vlplot(b::Bandstructure{2}; kw...) + vlplot(b::Bandstructure{1}; kw...) Plot the 1D bandstructure `b` using VegaLite in 2D. @@ -120,6 +120,7 @@ shaders: `sitesize`, `siteopacity`, `sitecolor`, `linksize`, `linkopacity` or `l ## Bandstructure-specifc - `points = false`: whether to plot points on line plots - `bands = missing`: bands to plot (or all if `missing`) + - `thickness = 2`: thickness of band lines ## Hamiltonian-specific - `axes = (1,2)`:lattice axes to project onto the plot x-y plane @@ -152,7 +153,7 @@ maximum visual value of site diameter and link thickness is given by `maxdiamete pixels) and `maxthickness` (as a fraction of `maxdiameter`) """ function VegaLite.vlplot(b::Bandstructure; - labels = ("φ", "ε"), scaling = (1, 1), size = 640, points = false, xlims = missing, ylims = missing, bands = missing) + labels = ("φ", "ε"), scaling = (1, 1), size = 640, points = false, xlims = missing, ylims = missing, bands = missing, thickness = 2) labelx, labely = labels table = bandtable(b, make_it_two(scaling), bands) sizes = make_it_two(size) @@ -160,7 +161,7 @@ function VegaLite.vlplot(b::Bandstructure; plotrange = (xlims, ylims) (domainx, domainy), _ = domain_size(corners, size, plotrange) p = table |> vltheme(sizes, points) + @vlplot( - mark = :line, + mark = {:line, strokeWidth = thickness}, x = {:x, scale = {domain = domainx, nice = false}, title = labelx, sort = nothing}, y = {:y, scale = {domain = domainy, nice = false}, title = labely, sort = nothing}, color = "band:n", @@ -168,11 +169,27 @@ function VegaLite.vlplot(b::Bandstructure; return p end -function bandtable(b::Bandstructure{2}, (scalingx, scalingy), bandsiter) +# We optimize 1D band plots by collecting connected simplices into line segments (because :rule is considerabley slower than :line) +function bandtable(b::Bandstructure{1,T}, (scalingx, scalingy), bandsiter) where {T} bandsiter´ = bandsiter === missing ? eachindex(bands(b)) : bandsiter - bnds = bands(b) - table = [(x = v[1] * scalingx, y = v[2] * scalingy, band = i, tooltip = string(v)) - for i in bandsiter´ for v in vertices(bnds[i])] + NT = typeof((;x = zero(T), y = zero(T), band = 1, segment = 1)) + table = NT[] + for (nb, band) in enumerate(bands(b)) + segment = 0 + verts = vertices(band) + sinds = band.simpinds + s0 = first(sinds) + for s in sinds + if first(s) == last(s0) + push!(table, (; x = verts[first(s)][1] * scalingx, y = verts[first(s)][2] * scalingy, band = nb, segment = segment)) + else + push!(table, (; x = verts[last(s0)][1] * scalingx, y = verts[last(s0)][2] * scalingy, band = nb, segment = segment)) + segment += 1 + s0 = s + push!(table, (; x = verts[first(s)][1] * scalingx, y = verts[first(s)][2] * scalingy, band = nb, segment = segment)) + end + end + end return table end @@ -184,7 +201,7 @@ function VegaLite.vlplot(h::Hamiltonian{LA}, psi = missing; linksize = maxthickness, linkopacity = 1.0, linkcolor = sitecolor, colorrange = missing, colorscheme = "redyellowblue", discretecolorscheme = "category10") where {E,LA<:Lattice{E}} - psi´ = maybe_unflatten_or_missing(psi, h) + psi´ = unflatten_or_reinterpret_or_missing(psi, h) checkdims_psi(h, psi´) directives = (; axes = axes, digits = digits, @@ -272,8 +289,8 @@ sanitize_colorrange(r::Tuple, table) = [first(r), last(r)] needslegend(x::Number) = nothing needslegend(x) = true -maybe_unflatten_or_missing(psi::Missing, h) = missing -maybe_unflatten_or_missing(psi, h) = maybe_unflatten(psi, h) +unflatten_or_reinterpret_or_missing(psi::Missing, h) = missing +unflatten_or_reinterpret_or_missing(psi, h) = unflatten_or_reinterpret(psi, h) checkdims_psi(h, psi) = size(h, 2) == size(psi, 1) || throw(ArgumentError("The eigenstate length $(size(psi,1)) must match the Hamiltonian dimension $(size(h, 2))")) checkdims_psi(h, ::Missing) = nothing @@ -341,6 +358,10 @@ function _corners(table) for row in table min´ = min.(min´, (row.x, row.y)) max´ = max.(max´, (row.x, row.y)) + if isdefined(row, :x2) + min´ = min.(min´, (row.x2, row.y2)) + max´ = max.(max´, (row.x2, row.y2)) + end end return min´, max´ end diff --git a/src/tools.jl b/src/tools.jl index 7b33ce57..7daee2a9 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -86,7 +86,7 @@ end padright(sv::StaticVector{E,T}, x::T, ::Val{E}) where {E,T} = sv padright(sv::StaticVector{E1,T1}, x::T2, ::Val{E2}) where {E1,T1,E2,T2} = - (T = promote_type(T1,T2); SVector{E2,T}(ntuple(i -> i > E1 ? x : convert(T, sv[i]), Val(E2)))) + (T = promote_type(T1,T2); SVector{E2,T}(padright(Tuple(sv), x, Val(E2)))) padright(sv::StaticVector{E,T}, ::Val{E2}) where {E,T,E2} = padright(sv, zero(T), Val(E2)) padright(sv::StaticVector{E,T}, ::Val{E}) where {E,T} = sv padright(t::NTuple{N´,Any}, x, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? x : t[i], Val(N)) @@ -200,32 +200,32 @@ 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 +# 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} @@ -264,10 +264,6 @@ function append_slice!(dest::AbstractArray, src::AbstractArray{T,N}, Rsrc::Carte return dest end -dualarray(a::DenseMatrix) = map(x->Dual.(x, 0), a) -# Need to preserve stored zeros, so we have to treat sparse case as special -dualarray(s::SparseMatrixCSC) = SparseMatrixCSC(s.m, s.n, s.colptr, s.rowval, map(x->Dual.(x, 0), s.nzval)) - ###################################################################### # convert a matrix/number block to a matrix/inlinematrix string ###################################################################### diff --git a/test/test_bandstructure.jl b/test/test_bandstructure.jl index 3383b646..1a17033f 100644 --- a/test/test_bandstructure.jl +++ b/test/test_bandstructure.jl @@ -1,51 +1,59 @@ +using Quantica: nbands, nvertices, nedges, nsimplices + @testset "basic bandstructures" begin h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) - b = bandstructure(h, points = 13) - @test length(bands(b)) == 1 + b = bandstructure(h, subticks = 13) + @test nbands(b) == 1 + + h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(2) + b = bandstructure(h, cuboid(0.99 .* (-π, π), 0.99 .* (-π, π), subticks = 13)) + @test nbands(b) == 3 h = LatticePresets.honeycomb() |> hamiltonian(onsite(0.5, sublats = :A) + onsite(-0.5, sublats = :B) + hopping(-1, range = 1/√3)) - b = bandstructure(h, points = (13, 15)) - @test length(bands(b)) == 2 + b = bandstructure(h, subticks = (13, 15)) + @test nbands(b) == 2 h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2) - b = bandstructure(h, points = (5, 7, 5)) - @test length(bands(b)) == 8 + b = bandstructure(h, subticks = (5, 7, 5)) + @test nbands(b) == 1 - b = bandstructure(h, :Γ, :X; points = 4) - @test length(bands(b)) == 8 + b = bandstructure(h, :Γ, :X; subticks = 4) + @test nbands(b) == 1 - b = bandstructure(h, :Γ, :X, (0, π), :Z, :Γ; points = 4) - @test length(bands(b)) == 8 + b = bandstructure(h, :Γ, :X, (0, π), :Z, :Γ; subticks = 4) + @test nbands(b) == 1 + @test nvertices(b) == 73 - b = bandstructure(h, :Γ, :X, (0, π), :Z, :Γ; points = (4,5,6,7)) - @test length(bands(b)) == 8 + b = bandstructure(h, :Γ, :X, (0, π), :Z, :Γ; subticks = (4,5,6,7)) + @test nbands(b) == 1 + @test nvertices(b) == 113 end @testset "functional bandstructures" begin hc = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, sublats = :A=>:B, plusadjoint = true)) matrix = similarmatrix(hc, LinearAlgebraPackage()) hf((x,)) = bloch!(matrix, hc, (x, -x)) - m = mesh((0, 1)) + m = cuboid((0, 1)) b = bandstructure(hf, m) - @test length(bands(b)) == 2 + @test nbands(b) == 2 hc2 = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) hp2 = parametric(hc2, @hopping!((t; s) -> s*t)) matrix2 = similarmatrix(hc2, LinearAlgebraPackage()) hf2((s, x)) = bloch!(matrix2, hp2(s = s), (x, x)) - m2 = mesh((0, 1), (0, 1)) + m2 = cuboid((0, 1), (0, 1)) b = bandstructure(hf2, m2) - @test length(bands(b)) == 2 + @test nbands(b) == 1 end @testset "bandstructures lifts & transforms" begin h = LatticePresets.honeycomb() |> hamiltonian(onsite(2) + hopping(-1, range = 1/√3)) - mesh1D = mesh((0, 2π)) + mesh1D = cuboid((0, 2π)) b = bandstructure(h, mesh1D, mapping = φ -> (φ, -φ), transform = inv) b´ = transform!(inv, bandstructure(h, mesh1D, mapping = φ -> (φ, -φ))) - @test length(bands(b)) == length(bands(b´)) == 2 + @test nbands(b) == nbands(b´) == 1 @test vertices(bands(b)[1]) == vertices(bands(b´)[1]) h´ = unitcell(h) s1 = spectrum(h´, transform = inv) @@ -53,35 +61,35 @@ end @test energies(s1) == energies(s2) # no automatic mapping from 2D to 3D h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2) - @test_throws DimensionMismatch bandstructure(h, mesh((0, 2pi), (0, 2pi))) + @test_throws DimensionMismatch bandstructure(h, cuboid((0, 2pi), (0, 2pi))) end @testset "parametric bandstructures" begin ph = LatticePresets.linear() |> hamiltonian(onsite(0I) + hopping(-I), orbitals = Val(2)) |> unitcell(2) |> parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k = 2, p = [1,2])-> t - k*I + p'p)) - mesh2D = mesh((0, 1), (0, 2π), points = 15) + mesh2D = cuboid((0, 1), (0, 2π), subticks = 15) b = bandstructure(ph, mesh2D, mapping = (x, k) -> (x, (;k = k))) - @test length(bands(b)) == 4 + @test nbands(b) == 4 b = bandstructure(ph, mesh2D, mapping = (x, k) -> ((x,), (;k = k))) - @test length(bands(b)) == 4 + @test nbands(b) == 4 b = bandstructure(ph, mesh2D, mapping = (x, k) -> (SA[x], (;k = k))) - @test length(bands(b)) == 4 + @test nbands(b) == 4 b = bandstructure(ph, mesh2D, mapping = (k, φ) -> (1, (;k = k, p = SA[1, φ]))) - @test length(bands(b)) == 4 + @test nbands(b) == 1 ph = LatticePresets.linear() |> hamiltonian(onsite(0I) + hopping(-I), orbitals = Val(2)) |> unitcell(2) |> unitcell |> parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k = 2, p = [1,2])-> t - k*I + p'p)) b = bandstructure(ph, mesh2D, mapping = (k, φ) -> (;k = k, p = SA[1, φ])) - @test length(bands(b)) == 4 + @test nbands(b) == 1 b = bandstructure(ph, mesh2D, mapping = (k, φ) -> ((;k = k, p = SA[1, φ]),)) - @test length(bands(b)) == 4 + @test nbands(b) == 1 @test_throws UndefKeywordError bandstructure(ph, mesh2D, mapping = (k, φ) -> ((;p = SA[1, φ]),)) end @testset "unflatten" begin h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = (Val(2), Val(1))) |> unitcell(2) |> unitcell sp = states(spectrum(h))[:,1] - sp´ = Quantica.maybe_unflatten(sp, h) + sp´ = Quantica.unflatten_or_reinterpret(sp, h) l = size(h, 1) @test length(sp) == 1.5 * l @test length(sp´) == l @@ -91,7 +99,7 @@ end h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = Val(2)) |> unitcell(2) |> unitcell sp = states(spectrum(h))[:,1] - sp´ = Quantica.maybe_unflatten(sp, h) + sp´ = Quantica.unflatten_or_reinterpret(sp, h) l = size(h, 1) @test length(sp) == 2 * l @test length(sp´) == l @@ -99,6 +107,6 @@ end h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = Val(2)) |> unitcell(2) |> unitcell sp = states(spectrum(h))[:,1] - sp´ = Quantica.maybe_unflatten(sp, h) + sp´ = Quantica.unflatten_or_reinterpret(sp, h) @test sp === sp end \ No newline at end of file diff --git a/test/test_mesh.jl b/test/test_mesh.jl index c305b4c9..a196894f 100644 --- a/test/test_mesh.jl +++ b/test/test_mesh.jl @@ -1,7 +1,7 @@ using Quantica: nvertices, nedges, nsimplices @testset "meshes" begin - m = mesh((0,1), (0,2), points = (10, 20)) - @test nvertices(m) == 200 && nedges(m) == 541 && nsimplices(m) == 342 - @test_throws MethodError mesh((SA[1,2], SA[2,3])) + m = cuboid((0,1), (0,2), subticks = (10, 20)) + @test nvertices(m) == 200 + @test_throws MethodError cuboid((SA[1,2], SA[2,3])) end \ No newline at end of file