Redesign of band cuts

[BacAmorim]

I was trying to implement a way to plot a bandstructure along a 1D path in the Brillouin zone (e.g.: plot the band structure of graphene along the path Gamma -> K -> M -> Gamma as a function of the arclength of the path). I managed to implement a function linemesh(verts; npts = 20, closed = false) that given vertices verts generates a piecewise linear mesh that connects the vertices. I also defined a function arclength(linepath) that returns a vector with the arclength measured along the 1D path.

So to compute the band structure along the path, I would like to do something like:

kpath = linemesh([Gamma, K, M], npts = 100, closed = true) # closed = true, because I want the closed path Gamma -> K -> M -> Gamma

bs = bandstructure(graphene, kpath)

Now the problem comes when I want to plot the band structure. Ideally I would like to do something like

using PyPlot # sorry, it is what I use

arc = arclength(kpath)
plot(arc, [x[end] for x in vertices(bs, 1)]) # plot the 1st band
plot(arc, [x[end] for x in vertices(bs, 2)]) # plot the 2nd band

The problem lies in the fact that in vertices(bs, i) the k points are reordered in a way that is unrelated to the original order in kpath (as far as I can tell).

How could we work arround this?

A possible solution would the to modify bandstrucure such that the cut argument also accepts arrays: such that it would be possible to make:

bs = bandstructure(graphene, arc, cut = kpath)

where kpath and cut must have the same length. Therefore, the vertices in the band structure would correspond to he 1D arc mesh and plotting should work with no problems by simply doing:

plot([x[1] for x in vertices(bs, 1)], [x[end] for x in vertices(bs, 1)]) # plot the 1st band
plot([x[1] for x in vertices(bs, 1)], [x[end] for x in vertices(bs, 2)]) # plot the 2nd band

since [x[1] for x in vertices(bs, 1)] contains the arclength points.

Would this be possible and/or desirable?

[BacAmorim]

After digging a bit in the bandstructure code, I belive my proposal would require a rewriting of the internal function

_bandstructure(matrixf::Function, matrix´::AbstractMatrix{M}, mesh::MD, d::DiagonalizeHelper)

to something like

_bandstructure(matrixf::Function, matrix´::AbstractMatrix{M}, mesh::MD, comesh::Mesh, d::DiagonalizeHelper)

The previous functionally with cut would be obtained by building

comesh = Mesh([cut(v) for v in vertices(mesh)], mesh.adjmat)

[pablosanjose]

Hi @BacAmorim, thanks for the issue. You might want to check #28 and #52, which are a followup of the discussion we had on #17. Check out the help docstring for bandstructure too. We can now perform cuts like the ones you need in a bandstructure call, using the bandstructure(h, mesh; cut= (function)) syntax. The key difference is that mesh is a discretization of the cut itself, and cut just maps coordinates from said mesh to the overall Brillouin zone + parameter space. Do let me know if this solves your issue, or if any extension is needed, please.

A future extension I had planned is to accept a piecewise 1D cut specification in terms of points in the Brillouin zone, so that one could do things like bandstructure(h, mesh, cut = (point1, point2, point3...)), which would then build a cut function using appropriate arclengths

[pablosanjose]

Ah, my apologies, I had misunderstood your issue, you are of course using the cut functionality, but have problems with the output ordering.

[BacAmorim]

My proposal would be to modify the function _bandstructure to something that takes a mesh and a comesh:

function Quantica._bandstructure(matrixf::Function, matrix´::AbstractMatrix{M}, mesh::MD, comesh::MDco, d::Quantica.DiagonalizeHelper) where {M,D,T,MD<:Quantica.Mesh{D,T}, MDco<:Quantica.Mesh}
    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

    dimh = size(matrix´, 1)
    nk = Quantica.nvertices(mesh)
    # function to apply to eigenvalues when building bands (depends on momenta type)
    by =  Quantica._maybereal(T)

    p =  Progress(nk, "Step 1/2 - Diagonalising: ")
    for (n, ϕs) in enumerate(vertices(comesh))
        matrix = matrixf(Tuple(ϕs))
        # (ϵk, ψk) = diagonalize(Hermitian(matrix), d)  ## This is faster (!)
        (ϵk, ψk) =  Quantica.diagonalize(matrix, d.method)
         Quantica.resolve_degeneracies!(ϵk, ψk, ϕs, matrix, 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, dimh, nϵ, nk)
        end
         Quantica.copyslice!(ϵks, CartesianIndices((1:nϵ, n:n)),
                   ϵk,  CartesianIndices((1:nϵ,)), by)
         Quantica.copyslice!(ψks, CartesianIndices((1:dimh, 1:nϵ, n:n)),
                   ψk,  CartesianIndices((1:dimh, 1:nϵ)), identity)
         ProgressMeter.next!(p; showvalues = ())
    end

    p =  Progress(nϵ * nk, "Step 2/2 - Connecting bands: ")
    pcounter = 0
    bands = Quantica.Band{M,Vector{M},Quantica.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 = CartesianIndex{2}[]
    sizehint!(pending, nk)
    while true
        src = findfirst(iszero, vertindices)
        src === nothing && break
        resize!(pending, 1)
        pending[1] = src # source CartesianIndex for band search
        band =  Quantica.extractband(mesh, pending, ϵks, ψks, vertindices, d.minprojection)
        nverts =  Quantica.nvertices(band.mesh)
        nverts > D && push!(bands, band) # avoid bands with no simplices
        pcounter += nverts
        ProgressMeter.update!(p, pcounter; showvalues = ())
    end
    return Quantica.Bandstructure(bands, mesh)
end

We would have two methods for bandstructure: one that takes a mesh and a comesh:

function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh, comesh::Mesh;
                       method = defaultmethod(h), transform = missing, kw...)
    # ishermitian(h) || throw(ArgumentError("Hamiltonian must be hermitian"))
    length(mesh.vertices) == length(comesh.vertices) || throw(ArgumentError("Number of vertices in mesh and comesh must be the same"))
    matrix = similarmatrix(h, method)
    matrixf(φs) = bloch!(matrix, h, φs)
    codiag = codiagonalizer(matrixf, matrix, comesh)
    d = DiagonalizeHelper(method, codiag; kw...)
    b = _bandstructure(matrixf, matrix, mesh, comesh, d)
    transform === missing || transform!(transform, b)
    return b
end

And the current method, which takes a cut function

function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
                       method = defaultmethod(h), cut::Function, transform = missing, kw...)
    
    copts = [cut(kpt) for kpt in mesh.vertices]
    comesh = Quantica.Mesh(copts, mesh.adjmat)
    
    bandstructure(h, mesh,  comesh; method = method, transform = transform, kw...)
end

This would solve the problem for plotting of band structure along a path, by creating a mesh with k-points, kmesh and the corresponding 1D mesh of arclength measured along the path, arcmesh and calling:

bandstructure(h, arcmesh, kmesh)

Edit: typos fixed.

[BacAmorim]

Of course, the use case with no comesh/cut, would be implemented with the method:

function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
                       method = defaultmethod(h), transform = missing, kw...)
  
    bandstructure(h, mesh, mesh; method = method, transform = transform, kw...)
end

[pablosanjose]

This proposal seems to be intended to replace the current approach for cut, right? The idea is to have two meshes, one defined in L dimensions of the original BZ and another in the L´ < L dimensions of the cut manifold, yes?

I think that requiring to build two meshes is unnecesarily complicated. I did consider this when designing the current cut, but I decided in favor of requiring a single mesh (for the cut manifold) plus a mapping function to the original BZ.

If I understand correctly, the core problem you are having is simply the ordering of the vertices. Can you provide a self-contained example?

Wouldn't this approach work for your case?

julia> using Plots, Quantica

julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1,range=1/sqrt(3))) |> unitcell(3);

julia> mesh1D = marchingmesh(range(0, 3, length = 13));

julia> function cutfunc(x)
           kn = if x<1
               x * SA[1/3,1/3]
           elseif 1<=x<2
               SA[1/3,1/3] + (x-1)*(SA[1/2,0]-SA[1/3,1/3])
           else
               SA[1/2,0] + (x-2)*(SA[0,0]-SA[1/2,0])
           end
           return 2π * kn
       end
cutfunc (generic function with 1 method)

julia> bs = bandstructure(h, mesh1D; cut = cutfunc)
Bandstructure{1}: collection of 1D bands
  Bands        : 18
  Element type : scalar (Complex{Float64})
  Mesh{1}: mesh of a 1-dimensional manifold
    Vertices   : 13
    Edges      : 12

julia> plot([x[1] for x in vertices(bs, 1)], [x[end] for x in vertices(bs, 1)])

julia> plot!([x[1] for x in vertices(bs, 2)], [x[end] for x in vertices(bs, 2)])

If you need to have each segment have a different length, there are a number of ways to do this, either by having mesh1D have varying vertex densities or transorming the x[1] in the plot command. Ideally, what I would do is create a new method specialized for 1D cuts bandstructure(h; cut = (p1, p2, ...)) that builds the 1D mesh and the cutfunc for you from the points pi.

[pablosanjose]

By the way, you might be aware of this already, but if you don't like using Makie to plot the bands we now also have support for the beautiful VegaLite package. The interface is simple

julia> using VegaLite

julia> vlplot(bs)

[BacAmorim]

This proposal seems to be intended to replace the current approach for cut, right? The idea is to have two meshes, one defined in L dimensions of the original BZ and another in the L´ < L dimensions of the cut manifold, yes?

That is correct.

The problem with your approach is that we would have to define a cut function for each different path. I think that plotting a band structure along a 1D path is a frequent enough task that it deserves a prebuilt solution. As I said before, I also purpose to define a function linemesh and arclength:

function linemesh(pts; npts = 20, closed = false)
    
    verts = [p for p in pts]
    if closed
        push!(verts, pts[1])
    end
    
    nverts = length(verts)
    nedges = nverts - 1
    
    edgelengths = [norm(verts[n+1]-verts[n]) for n in 1:nedges]
    totallength = sum(edgelengths)
    
    nptsedge = [Int(div((npts-1)*edgelengths[n], totallength)) for n in 1:nedges]
    remedge = [npts*edgelengths[n]/totallength - nptsedge[n] for n in 1:nedges]
    perm = sortperm(remedge, rev = true)
    r = npts - 1 - sum(nptsedge)
    p = 1
    while r>0
        nptsedge[perm[p]] += 1
        r -= 1
        p += 1
    end
    
    kptsmesh = typeof(verts[1])[]
    
    for n in 1:nedges
        for m in 0:nptsedge[n]-1
            push!(kptsmesh, verts[n] + (verts[n+1]-verts[n])*m/nptsedge[n])
        end
    end
    
    push!(kptsmesh, verts[end])
    
    I = Int[]
    J = Int[]
    V = Bool[]
    for n = 1:npts-1
        push!(I, n)
        push!(J, n+1)
        push!(V, true)
        push!(I, n+1)
        push!(J, n)
        push!(V, true)
    end
    
    if closed
        push!(I, npts)
        push!(J, 1)
        push!(V, true)
        push!(I, 1)
        push!(J, npts)
        push!(V, true)
    end
    
    return Quantica.Mesh(kptsmesh,  sparse(I, J, V))
    
end

and

function arclength(mesh1D)
    
    accu = zero(norm(mesh1D.vertices[1]))
    arcvec = SVector{1, typeof(accu)}[]
    push!(arcvec, SVector(accu))
    
    for n in 2:length(mesh1D.vertices)
        
        accu += norm(mesh1D.vertices[n]-mesh1D.vertices[n-1])
        push!(arcvec, SVector(accu))
        
    end
    
    return Quantica.Mesh(arcvec, mesh1D.adjmat)
end   

With these functions and my modification to bandstruture, we could simply do:

a1 = SVector(0.5, sqrt(3)/2)
a2 = SVector(-0.5, sqrt(3)/2)

lat = lattice(
bravais(a1, a2), 
sublat(SVector(0.0, 0.0), name = :A), 
sublat(SVector(0.0, 1/sqrt(3)), name = :B)
)

graphene = hamiltonian(
lat, 
hopping(-1.0, sublats = :A => :B, range = 1) + hopping(-1.0, sublats = :B => :A, range = 1)
)


Γ = SVector(0.0, 0.0)
K = 2pi*SVector(1/3, -1/3)
M = 2pi*SVector(1/2, 0)

kptsmesh = linemesh((Γ, K, M), npts = 100, closed = true)
arcmesh = arclength(kptsmesh)

gbands = bandstructure(graphene, arcmesh, kptsmesh)

using PyPlot

scatter([x[1] for x in vertices(gbands, 1)], [x[end] for x in vertices(gbands, 1)])
scatter([x[1] for x in vertices(gbands, 2)], [x[end] for x in vertices(gbands, 2)])

[BacAmorim]

Ideally, what I would do is create a new method specialized for 1D cuts bandstructure(h; cut = (p1, p2, ...)) that builds the 1D mesh and the cutfunc for you from the points pi.

I think that would be a nice interface, yes!
But I am not sure that having a method that buils the cutfuntion would be a simpler way of doing things. I know that in julia functions are fist class, but it seems to me that creating a cut mesh instead of a cut function, would be a cleaner way of doing things. Besides, I already have the methods implemented, so we just need to settle on the interface.

Is there any advantage in _bandstructure accepting a function instead of my proprosal of it taking a mesh and a comesh?

[BacAmorim]

Some more info:

The band extraction algorithm goes a bit crazy when I consider a closed path: (same example as before):

kptsmesh = linemesh((Γ, K, M), npts = 100, closed = true)
arcmesh = arclength(kptsmesh)

gbands_closed = bandstructure(graphene, arcmesh, kptsmesh)

using PyPlot

plot([x[1] for x in vertices(gbands_closed, 1)], [x[end] for x in vertices(gbands_closed, 1)])
plot([x[1] for x in vertices(gbands_closed, 2)], [x[end] for x in vertices(gbands_closed, 2)])

With an open path it words fine:

kptsmesh = linemesh((Γ, K, M, Γ), npts = 100, closed = false)
arcmesh = arclength(kptsmesh)

gbands_open = bandstructure(graphene, arcmesh, kptsmesh)

using PyPlot

plot([x[1] for x in vertices(gbands_open, 1)], [x[end] for x in vertices(gbands_open, 1)])
plot([x[1] for x in vertices(gbands_open, 2)], [x[end] for x in vertices(gbands_open, 2)])

[pablosanjose]

The advantage of a function-based approach is merely for efficiency. If you have a large mesh and comesh, allocating both in memory is redundant, and could be costly. You actually only need memory for one of them. Basically the two meshes in your proposal are just used to encode a mapping (from cut coordinates to BZ coordinates), which is more natural to encode in a function (that's what the mapping is!).

Having said that it is completely true that adding support for 1D cuts in terms of points in the API is a good idea. What I don't see the advantage of is requiring the user to build two meshes. The information we need from the user is the sequence of points in the BZ, and perhaps the number of sampling points for each segment. Requiring anything beyond this from the user is unnecessary. This (bandstructure(h; cut = (pts...))) would be analogous to the bandstructure(h; resolution = 13) API, which you can use without bothering with meshes.

[pablosanjose]

The band extraction algorithm goes a bit crazy when I consider a closed path: (same example as before):

Hm, interesting, I think I know what is going on here. In general the band-extraction algorithm returns vertices in the order resulting from an adjacent neighbor search, starting from the first vertex in the mesh. If the mesh is compact (i.e. closed), the first vertex is connected to two neighbors, so the search goes both backwards and forwards. This is an unavoidable situation in dimensions greater than 1. But in 1D there is a natural order to the band vertices. In this case we can make sure we sort the vertices according to their lexicographic order, which is trivial to do. For the moment try

let vs = sort(vertices(gbands_closed, 1)); plot([x[1] for x in vs], [x[end] for x in vs]); end

[pablosanjose]

I moved the closed mesh issue to #67 to keep it separate from the comesh/cutfunc discussion

[BacAmorim]

The advantage of a function-based approach is merely for efficiency. If you have a large mesh and comesh, allocating both in memory is redundant, and could be costly. You actually only need memory for one of them. Basically the two meshes in your proposal are just used to encode a mapping (from cut coordinates to BZ coordinates), which is more natural to encode in a function (that's what the mapping is!).

Hmm. I was thinking that you might be worried about the size of the mesh and having to duplicate it. I am starting to think that things could be simplified if _banstructure was split in two functions: (i) the diagonalization part _band_diagonalize, and (ii) the band extraction part _band_extract. This would propably be more flexible. Then we can have two methods for _band_diagonalize one that takes a a mesh and a cut function and other that takes the a mesh and a comesh. Or would this create to much code duplication?

What I don't see the advantage of is requiring the user to build two meshes.

No need for that actually. We can just define:

function Quantica.bandstructure(h::Union{Quantica.Hamiltonian,Quantica.ParametricHamiltonian};
                       method = defaultmethod(h), cut::Tuple{NTuple{N, SV}, Int}, transform = missing, kw...) where {N, SV}
    
    kptsmesh = linemesh(cut[1], npts = cut[2], closed = false)
    arcmesh = arclength(kptsmesh)
    bandstructure(h, arcmesh, kptsmesh; method = method, transform = transform, kw...)
end

such that the user only writes: gbands_open = bandstructure(graphene, cut = ((Γ, K, M, Γ), 100))

[pablosanjose]

(i) the diagonalization part _band_diagonalize, and (ii) the band extraction part _band_extract.

That refactor would indeed make sense if we decided to use comeshes

No need for that actually. We can just define:

In reality I'm most worried about good API design and excessive number of options for the same thing (#19). I'd like to have a conceptually simple bandstructure API, just bandstructure(h, [mesh]; opts...). It is easy to explain what mesh is (a discretization of the manifold for the bands to compute). So if we want to do cuts through a mesh-comesh combination I would prefer to do

bandstructure(h, mesh; cut = mappedmesh)

Now, this is a possibility, indeed, but I think it is the wrong abstraction. The cut is the mapping, not the mapped mesh. If we allow to provide an mappedmesh as a cut, it should be because the freedom offered by mappedmesh makes sense. What if mappedmesh is incompatible in any way with mesh? That would have to throw an error. This suggests to me that the correct abstraction is not the mapped mesh, but the mapping itself. The performance/allocation issue is also there, but it is quite secondary. It is mostly about a conceptually sound abstraction for what a cut is.

So what is a cut, precisely? It is the function that, applied to the vertices in mesh, produces mappedmesh. That's what the current API exposes. The fact that the mesh is open or closed, or its number of sample points, belongs in mesh, not cut.

If we agree on that, the only issue that remains is how to make it easier for the user to build a 1D cut as a polygon through a set of points. The points themselves should be the input, and there should be no need to provide a mesh in this case, unless we want to set its resolution on each segment, its open/close property or whether each segment has an equal length or a length consistent with its original length in the BZ. If we want to support a mesh-free API for 1D cuts we would then need to be able to accept input for all these extra properties. This suggests a possible API like

bandstructure(h, mesh; cut::Function)
bandstructure(h; cut = Linecut(pts...; resolution, closed::Bool, samelength::Bool)

where Linecut is a type representing a polygon cut. We could even enrich this type to support common symbolic pts, such as :Γ, :K or :X, thus encapsulating all this functionality separately from the basic and general bandstructure(h, mesh; cut::Function) codebase.

I'll try to make a PR hashing out all this soon.

[BacAmorim]

I see your point. Your idea of using a Linecut might indeed solve the problem nicelly. I think this would work as you intend:

struct LineMesh{D, T}
    vertices::Vector{SVector{D, T}}
    nverts::Int
    npoints::Int
    edgelengths::Vector{T}
    intervals::Vector{T}
end

# the piecewise linear path is obtained from the Linecutby  s: [0, 1] -> line(s)
function (line::LineMesh)(s) 
    @assert 0 <= s<= 1
    arclength = s*line.intervals[end]
    n=2
    while arclength>line.intervals[n]
        n +=1
    end
    
    return line.vertices[n-1] + (line.vertices[n]-line.vertices[n-1])*(arclength - line.intervals[n-1])/(line.intervals[n]-line.intervals[n-1])
end

function linecut(verts, npts)
    nverts = length(verts)
    nedges = nverts - 1
    
    edgelengths = [norm(verts[n+1]-verts[n]) for n=1:nedges]
    intervals = cumsum(edgelengths)
    
    return LineMesh(verts, nverts, npts, edgelengths, pushfirst!(intervals, zero(intervals[1])))
end

Then doing something like

kline = linemesh([Γ, K, M, Γ], 100)
scatter([kline(s)[1] for s=0:0.01:1], [kline(s)[2] for s=0:0.01:1])

would give:

So maybe the general interface shoud be something like

bandstructure(h; mesh = MarchingMesh(..., cut = ...)) # previous case. The MarchingMesh accepts as an option a cut function to uplift the lower dimensional mesh, to the higher one
bandstructure(h; mesh = LineMesh(verts, npoints))
```?

Since we are on the topic of API: is cut a good name? It implies that we have a higher dimensional lattice and taking a slice of it. While what we are actually doing is a pushforward (or uplifting) from the lower dimensionality mesh to the higher one.

[pablosanjose]

So maybe the general interface shoud be something like

Mmm, this is interesting! So, bandstructure only really needs a Hamiltonian (or a function that produces a matrix over a parameter space). The mesh is optional, defaulting to a marchingmesh over the full Brillouin zone with a default resolution = 13.

I have two problems with this however

• In the case of bandstructure(::Function) we need the mesh as input, because we cannot build it from the function (we cannot know the dimensionality).
• I also worry about incorporating the cut (or let's call it lift, is that a good name?) inside the Mesh type. It is a bit strange to add such a field to a Mesh, it seems out of place. But there is also a non-philosophical issue. The problem is that for the LineMesh case, I need h to define the lift field, as it requires once more that we know the number of parameters or lattice dimensions. Hence, I cannot build the mesh from just LineMesh(verts, npoints)

The alternative is to have lift as an extra kwarg, separate from the mesh, that can be either a function (easy case), or a PolygonLine (better name than Linearcut?) for the 1D cut use case.

Please do give this some critical thought @BacAmorim, this kind of input is very valuable.

[BacAmorim]

So maybe we can define

abstract type MeshBuilder end

struct MarchingMesh <: MeshBuilder 
...
end

struct LinearMesh <: MeshBuilder 
...
end

Then we can use multiple dispatch for the interface:

bandstructure(h::Hamiltonian, mesh::MarchingMesh): does what we currently have
bandstructure(h::Hamiltonian, mesh::LinearMesh): solves the linear path problem

If we just call bandstructure(h::Hamiltonian) it defaults to bandstructure(h::Hamiltonian, MarchingMesh(resolution=13))

We also define: bandstructure(f::Function; mesh::MarchingMesh) for which we must pass a mesh argument, otherwise it fails.

lift would be a good name!

However, I believe it is better to encode the lift function in the MeshBuilder, rather than as an extra argument to bandstructure. I think it is a better organization of arguments: we have some arguments for the mesh builder, some for the eigensolver and so on.

In this way we don't even have to change the underlying _bandstructure function: It always takes a concrete mesh and a lift function. The particular lift is specified by the MeshBuilder, either by using functors in my proposal above or by defining: lift(::MeshBuilder, x), used in the user facing bandstructure. The user always has to provide a function/hamiltonian and a mesh builder.

For the particular case of a hamiltonian, if no mesh is provided it can use the default MarchingMesh with predetermied parameters. Or maybe we always force the user to provide a mesh and have a predefined MarchingMesh() maybe using Parameters.jl

What do you think?

[pablosanjose]

I'm a bit confused, let me see. So you propose the user never to input a Mesh object, just MeshBuilder objects, which would then produce a Mesh internally?

I'm currently working on the linecut branch, where I am experimenting with something similar, but that still keeps Mesh accessible to the user. I just add one new method for bandstructure and a Linecut builder

bandstructure(h::Hamiltonian, lc::Linecut; resolution, ...)
linecut(pts...; closed = false, samelength = false)

This Linecut is not a Mesh, and is not a lift function. It is a different type that takes care of the logic of writing things line linecut(:Γ, :K, :M, :Γ). Once we have a lc::Linecut, we combine lc and h to produce a 1D Mesh and an appropriate piecewise lift function, which are then passed to bandstructure(h, mesh; lift=...). This is somewhat similar to your MeshBuilder, if I understand correctly, but without compounding it with other types of Meshes (we could add more complexity in the Mesh type in the future, independently from this). So in my view the user can still create a specific mesh, transform it, refine it, whatever, and then pass it to bandstructure, quite independently from the Linecut stuff.

[BacAmorim]

I'm a bit confused, let me see. So you propose the user never to input a Mesh object, just MeshBuilder objects, which would then produce a Mesh internally?

Exactly! A general MeshBuilder would store information to be able to do two things:

• generate a mesh
• create a lift function if necessary (or this can be just the identity)

This Linecut is not a Mesh, and is not a lift function. It is a different type that takes care of the logic of writing things line linecut(:Γ, :K, :M, :Γ). Once we have a lc::Linecut, we combine lc and h to produce a 1D Mesh and an appropriate piecewise lift function, which are then passed to bandstructure(h, mesh; lift=...). This is somewhat similar to your MeshBuilder,

I agree. If I understand correctly, this Linecut is essentially what I was proposing to be a LineMesh , which can build a 1D mesh and provide a lift function (implemented as a functor here):

struct LineMesh{D, T}
    vertices::Vector{SVector{D, T}}
    nverts::Int
    npoints::Int
    edgelengths::Vector{T}
    intervals::Vector{T}
end

# the piecewise linear path is obtained from the Linecutby  s: [0, 1] -> line(s)
function (line::LineMesh)(s) 
    @assert 0 <= s<= 1
    arclength = s*line.intervals[end]
    n=2
    while arclength>line.intervals[n]
        n +=1
    end

function linecut(verts, npts)
    nverts = length(verts)
    nedges = nverts - 1
    
    edgelengths = [norm(verts[n+1]-verts[n]) for n=1:nedges]
    intervals = cumsum(edgelengths)
    
    return LineMesh(verts, nverts, npts, edgelengths, pushfirst!(intervals, zero(intervals[1])))
end

Now, in truth, it is not necessary to create an aditional MeshBuilder type. We could just dispatch on the LineMesh/Linecut type as you purpose. I just think it would be a more unified approach, where what is passed to bandstructure is always hamiltonian + mesh builder (which contains all the optitions regarding the mesh).

[pablosanjose]

I have grown to like your unifying proposal, and have drafted #70 along those lines. Please let me know what you think! I'm particularly interested to know if the unified API is sufficiently natural, or feels over-engineered.