Support Dirac points in bandstructures

[pablosanjose]

The isse #122 proposes to introduce the struct Simplices that contains a svecs::Vector{NTuple{D´,S}}, where S<:SubArray. So, vecs in each simplex are views of eigenvector matrices produced by diagonalization. This saves a lot of space and indirections, but it also opens the door to solving a long-standing issue in an elegant way: the correct formation of simplices with Dirac-point degeneracies at a vertex. As a bonus, we also solve the case of random non-Abelian connection in multiply degenerate bands.

The issue with Dirac points in a vertex is that, however you choose the basis of the degenerate eigenstate in such a vertex, you will find a frustrated connectivity of band vertices, as defined by eigenstates. The Dirac point is actually a dislocation in connectivity. The issue cannot be solved as long as eigenstates remain uniquely associated to vertices. In the svecs::Vector{NTuple{D´,S}} implementation we have the freedom, however, to have two different eigenstate bases for the same Dirac-point vertex within two adjacent simplices that share said vertex. Each simplex needs to keep an independent copy of the degenerate view and rediagonalize independently. Exploiting this freedom, the Dirac point problem is solved.

This idea is cool, but an actual implementation is far from trivial. It involves changing how we deal with degenerate eigenstates in the first place. We don't need to codiagonalize them as we diagonalize vertices. We do it on a simplex-by-simplex basis when extracting simplices.

This issue is a stab at designing such implementation. It involves dropping the whole codiagonalization code altogether, and has the potential of increasing bandstructure build performance. Since it is important to document the details for the future, we describe the algorithm is full depth here.

Revised Bandstructure algorithm

The redesign of Bandstructure would be, according to #122

struct BaseMesh{D,V} <: AbstractMesh{D}  # Was CuboidMesh, implements AbstractMesh API (edges, edgedest, vertices, simplices...)
    axesticks::NTuple{D,V}  # Here V could be any container, although tuples are fastest for indexing
end

struct BandMesh{D´,T<:Number} <: AbstractMesh{D´}  # Was Mesh in master, D´ = D+1 only
    verts::Vector{SVector{D´,T}}
    adjmat::SparseMatrixCSC{Bool,Int}
    simpinds::Vector{NTuple{D´,Int}}    # Needed for Makie plotting. Could be computed each time
end

struct Simplices{D´,T,S<:SubArray}
    sverts::Vector{NTuple{D´,SVector{D´,T}}}
    svecs::Vector{NTuple{D´,S}}
    sptrs::Array{UnitRange{Int},D´}  # range of indices of verts and vecs for each base simplex index
end

struct Bandstructure{D,T,M<:BaseMesh{D},D´,B<:BandMesh{D´,T},S<:Simplices{D´,T}}   # D is dimension of base space, D´ = D+1
    base::M           # We don't allow a generica AbstractMesh here, only BaseMesh to enable fast access
    bands::Vector{B}
    simplices::S
end

The construction of Bandstructures would have three steps

Step 1: diagonalize

The base mesh is assumed a BaseMesh{D}, following #122. We create a evals::Array{Vector{T},D} and evecs::Array{Vector{SubArray{2,T´}},D}, where T´ is the eltype of eigenvectors. These are left uninitialized. The size of the Arrays is the same as vertices(base).

Iterating over base vertices, (n, ϕs) in enumerate(vertices(base)) (here n is a CartesianIndex), we obtain ϵs, ψs = diagonalize(h(ϕs...), ...) on each. ϵs should be sorted. We store evals[n] = ϵs, making sure to set any degenerate ϵs[n]≈ϵs[n´] to be exactly equal ϵs[n]==ϵs[n´].

Regarding the eigenstates, we allocate and store in evecs[n] a vector of eigenstate views of the ψs matrices

evecs[n] = [stateview(ψs, rng, shift) for (rng, shift) in approxruns(ϵs)]

where

function stateview(ψs, rng, shift)
    v = view(ψs, :, rng)
    if !iszero(shift)
	    r = circshift(v, (0, shift))  # creates a copy
		v = view(r, :, size(r, 2))
	end
	return v
end

The output of approxruns(ϵs) should be (rng, shift) where rng is the range of indices j of ϵ around i in eachindex(ϵ) such that all ϵ[j] are approximately equal (≈). shift should be the offset of j into rng, i.e. shift = j - minimum(rng)

The interpretation of the above is the following:

• Each non-degenerate eigenvector in evecs[n] will correspond to a 1-column view into its corresponding eigenstate in ψs
• Each d-degenerate eigenvalue in evecs[n] will correspond to a d-column view into the corresponding subspace in ψs
• Each of the d repeated eigenvectors is a view into an independent copy of the full degenerate subspace in ψs, and with a different first eigenstate (thanks to the circshift). This is important for step 3.

Step 2: extract bands

We now call a version of extract_band adapted to the format of new evecs and evals.

It proceeds by incrementally classifying all evals and evecs into different bands, while at the same time pushing them into flat Vectors of BandMesh vertices and creating the corresponding adjacency matrix. We now revisiting how extract_band works, and adapt it to evecs and evals

We start by creating an array vertexindices::Array{Vector{Int},D} similar to evals that serves double duty, as a registry of unprocessed vertices (0), already processed vertices (-1) and in the process of being added to a new BandMesh (vertex index within the band). We initalize it with zeros.

We create an empty list of band mesh vertices verts and reuse I,J to construct sparse adjmat that will be populated to create a BandMesh.

We also create a first-in-last-out queue pending = [(srcidx, dstidx),...] that contains vertex indices of connected vertex pairs that have still not been added to the band. srcidx and dstidx are actually of type Tuple{CartesianIndex{D},Int}, and point to specific entries in vertexindices, evecs and evals. The pending queue starts empty.

We need to keep also a processed::Vector{Tuple{CartesianIndex{D},Int}} which stores the evals and evecs indices of each processed vertex in order. This processed is necessary to construct the overall set of Simplices in step 3.

extract_band algorithm:

• Create an empty processed, pending, I and J.

• (START) We identify the first non-zero element in vertexindices, with index seedidx. If none exist, break.

• Resize reusable pending, I and J to zero length. Create an empty verts vector.

• We push!(pending, (seedidx, seedidx))

• processed_offset = length(processed)

• While !isempty(pending)

  • (srcidx, dstidx) = pop!(pending)
  • We push! a new vertex dstidx into verts
  • We push!(processed, dstidx)
  • We register new vertex index with vertexindices[dstidx] = length(verts)
  • If srcidx != dstidx (i.e. not seed) we push! the adjacent vertexindices[srcidx] and vertexindex[dstidx] into I,J
  • We iterate over neighbors dstbase of srcbase = first(srcidx) in base
  • (proj, m´) = findmostparallel(evecs, dstbase, srcidx, vertexindices[dstbase]) will select the most-parallel target vertex evals[dstbase][m´]
  • If proj is above threshold and m´ != 0 we push!(pending, (srcidx, (dstbase, m´)))

• After this concludes (pending is empty), construct adjmat, build simpindices from adjmat, and push a BandMesh object into the bands vector.

• Set all non-zero vertexindices to -1

• Go back to (START)

I think there is currently a type-instability in the master version of this extract_band algorithm (srcidx is not type-stable?), so I wouldn't be surprised if this version would run faster. We also assume in master that all base vertices have the same number of eigenvalues, which need not be true if we are using Arnoldi methods with some libraries.

A key modification needs to be made to findmostparallel to handle Dirac points. We need to ensure that if a Dirac point occurs at a given base vertex (i.e. we have two copies of the Dirac 2D degenerate subspaces), the same copy will be selected from different angles within the same band. If a Dirac point occurs within a simplex, that will lead to a dislocation, and the corresponding simplex will be omitted (to be checked). So how do we do that?

We have multicolumn views evecs[n][m] for degenerate subspaces, so proj(ψdst, ψsrc) = sum(abs2, ψdst' * ψsrc). m´ should be the state with maximum overlap, unless its vertexindex is non-zero, in which case we return m´ = 0 . This sentinel value means "do not add any vertex from this dstbase column to the band". In the presence of degeneracies, however, different copies of degenerate vertices will exhibit the same overlap. If the maximum overlap has ties, an m´ = 0 should returned if either all tied vertexindices are -1 (dislocation) or if any of them is positive (already classified in this band). Otherwise, return an m´ with zero vertexindex.

Step 3: collect simplices

At the end of step 2 we will have constructed the list of base indices processed corresponding to band.verts for band in bands. We now need to use this, in combination with band.simpinds for band in bands to build sverts (simplex vertices), svecs (simplex states) and sptrs (mapping of base simplices to range of bandstructure simplices) that go into Bandstructure's simplices::Simplices.

The strategy is the following

• Create temporary arrays sverts::Array{Vector{NTuple{D´,V}}, D´} and svecs::Array{Vector{NTuple{D´,S}}, D´}, where V = eltype(first(evals)) and S = eltype(first(evecs)). We will collect all simplex vertices and vectors from all bands here.

• For each is::NTuple{D´,Int} in band.simpinds we obtain their band indices inds = getindex.(Ref(processed), is), where each element of inds is a (n,m)::Tuple{CartesianIndex{D},Int}, i.e. a base index and a column index.

• We extract the base vertex indices basevinds = first.(inds). These should correspond to a simplex in the base mesh, but with shuffled order. We do basesinds = (n0..., s) = basesimplex(basevinds) to get the indices of the simplex in the base mesh, and inds´ = sortsimplex(inds) to get full (n,m) indices inds into canonical order.

• We do svert = NTuple{D´}(vcat(base[n]..., evals[n][m]) for (n,m) in inds´) and push!(sverts[basesinds...], svert)

• The same with vectors: evec = NTuple{D´}(evecs[n][m] for (n,m) in inds´) and push!(svecs[basesinds...], evec)

• We do the above with all bands. Finally, we can count how many svecs and sverts, and copy them to contiguous vectors sverts´::Vector{Ntuple{D´,V}} and svecs´::Vector{Ntuple{D´,S}}, filling sptrs::Array{UnitRange{Int}, D´} along the way.

There is one important detail that is missing above. We have simply copied evecs[n,m] to svecs´, but these are SubMatrix, since they may come from degenerate subspaces. We need to project these onto SubVectors, with the specific direction in the subspace fixed simplex by simplex. Here is the core of the Dirac point resolution. The resulting svecs´ will thus be of type svecs´::Vector{NTuple{D´,S´}} where S´<:SubArray{T´,1}.

The first step is to take the least-degenerate subspaces in the simplex, and select the first column of each as the SubVector for that vertex. We denote by projection_inds their indices within svecs´. No copy is involved here, just a view of a view. Due to the circshift in step 1 we are sure that these SubVectors will all be orthogonal to other copies atop the same base vertex, so no redundacies should occur. Then, the orientation of the remaining vertices with higher degeneracy (if any) will be decided by the (now 1D) SubVector vertices. To avoid aliasing, these remaining states will be first copied, then projected, and then a view will be made of them, to preserve the eltype of svecs´. The projection will be made onto vsum = sum(i -> svecs´[:,i], projection_inds) |> normalize!, so as to maximize overlap with all projection vectors.

[pablosanjose]

There is a very important issue that can arise in the approach above: the interpolated set of eigenstates at a given ϕs does not form an orthogonal basis. I can imagine this could give problems when computing observables.

This doesn't mean that the general algorithm is invalid, it just might need to be modified. It essentially consists of two distinct elements: (1) degenerate eigenvalues carry a full basis of their degenerate space for purposes of connecting band vertices, and (2) simplices have eigenstates that do not need to be vertex-unique when compared with neighboring simplices.

To solve the Dirac-point issue, what we actually need is (1), which resolves frustrations in the connectivity. We could dispense with (2), in fact. Once connected, we can actually choose vertex-unique states in degenerate subspaces. They might simply not be too similar to the other vertex states in a given simplex with vertex degeneracies. Interpolation of eigenstates will thus not be very smooth. This is inevitable in Dirac points, but not in trivial band crossings, where velocity codiagonalization is optimal.

So I see two options to achieve a sane bandstructure with simplex-unique states: keep (1) and the current codiagonalization pass, dropping the last two paragraph in the OP above, or keep (1) and adapt the projection step, dropping the codiagonalization. I'm inclined towards the latter.

How would it go? First we would never do copies of degenerate vertex subspaces. This way we have a built-in guarantee that all vertices end up with unique eigenstates, because all transformations we do of subspace states are shared between vertices. We can transform these subspaces in any way we want and, as long as in the end we have different degenerate vertices pick at the end a different column of their space (which will remain orthogonal) we're good. So the crucial change, apart from never doing copies of states, is to know how to transform subspaces at simplex-collection time. The ideal would be to be able to have views with shuffled columns, and use the first one as the axis to align. Views don't support that, but we can create a new wrapper for views that holds the desired alignement column for each vertex (i.e. store the shift of the circshift above together with the view, without copying.).

So all steps above will proceed the same, except for replacing the circshift with a wrapper of a view basis and the shift, renamed to axis (and adapting proj to unwrap basis, of course). At the end of step 3, however, we wish to do the following:

• When populating svecs we will be processing simplices in the same order they appear in each band, which implies that they are order by the first vertex.
• Take the first least-degenerate subspace and reduce its basis to a view of the column given by its axis (no copies). This will be the reference vector vref to align the simplex
• We now iterate over the remaining subspaces (vertices in the simplex). For a given subspace, we compute the projection vproj of vref on its basis
• We build a permutation P_axis that reorders the basis, so that basis*P_axis simply switches the axis vector to the first column in basis
• We build a matrix A = hcat(vproj, basis*P_axis) whose first two columns are vproj and the axis vector. Its rank is the same as basis, since the added vproj is a linear combination of basis.
• We do a QR decomposition of a matrix A. The Q matrix obtained is the new orthonormal basis of the subspace we're after. Its first vector is proportional to vproj, and the rest are orthonormal, but not randomly so.
• We now just need to rewrite basis with Q, but first we have to move the first column back to position axis, which is done with the transpose of P_axis. We do basis .= Q * P_axis'. Note that this modifies the parent array that basis points to. We don't need to do anything else, as the basis, axis wrapper for the vertex in question is not modified, only the parent array is.
• Finally, we take a view(basis, :, axis) as the 1D SubArray for the vertex. Note that this is equal to vproj but the preceding steps are necessary to update also the rest of the basis, to keep it orthonormal for all vertices in the same vertex stack.

This is quite a nice method, because it guarantees a unique (and maximally continuous?) sampling of the eigenstate manifold, even in the non-Abelian case of degenerate . This might still have dislocations and jumps, but they will be unavoidable. The one loose end I see is the choice of vref on each simplex. I think it is fine, but I worry about the edge case in which the first vertex of a simplex has higher degeneracy than the rest, which are still degenerate. Then, vref might have a completely random orientation (if that vertex has not been visited and aligned before). Its orientation will thus propagate to any later simplex connected to the first, even if at a later step the random vref gets aligned with some other lower-degeneracy vertex. But I'm not sure what the consequences of this might be. I'll have to check.

[pablosanjose]

While implementing this proposal I've realized that the complexity of organizing band simplices into groups with the same base simplex (the sptr array) is tricky to implement clearly and efficiently. It is much easier to classify them simply by base minicuboid (i.e. little D-dimensional cubes made of D! base mesh simplices). This implies making sptr::Array{UnitRante{Int},D} instead of sptr::Array{UnitRante{Int},D+1}.

Marching meshes have the property that the minicuboid a simplex belongs to can be identified by the "lower left corner". That is, all simplices in a minicuboid share a vertex that can be found as the first in the simplex list of vertices when sorted lexicographically (this is not true of other meshings). Then, sorting simplices into minicuboids is straightforward. Accessing them with bs[(ϕs...)] is also straightforward, as one just needs to find the mincuboid that contains (ϕs...), which is trivial. The price is that one will need to check D! times more simplices to find the ones contributing to bs[(ϕs...)], but that is O(1), and it might actually be faster because, crucially, we don't need to keep and search an ImmutableDict to do translations of vertex groups to simplex indices in a given minicuboid. We do a linear lookup of all simplices in the minicuboid instead, which for D<=3 might be more efficient.