pablosanjose / quantica.jl Goto Github PK

When two modifiers touch the same element the last overwrites the first, e.g. parametric(h, @onsite!(o->o+1), @onsite!(o->o-1)) == parametric(h, @onsite!(o->o-1))

Currently, we can visualize a field, such as an eigenstate or the LDOS of a finite or unbounded system using e.g. vlplot(h, field), where field is a vector. Since an infinite system has a uniform LDOS (each unit cell has the same LDOS at a given energy), this is enough visualize the LDOS in this case. However, it doesn't allow us to visualize the LDOS of a semi-infinite system, or a system with boundaries. The PR #150 now allows us to compute the spatial Greens function in semi-infinite quasi1D systems, so it would be good to enable a way to visualize fields in semi-infinite lattices too.

A flexible way would be to allow a method vlplot(h::Hamiltonian{LA,L}, field::AbstractVector{<:Pair{NTuple{L,Int}}}) that accepts things like field = [(1,1) => [amplitudes at cell (1,1)...], (1,2) => ...]. Then, vlplot would plot the cells contained in field with their specific fields.

One problem with this is that there is no way to indicate where the boundaries are, or where to represent the shaded cells that represent the unbounded direction. Perhaps that could be solved by passing a boundaries kwarg similar to the one for greens, and have vlplot compute a bounding box that includes the boundaries. Another, more natural way, would be to include cells to be shaded inside field in some way.

However, I think it would be best to not overcomplicate things and forego shading of open cells. Ultimately, we want this to encode information about the field, of which we know nothing in unspecified cells, so nothing should be plotted there, probably. At most, we could allow an opacity directive that depends on the cell variable

I am trying to study the massless Dirac equation in a 2D square lattice. The model is defined as

lat2D = lattice(
    sublat((0, 0), name = :A),
    bravais = ((1, 0), (0, 1))
)

σ = (
    SA{ComplexF64}[0 1; 1 0],
    SA{ComplexF64}[0 -im; im 0],
    SA{ComplexF64}[1 0; 0 -1]
)

masslessmodel2D = 
    hopping(-im*σ[1], dn = (1, 0), plusadjoint = true) +
    hopping(-im*σ[2], dn = (0, 1), plusadjoint = true) 

dirac2D_massless = hamiltonian(lat2D, masslessmodel2D, orbitals = (:up, :down))

Evaluating the bandstructure raises the warning:

bandstructure(dirac2D_massless)

┌ Warning: Unexpected band connectivity between [0.0, -0.5235987755982988] and [0.5235987755982988, 0.0]. Rank 1 in (1, 1) projector.
└ @ Quantica /home/bamorim/.julia/packages/Quantica/5Pt0g/src/bandstructure.jl:637

This is a simple discritization of the Dirac equation in 2D:

This model as 4 Dirac cones located at (0, 0), (pi, 0), (0, pi) and (pi, pi). I believe the warning is emitted due to the fact that the Dirac cones are located at the edges of the Brillouin zone. If I use a shifted sampling of the Brilouin zone, e.g.:

bandstructure(dirac2D_massive, cuboid((-pi+0.3, pi+0.3), (-pi+0.3, pi+0.3)))

No issue is raised. Adding a mass term

dirac2D_massive = hamiltonian(lat2D, masslessmodel2D + onsite(0.1*σ[3]), orbitals = (:up, :down))

also removes the warning

This problem is not specific to 2D. The same happens in 3D. Considering the lattice Dirac equation in 3D:

α = ntuple(i -> kron(σ[1], σ[i]), Val(3))
β = SArray{Tuple{4, 4}}(kron(σ[3], I(2)))

lat3D = lattice(
    sublat((0.0, 0.0, 0.0), name = :A), 
    bravais = (
        (1, 0, 0),
        (0, 1, 0),
        (0, 0, 1)
    )
)

masslessmodel3D = 
    hopping(-im*α[1]; dn = (1, 0, 0), plusadjoint = true) + 
    hopping(-im*α[2]; dn = (0, 1, 0), plusadjoint = true) +
    hopping(-im*α[3]; dn = (0, 0, 1), plusadjoint = true) 

MasslessDirac3D = hamiltonian(lat3D, masslessmodel3D, orbitals = :A => Val(4))
MassiveDirac3D = hamiltonian(lat3D, masslessmodel3D + onsite(0.2*β), orbitals = :A => Val(4))

The massless model has 8 Dirac cones: (0, 0, 0), (pi, 0, 0), (pi, pi, 0), (0, pi, 0), (0, 0, pi), (pi, 0, pi), (pi, pi, pi), (0, pi, pi). And we get the same "Unexpected band connectivity issue". Once again, adding a mass or using a shifted Brillouin zone sampling makes the warning go away:

bandstructure(MasslessDirac3D) # Warning: Unexpected band connectivity between
bandstructure(MassiveDirac3D) # No warning
bandstructure(MasslessDirac3D, cuboid((-pi + 0.1, pi + 0.1), (-pi + 0.1, pi + 0.1), (-pi + 0.1, pi + 0.1))) # No warning

Even when the warning is raised, a Bandstructure is still returned. So far I am not able to determine if the results are reliable.

Currently, if we have a sp = Spectrum(ham), we can do energies(sp) to obtain the eigenenergies and states(sp) to obtain a matrix with the eigenstates. If we have a bs = bandstructure(ham, mesh), we can do bands(bs) to obtain a vector of Band objects, energies(bs) to obtain a vector of all energies of all bands and states(bs, i) to obtain all the eigenvectors of band i.

This is a bit of a mess. It also doesn't allow us to directly access one specific eigenvalue of a spectrum or bandstructure. If we want to plot the spatial profile of an eigenstate close to a given energy at a given momentum of bs, there is no easy way.

I propose redesigning the way we extract eigenvalues and eigenvectors from Spectrum, Band and Bandstructure by defining a convenient getindex indexing syntax. The following possibility comes to mind:

Spectrum sp

sp[1] -> ϵ₁, ψ₁  # first eigenpair in sp, with ϵ₁::Number and ψ₁::Vector
sp[1:3] -> ϵv, ψm # tuple of ϵv::Vector and  ψm::Matrix of first three eigenvalues and eigenvectors, respectively
sp[around = 0.2] -> ϵᵢ, ψᵢ  # single eigenpair closest to energy 0.2
sp[around = (0.2, 10)] -> ϵv, ψm  # 10 eigenpairs (encoded as ϵv::Vector and  ψm::Matrix) closest to energy 0.2

Band b

b[(ϕᵢ...)] = ϵ, ψ  # interpolated eigenpair at Bloch phases (ϕᵢ...)

Bandstructure bs

bs[1] -> band₁  # equivalent to bs.bands[1]
b2[2:3] -> bandvector # equivalent to bs.bands[2:3]
bs[1, (ϕᵢ...)] -> ϵ, ψ  # equivalent to bs.bands[1][(ϕᵢ...)]
bs[2:3, (ϕᵢ...)] -> ϵv, ψm  # ϵv::Vector and  ψm::Matrix generalization of the above for bands 2 and 3
bs[(ϕᵢ...), around = 0.2] -> ϵ, ψ  # interpolated eigenpair of the band with interpolated energy closest to 0.2
bs[(ϕᵢ...), around = (0.2,10)] -> ϵv, ψm  # interpolated eigenpairs of the 10 bands with interpolated energies closest to 0.2

The trickiest bit here is b[(ϕᵢ...)], which requires identifying the simplices which contain (ϕᵢ...) inside, and then using linear interpolation of eigenpairs on each. With our current definition of multivalued bands, there can be several such simplices, and as a result b[(ϕᵢ...)] can yield several eigenenergies and eigenstates. In order to keep this type-stable, we might need to make b[(ϕᵢ...)] a Tuple{Vector,Matrix}, even though the energy vector will most often have a single element. As a result, also bs[1, (ϕᵢ...)] will yield a Tuple{Vector,Matrix}. The only methods that may return a single Tuple{Number, Vector} are sp[1], sp[around=0.2] and bs[(ϕᵢ...), around = 0.2]. It might be worth considering making all these method return Tuple{Vector,Matrix} instead.

I'd like to revisit the cancelled #105 and its successor #106.

In #105 we considered internalizing the similarmatrix(ham, blochtype) matrix into ham::Hamiltonian, so as to make it less cumbersome to produce Bloch matrices repeatedly without allocations. Currently we need to do

bloch!(similarmatrix(ham, blochtype), phis, [axis])
bloch!(similarmatrix(pham, blochtype), (phis, params), [axis])

But #105 turned out to be a bad design for the reasons stated in #106 (essentially, the choice of blochtype might depend on context, which forces one to regenerate ham to change blochtype), so #105 was cancelled. Instead #106 simply enhanced the powers of similarmatrix, but left the Bloch matrix as external.

However, there is perhaps a good alternative to #105 that we could consider, with additional user-facing advantages. We could define a new type

struct BlochMatrix{H<:Union{Hamiltonian, ParametricHamiltonian},M<:AbstractMatrix}
    matrix::M
    h::H
end

bloch(h::Union{Hamiltonian, ParametricHamiltonian}, blochtype = missing) = BlochMatrix(similarmatrix(h, blochtype), h)
bloch(blochtype = missing) = h -> bloch(h, blochtype)

(bm::BlochMatrix{Hamiltonian})(phis...) = bloch!(bm.matrix, bm.h, phis)

(bm::BlochMatrix{ParametricHamiltonian})(phis...; kw...) = bloch!(bm.matrix, bm.h(;kw...), phis)

This could allow us to do the following

m = lattice |> hamiltonian(model) |> bloch(flatten)
m(0.2, 0.3) ## produces the bloch matrix for phis = (0.2, 0.3) in-place
m = lattice |> hamiltonian(model) |> parametric(modifiers) |> bloch(flatten)
m(0.2, 0.3; params...) ## produces the bloch matrix for phis = (0.2, 0.3) and params in-place

The base spectral methods such as bandstructure, spectrum or even KPM could then take a BlochMatrix directly, which already knows how to produce the matrices needed. For convenience we could also define the current shortcuts in terms of Hamiltonians.

We could be tempted to make BlochMatrix<:AbstractMatrix and use the m object directly in any function that expects a matrix, but that would be inelegant, as the value of the actual matrix would depend on the last call to m(phis...), which can be highly non-local. Better to just use m(phi...) as explicit input instead.

Current is a slightly tricky thing. To compute the current carried by a state we need to use the continuity equation ∂t ρ(r) = -∇⋅j(r). So, the current operator depends on the Hamiltonian, which is what drives time evolution. Even more, the current operator depends on the gauge used to write the Hamiltonian. This can be seen from the form of the j current flowing along a hopping link kl in a lattice. Using the continuity equation with zero vector potential A=0 we find j_kl = imag(⟨k|psi⟩t_kl ⟨psi|l⟩). Note that if we change gauge, a vector potential A will arise that will affect the phase of t_kl (Peierls), so the current operator changes. The expectation value of the current, however, is unaffected by the change of gauge, because eigenstates psi also acquire a phase that compensates. The current is a gauge-independent observable.

The question I'd like to answer here is how to visualize current using vlplot or plot. One way is straightforward: weight the color/opacity of each link with j_kl. But in systems with lots of hoppings this can be quite a mess to understand. Another option would be to show the current flowing across each site. We could then say

vlplot(h, psi; sitecolor = current(axis))

where current(axis) should then resolve with h and produce something like (psi, i)-> sum_k (r[k][axis] - r[l][axis]) imag(t_kl psi[k]' psi[i]), and the sum is over sites k linked to site i

The problem is that currently vlplot site shaders only accepts functions of psi[i]. To plot the current, then, we have two options.
Option A: pass j instead of psi to vlplot, with something like

vlplot(h, current(psi, h); sitecolor = psi_vec -> psi_vec[axis])

Option B: allow sitecolor and other site shaders to accept a special object current(axis) that does the right thing once h is known

vlplot(h, psi; sitecolor = current(axis))

The advantage of B is that we could plot density and current simultaneously (using e.g. the siteopacity and sitecolor channels).

This current(axis) could be generalized so that current() represents the modulus of the current. We could also allow this quantity in link shaders. For completeness we could also define a density() object.

When dealing with very large systems, it might be undesirable to materialize the system in memory when using the KPM functionality. In these cases other specialized packages (e.g. Kite) use matrix-free approaches, by implementing matrix-vector multiplications without storing the matrix in memory.

We have the underpinnings to do this quite easily in Quantica. We can use a Hamiltonian{<:Superlattice} for this.

The idea here could be something like this

• Move the 'modifiers' functionality currently in unitcell into hamiltonian(::Superlattice), so that a Hamiltonian{<:Superlattice} can incorporate the modifiers without materializing into a Hamiltonian{<:Lattice}
• Create a new type of lazy LazyBlochMatrix::AbstractArray, to be passed to bloch!(::LazyBlochMatrix, h::Hamiltonian{<:Superlattice}, phases), that incorporates the base harmonics of h, any modifiers, any bloch phases if present, and the superlattice
• Implement a 5-argument mul! method for LazyBlochMatrix that implements the product on a vector as efficiently as possible, preferably using threads and SIMD. This can potentially even be faster than the generic matrix-vector multiplication because of the fixed banded structure of all Bloch Hamiltonians of superlattices.

The onsiteselector!(model; ...) -> model and hoppingselector!(model; ...) API was introduced in #18. This issue is prompted by some problems with this syntax

1. Since model::TightbindingModel is actually not mutable, what is being done here is to contruct a new model with its onsite or hopping terms modified according to the specified selectors. Hence the ! is misguided.

2. There is a case to be made that adding yet another name to the (quite large) API is a bad idea. Could we reuse some other API function for this?

A possibility for 2: do onsite(model; kw...) instead of onsiteselector!(model; ...). @fernandopenaranda comments on Slack that this usage could potentially be confused with onsite!, which is a ElementModifier. However, currently onsite always returns a model. This new method would keep that property.

We need a way to visually represent eigenstates on the lattice, both using Makie and VegaLite.

An appropriate API would be plot(h, siteradius = abs2.(eigenstate)), or plot(h, sitealpha = abs2.(eigenstate) or plot(h, sitecolor = abs2.(eigenstate), to be able to fine-tune the result.

An additional method plot(h, eigenstate) could be an alias of plot(h, sitealpha = abs2.(eigenstate) or whichever looks nicer.

Question: what should eigenstate be? I think it should be simply a Vector, and could be extracted from a bandstructure or spectrum with states(::Spectrum, level = level_number) and states(::Bandstructure, ϕs, band = band_number) (this API does not yet exist). We would then need to check if the eltype and dimension matches h. Another possibility is to bundle h into bandstructure and a new Ket struct (we currently have a KetModel but not a Ket). However, it probably makes everything too inaccessible. Having eigenstates as Vectors allows the user more familiar ways to manipulate them directly.

We wants it.

Currently, SiteSelector and HopSelector allow sublats, dn and region keywords to specify sites.

When we want to specify individual sites, this is not very flexible. Possible use cases: select exactly two sites close to a point to find the local DOS in them, or a single one to add an adatom, etc...

We could introduce a kwarg siteindices to SiteSelector and HopSelector, like

siteselector(; region = ..., siteindices = 1:4)  # select among sites with index 1 to 4 in region
hopselector(; region = ..., siteindices = ((1,3), (2,4)))  # select among site pairs in region

Here siteindices should be an integer or interable of integers, and those integers should refer to indices in sites(lat) directly. Incidentally, we should rename sites(lat) to be sitepositions(lat) and extend it to allow a selector too, sitepositions(lat; kw...) = sitepositions(lat, siteselector(;kw...))

To aid in providing siteindices to a siteselector we should also introduce a siteindices function (it can have the same name as the kwarg, that's not a problem), analogous to sitepositions, that returns the site indices of a selector siteindices(lat; kw...) = siteindices(lat, siteselector(;kw...)). We should guarantee that the order of the returned indices matches the positions returned by sitepositions.

We should also have siteindices and sitepositions defined on Hamiltonians, not only Lattices

With this we could do things like [ket(1; siteindices = s) for s in siteindices(h, region = myregion)] to obtain a set of single-site kets in myregion.

Checking whether a very large Hamiltonian is Hermitian is inefficient, particularly for compressed-column sparse Hamiltonians which do not favor row access. Currently, we have no way to encode Hermiticity, and this is unfortunate.

We could add a field to a Hamiltonian that signals hermiticity when the original model is hermitian, and then force setindex! to preserve hermiticity. Alternatively we could wrap it in Hermitian. The latter would be the more Julian solution, but then we should make Hamiltonian satisfy the AbstractMatrix interface, which is probably not possible, as it is a more general object.

This issue is a reminder and an invitation to discuss possible designs.

A bandrange calculation, for instance in momentaKPM() (when bandrange kwarg is missing), using bandrangeKPM() can't be performed if the hamiltonian h isn't flat. _partialshur() in ArnoldiMethod throws the following error:
julia> momentaKPM(h)
ERROR: MethodError: no method matching complex(::Type{SArray{Tuple{4,4},Complex{Float64},2,16}}

I guess that the easiest fix without loss of code generality/flexibility is to simply introduce a h |> flatten() in bandrangeKPM(). What do you think?

I've found that under some circumstances (I cannot be more precise at the moment) we cannot generically access all subspaces or subarrays of a given bandstructure by means of the built-in indexing syntax. I provide here a not-so-minimal failing example:

using Arpack
const σyτz = @SMatrix[0 -im 0 0; im 0 0 0; 0 0 0 im; 0 0 -im 0]
reg = Quantica.RegionPresets.Region{2}(r ->  0 < r[2] < 100)
lat = unitcell(LatticePresets.honeycomb(), (1,1), region = reg)
h = lat |> hamiltonian(onsite(I) + hopping(I) + onsite(σyτz), orbitals = (:eup,:edown,:hup,:hdown))
b = bandstructure(h, 
                   cuboid((-.5, .5); subticks = 5), 
                   mapping = x -> 2pi*x, 
                   method = ArpackPackage(nev=16, maxiter = 300, sigma=0.001im));
b[(0.25, ), around = 0]

b[(ϕ, ), kwargs...] fails to select the subpace with ϕ = 0.25. Moving away from the vertices by a small amount, ϵ, as @pablosanjose suggested in Slack does not generally fix this issue.

It seems that the type of aParametricHamiltonian field cannot be inferred from the input arguments in parametric(). For instance:
julia> @code_warntype LatticePresets.honeycomb() |> hamiltonian(onsite(0)) |> parametric(@onsite!((o;μ) -> o- μ)),
reveals a type instability.

As a result, the benchmark of parametric() leads to what I think are avoidable allocations (correct me if I'm wrong). This extra memory allocation due to the dynamic dispatch increases with the system size, being negligible in small to moderate size systems, but significant when large hamiltonians (>1 million sites) are built.

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?

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.

In the present documentation of the bravais function, the following example is presented:

Examples
julia> bravais((1.0, 2), (3, 4))
Bravais{2,2,Float64} : set of 2 Bravais vectors in 2D space.
  Vectors     : ((1.0, 2.0), (3.0, 4.0))
  Matrix      : [1.0 3.0; 2.0 4.0]

However, there is no method of bravais that takes as arguments tuples, vectors or static vectors. bravais has the following methods:

# 3 methods for generic function bravais:
bravais(lat::Quantica.AbstractLattice) in Quantica at /home/bamorim/.julia/packages/Quantica/5Pt0g/src/lattice.jl:171
bravais(h::Quantica.Hamiltonian) in Quantica at /home/bamorim/.julia/packages/Quantica/5Pt0g/src/hamiltonian.jl:915
bravais(ph::Quantica.ParametricHamiltonian) in Quantica at /home/bamorim/.julia/packages/Quantica/5Pt0g/src/parametric.jl:191

Is this intended? For me it would make sence to use bravais to define a bravais basis, which can then be used to define an Hamiltonian.

When computing a KPMobservable, in addition to the stochastic procedure the current version of KPM.jl allows to specify the subset of kets over which the trace is evaluated through a kwarg called ket. This is very handy in the following situations:

1) If some kets won't give a finite contribution to the trace and, therefore, we should restrict the randomket selector to the subspace of non-null contributions. As a result, a performance enhancement is achieved.

2) If we want to replace the random ket selector procedure by an exact evaluation of the trace over all basis kets. This may be necessary if the stochastic approximation leads to poor precision, for instance, due to a small system size.

3) If we combine an exact trace evaluation ( 2)) over the non-null contribution subspace (1)).

The KPM techniques already implemented in src/KPM.jl allow to cover the previous situations by passing a ketset::AbstractMatrix to momentaKPM() with the pointers to the basis ket subspace. In my opinion the API could be greatly simplified if for instance point 1) above is addressed as a default in the computation of momentaKPM() and if a new kwarg let's call it method is added. method may either be random (default) or exact to deal with the cases of 2)). Does this API makes sense?

As pointed out by @pablosanjose in Slack, it might be worth rethinking the way a ketset is passed to momentaKPM().

We currently have no way to constrain hopping(t; range =...) with a lower bound on range [EDIT: we can do it with region]. Allowing range = (min, max) could be useful. Currently, the NearestNeighbor library doesn't allow specifying a minimum range, so we would have to rely on filtering out what NearestNeighbors.inrange returns using min range.

CC: @BacAmorim, who came up with the idea

Currently the @onsite! and @hopping! macros don't understand a , after the function. Could be improved

julia> @hopping!(t -> 2t, sublats = :A=>:A)
ERROR: LoadError: ArgumentError: Only @hopping!(args -> body; kw...) syntax supported Mind the `;`.

It does work with a ; of course

julia> @hopping!(t -> 2t; sublats = :A=>:A)
Quantica.HoppingModifier{...}

There is a serious design problem that has been inadvertently introduced in #52. By combining parameters and Bloch phases into a single (params..., ϕs...) tuple, to be used in bloch, bandstructure, mesh etc.., I did not consider a use case where params are not Numbers (ϕs are always Reals).

Indeed, it is not unlikely that a ParametricHamiltonian has parameters that are vectors, matrices, or more complicated types. In such case, (params..., ϕs...) is not an NTuple{N,Number}, and hence cannot be converted to a SVector of concrete eltype (which is what we internally use for Bloch phases) if there are both params and ϕs.

The problem runs deep, because when building bandstructures, we need to build a mesh where each vertex is SVector(params..., ϕs...). We further need a norm to be defined on these vertices, because we use distances in mesh space for co-diagonalization.

A possible solution is to only allow non-numeric params in the case where there are no ϕs (i.e. in zero-dimensional ParametricHamiltonians). A less drastic approach is to allow them for creating ParametricHamiltonian, but not when calling bloch or buildmesh. Even with such restriction, a notion of distance is needed for non-numeric types, which probably cannot be defined in general.

Another possible solution is to not conflate parameters and Bloch phases at all, and keep them separate somehow. This entails a thorough revision of recent PRs that enabled bandstructures as a function of parameters (i.e. a serious amount of work).

This might be useful for KPM calculations:
https://github.com/jagot/ThreadedSparseArrays.jl

While not as efficient as a specialized method as in #66, it might be of use.

This is a followup of the discussion in the cancelled PR #76. I summarize the relevant posts here:

@BacAmorim said:

... the prelift/postlift destinction (even if not directly exposed to the user) still seems unnecessary to me. For this particular example, it would seem more natural for the user to just specify a linearmesh(0, 1) and then provide a lift = x -> (SVector(x[1], x[1]), x[1], x[1]).

Maybe my confusion is due to the name lift as I would expect this to map a lower dimensional space to a higher dimensional one. With the present proposal, this might not be case, as we can have just a repackaging of the params/bloch phases to a form suitable to the parametric Hamiltonian. For example, going back to the case of 2D system with electric field. We want to obtain the bandstructure in the 4D (Ex, Ey, kx, ky) space. We would use marchingmesh((0, 1), (0, 1), (0, 1), (0, 1)) to create a 4D mesh. With this proposal we would then write lift = (Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky), so a map between two 4D spaces. So maybe a more descriptive name would be parametrize = (Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky) or mesh2params = (Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky) ?

@pablosanjose replied

As usual, you make excellent points, so I'm going to delay merging this, and give us some time to elaborate.

In particular, you are implicitly proposing an alternative to the whole cut/lift approach. You say, have a single kind of mesh, parametrized by dimension (we could just make do with marchingmesh, and call it simply mesh(...; type = :marching). We would then implement cuts fully through a single parametrize or mapping or mesh2params kwarg interface (I personally like the name mapping). Then, all bandcut functionality is concentrated on the mapping, and Issue #75 is solved through a single mapping function. So, following your reasoning, these could be some usage examples

Full band, with an L-dimensional h::Hamiltonian

bandstructure(h)
bandstructure(h, mesh((0, 2pi),... ::Vararg{Tuple,L}; points = ...)) 

Straight bandcut, with an L-dimensional h::Hamiltonian

bandstructure(h, mesh(0,1), mapping = k -> (k, k,... )::NTuple{L})

Piecewise bandcut along a 1D cut from k = (0,0) to k = (0,pi), with an 2-dimensional h::Hamiltonian

bandstructure(h, mesh(0,1), mapping = piecewise((0,0), (0,pi)))

Piecewise bandcut along a 1D cut with symbolic nodes

bandstructure(h, mesh(0,2), mapping = piecewise(:Γ, :K, (pi,pi)))

General parametric Hamiltonian ph, hand-written mapping

bandstructure(ph, mesh((0,1), (0,2)), mapping = (x, y) -> (SA[x,x+y], y, y^2))

Parametric mapping from a 1D mesh, with the ph(; E) in the OP, and a sweep from E = SA[0,0] to E = SA[1,1]

bandstructure(ph, mesh(0,2), mapping = piecewise((SA[0,0], 0, 0), (SA[1, 1], 0, 0)))

The piecewise function could produce a function that linearly interpolates between arguments (using broadcast .+ and .*), as the arguments. If given purely symbolic nodes, it could build a lazy object that only resolves into a concrete function when the h becomes known.

Please let me know how you feel about this alternative. It does indeed seem simpler... If you agree, I might close this and start a new PR where we can polish the details.

@BacAmorim replied

I really like this! I also like the name mapping!
I think that this is single unifying API that appears to be quite natural and solves several problems in one go. I imagine that someone that is trying to learn Quantica (besides me) will have no troubles understading this feature.

One question: with the piecewise mapping, now are you going to split the points of the 1D mesh through the different segments? As equally as possible between the different segments, based on the length of each segments? (While I like this and have proposed it before, this would require defining a metric for arbitrary parameters, which was one of the initial problems here). linearmesh solved this by allowing explicitly specify the number of points per segment. How would this be addressed here? Maybe keep linearmesh and marchingmesh? Also, with the piecewise option, there is really no difference in writing

mesh(0,1), mapping = piecewise((SA[0,0], 0, 0), (SA[1, 1], 0, 0))

vs

mesh(0,2), mapping = piecewise((SA[0,0], 0, 0), (SA[1, 1], 0, 0))

The only thing that matters is the number of points in the segment. Correct?

I would suggest keeping in the API both linearmesh and marchingmesh (maybe just calling it mesh as the the marching part is an implementation detail). But I I am sure you will come up with a more elegant solution.

The algorithm in #123 has a deep flaw: not all possible simplices from a given band-mesh adjacency matrix are valid simplices if there are vertex degeneracies. There is a critical property that is not always satisfied by #123: at any given base mesh point, the bandstructure of a Hamiltonian h must have at most N = size(h, 1) simplices. This is briefly mentioned in #134

Here's a sketch of a case with N=3 where the adjacencies are correct but lead to 4 simplices atop a single base-mesh simplex

The number by each vertex is the degeneracy.

The following test, currently broken, showcase this issue

    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(2)
    bs = bandstructure(h, subticks = 13)
    @test_broken sum(degeneracy, bs[(0.2,0.3)]) == size(h,1)

Fixing this will likely require reengineering the whole bandstructure_knit approach, since it is fundamentally impossible to build the correct simplices knowing only the adjacency matrix. We will probably need to turn #123 upside down: compute the N simplices atop each base simplex first using eigenstate knowledge, and then build the adjacency matrix from the simplex set... perhaps. Computing simplices column by column like this will probably also accelerate the current bottleneck of band knitting, which is the many projections between subspace basis, which are small views of the full eigenstate matrices. Doing full-eigenstate matrix projections instead is likely to benefit much more from cache locality and BLAS optimizations.

As noted by @BacAmorim, it is good API design to have a single way to do each basic task. Simpler to learn, simpler to remember. This issue is a list of possible redundancies in the Quantica API that it would be better to simplify. I will update this post as these cases are found

• hopping(1, sublats = :A) == hopping(1, sublats = (:A,:A)) == hopping(1, sublats = :A => :A)

  I would favor to keep only the last form, or a tuple like sublats = (:A=>:B, :B=>:C) if there are more than one.

  Fixed in #29

• onsite(1, sublats = :A) == onsite(1, onsiteselector(sublats = :A))

  And similarly for hopping. I would favor keeping only the first in the API, and not export onsiteselector or hoppingselector

  Fixed in #23

• bloch(h, (1,2)) == bloch(h, 1, 2) == h(1,2)

  I would favor to keep only the first form. It involves one less splatting operation when sampling a list of ϕs, and is also more extensible in case we need to add other arguments (e.g. a method to generate the preferred similarmatrix(h, method))

  Fixed in #27

• bandstructure(h, npoints = 13) == bandstructure(h, npoints = (13, 13, 13)) for a h::Hamiltonian in 3D

  (Note that the npoints API is the one in #17. In current master it is called resolution). I would favor the second form, perhaps only leaving the first as a convenience for the 1D case.

  Fixed in #17: this method now accepts only a resolution = r kwarg, which is equivalent to npoints = (r, r, r). The npoints kwarg is removed. For uneven resolutions or arbitrary ranges, the user now needs to build the mesh explicitly.

Since in #17 we consider the possibility of parameter spaces different from Brillouin zone, it no longer makes sense to default to intervals -π:π for meshes. We should make marchingmesh more agnostic, accepting only ranges, and perhaps have bandstructure kwargs specify ranges also.

This issue aims to brainstorm the design of coupled Hamiltonians (let's call them BlockHamiltonian for the moment) so that it can cover a wide range of possible geometries. It's sometimes useful to start from the end: work out what you want to solve in sufficient detail, and then build backwards towards the fundamental design. In this first post I'll try to make a list of such problems.

(1) Standard multiterminal scattering problem

One 0D hamiltonian (central) coupled to a number N of 1D hamiltonians (leads). The leads need to be semibounded, with the zero cell located at an arbitrary point. The coupling is only between each lead and the central part, not between leads

This suggests we need to keep semibounded as a property of Bravais, and to add an origin field in Unitcell, and a function cellposition(lat, cell::SVector{L,Int}) = origin + bravais*cell

(2) Multiterminal scattering problem with coupled leads

A variation of (1) in which we allow couplings between leads. Arguably one can always reduce this problem to (1) by enlarging the central part sufficiently, but that may be very inefficient if the central system needs to grow too much.

One possible example is sketched here, where number n denote a cell of subsystem n and ⋱ denotes repetition:

11
11 22222222222222⋱
11         3
           3
           3
           3
           ⋱

This suggests we need a model for BlockHamiltonian in which all blocks h_ij (between system parts i and j) can be made non-zero, either a Hamiltonian (for i == j) or a Coupling (for i != j).

This also suggests that Couplings between Hamiltonians need to encode both a source cell and a destination cell, as these may or may not be the zero cell.

(3) Multiterminal scattering problem with higher-dimensional leads

Leads need not be one-dimensional. In general we may want to couple N Hamiltonians with arbitrary dimensions L_i. As an example, we might want to couple a semibounded nanorribon (L_1 = 1) to a half-plane (L_2 = 2).

This suggests that giving one of the Hamiltonians a privileged status as central is a bad idea. Our general goal will be to compute the Green function between any pair of sites and cells, each of them located on any subsystem.

(4) Leads with "holes"

One type of system that does not fall within the class of "unbounded or semibounded hamiltonians with arbitrary dimension and arbitrary pair-wise coupling" is e.g. a semibounded nanotube in the z direction that merges into a 2D graphene in the x-y plane. The latter is not a standard unbounded lattice, because it has a hole where the nanotube is coupled. This might be an unusual problem to pose, but how would one describe it?

This suggests the need to have unbounded Hamiltonians that may implement holes at a given point. In other words, we need support for vacancies in Lattice. To make vacancies as modular as possible we could do a WithVacancies{<:AbstractLattice} <: AbstractLattice, and have Green functions implement the effect of vacancies at some level.

(5) Hamiltonians with unbounded couplings

Where does the case of twisted bilayer graphene or an infinite coupled ladder lattice fit in? These cases are characterized by conmensurate lattices between which an unbounded coupling is set up. All other examples (1-4) have couplings that are restricted to a finite set of cells.

This suggests that such case does not fall within the BlockHamiltonian category, but rather should be treated separately, as is currently done by combine (which checks for commensuration and then produces a new Hamiltonian, not a BlockHamiltonian). By separating these cases clearly we avoid conflating concepts and adding unnecessary complexity.

CC: @fernandopenaranda, @BacAmorim

Just a quick thought:
Would it make sense to add a method to the current combine function so it can accept two Quantica.RegionPresets.Region and return the intersection (&&) of both? Something like:
Quantica.combine(reg1::Quantica.RegionPresets.Region, reg2::Quantica.RegionPresets.Region) = Quantica.RegionPresets.Region{2}(r-> reg1.f(r) && reg2.f(r))
And then generalise it to an arbitrary number of regions as inputs.
Although this is not strictly necessary since there are other ways to do it, I think it might contribute to lighten a bit the syntax in a large project in which long inline definitions should be avoided.
What do you think?
(It could also be interesting to create, in addition, two new functions that perform the union and the symmetric difference of two or more Regions, respectively)

This issue is a follow-up discussion of #17.

Thanks to #25, we now have ph::ParametricHamiltonians that know their parameters. As suggested by @BacAmorim, this opens up the door to implementing, generally, a bloch(ph, ϕs) that merges Bloch phases and parameter values into the same Tuple/SVector ϕs. It would in turn enable to write a bandstructure(::ParametricHamiltonian) over a mixed Brillouin zone/parameter space with minimal code complexity.

Mesh cuts

Typically, we will likely also want to perform cuts in this space, fixing one/several parameter or phase, or doing a more complicated constrained sampling. This can now be done through bandstructure(::Function) (#17), but it is slightly cumbersome to do efficiently, as the user needs to worry about allocating the work matrix for bloch! manually and building a matrixf function that implements the cut.

We could then add a third method bandstructure(ph::ParametricHamiltonian[, mesh]; cut::Function = missing) that allows us to do e.g.bandstructure(ph, marchingmesh(range(-π,π,length=10), 0:0.1:1); cut = (p,ϕ) -> (p, ϕ, -ϕ)) on a single-parameter ph

EDIT: for consistency, we should then also allow cuts with h::Hamiltonians, e.g. banstructure(h, marchingmesh(range(0,π,length=10)); cut = ϕ -> (ϕ, 0))` for a Γ-X cut.

One aspect of the API that we would need to fix is the order of parameters and bloch phases in the ϕs of bloch(ph::ParametricHamiltonian, ϕs). My gut feeling is that having parameters first and phases second, ϕs = (p₁, p₂,..., ϕ₁, ϕ₂, ...), is the most natural, since parameters need to be applied first to ph to obtain a Hamiltonian, which then takes the Bloch phases to give a matrix to diagonalize.

About the ParametricHamiltonian call syntax

In any case I think I would not overload the call syntax ph(;params...) of ph::ParametricHamiltonian, only bloch/bloch!, so that to obtain the Hamiltonian matrix of ph we would need to be explicit, and do bloch(ph(p=1), (ϕ₁, ϕ₂,...)). I say this because there exists the tempting alternative to have h(ϕ₁, ϕ₂,...) return the Bloch matrix of Hamiltonian h (a form that we recently had but that I deprecated in #27 to reduce API redundancy) and equivalently do hp(ϕ₁, ϕ₂,...; p₁ = 1, p₂ = 2...) to get the Bloch matrix of a ParametricHamiltonian ph. However, that clashes with the hp(p₁ = 1, p₂ = 2...) syntax to obtain a Hamiltonian, which would become ambiguous in the 0-dimensional case.

Followup of #120 (comment)

A D-dimensional Bandstructure is currently just a collection of meshes (bands) in D+1 dimensions (k+energy vertices), a matrix of states for each vertex (flattened into a vector) and a vector of simplex pointers into vertices.

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

struct Band{D,M,A<:AbstractVector{M},MD<:Mesh{D},S<:AbstractArray{<:NTuple}}
    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`
end

struct Bandstructure{D,M,B<:Band{D,M}}   # D is dimension of ε + parameter space
    bands::Vector{B}
end

This has several problems. In particular, simplices do not keep track of their cartesian indices on the base k-space mesh, so finding a subset of simplices around a given momentum becomes costly (we don't have higher-dimensional interval trees in Julia yet).

In this issue I want to consider a different approach to store the bandstructure that (1) solves the above access problem, while (2) still allowing efficient iteration over all simplices.

A third consideration to keep in mind is how to solve the problem of Dirac points. A Dirac point is a 2D degenerate subspace that has no natural basis, since velocity eigenstates depend on direction. Hence, the per-vertex codiagonalization we now do cannot yield simplices all around the Dirac point (half of them are always missing). A solution considered in the past is point-splitting: duplicate the Dirac point in the mesh, and choose a different basis in each copy. An alternative solution, incorporated in the present proposal, is to associate eigenstates to simplices, instead of vertices. To avoid the overhead of keeping multiple copies of eigenstates (since simplices share vertices), we resort here to storing state views inside simplices. Since we can store any view we want, we just need to store two different state bases for a Dirac point, and keep views to them in the relevant simplices.

So, a sketch of the proposed alternative implementation would be

struct CuboidMesh{D,TE}  # marching tetrahedra mesh of a D-dimensional cuboid in (ϕs...) parameter space
    axesticks::NTuple{D,Vector{TE}} # axes tick points of the form ([ϕ1...], [ϕ2...], ... [ϕD...])
end

struct Simplex{D,TE,TS,S<:SubArray{TS,1}}  # a single simplex with D+1 vertices (energy eigenvalues vals, and eigenstates vecs)
    vals::Tuple{TE, Vararg{TE,D}}   # D+1 vertices for simplices with a D-dimensional base space
    vecs::Tuple{S, Vararg{S,D}}     # vecs is a tuple of views into `states` below, or transformed copies thereof for degenerate vals
end

struct Bandstructure{D,TE,TS,SI<:Simplex{D,TE,TS}}
    base::CuboidMesh{D,TE}
    energies::Array{Vector{TE},D}   # This stores the eigenvalue output of diagonalize itself, no copies
    states::Array{Matrix{TS},D}     # This stores the eigenvector output of diagonalize itself, no copies
    simplices::Array{Matrix{SI}, D} # a column of simplices on top of each base simplex
    connectivity::Matrix{Int}
end

The reason for the Matrix eltype of simplices above is that for each "minicuboid" in the D-dimensional base mesh there is a fixed number N = D! of marching tetrahedra simplices inside, and for each of these we store a column of N´ simplices that form the lifted (bandstructure) mesh. This N x N´ would form the Matrix. It would be worthwhile to think if this could be improved further. For example not all N base simplices need have the same number N´ of simplices in the column above them (since the number of eigenvalues on each vertex need not be the same in general).

Perhaps a better way would be to have simplices::Array{Vector{SI}, D+1}, with the extra dimension spanning the N=D! simplices in each minicubiod. But then D+1 should be an extra parameter of the type, since we don't have computed type parameters in Julia yet.

CC: @fernanopenaranda

In the same vein as the shaders in vlplot(h::Hamiltonian, psi), we need to be able to encode eigenstate information into vlplot(::Bandstructure) bandstructure plots. One of the few ways VegaLite allows to do this is to use the trail mark and encode the information into its thickness, as in https://vega.github.io/vega-lite/examples/trail_color.html

Additionally, we could enable pointcolor/pointopacity/pointsize shaders, similar to site shaders in Hamiltonian plots, by plotting dots at band vertices, on top of the band curves themselves, a bit like this

Its probably best to first finish with the fixes for #135 that involve a Bandstructure refactor. In particular we should store eigenstates at band vertices inside bands::Vector{<:BandMesh}, which we currently don't do (we just store them in sbases, collected by simplices, but not classified into bands)

Continuing the discussion in pablosanjose/Elsa.jl#22 (comment)

So two orbitals (say :px and :py) would be labeled as different i's, even if they belong to the same atom?

Correct, in the list of hoppings that is outputed by Wannier90, the orbitals :px and :py would be labelled as 1 and 2, for example, independently of whether they are centered at the same or different sites.

If that is the case and if it is not possible to extract the atom information from Wannier90, I guess the best approach would be to create one Elsa site per Wannier90 orbital, even if some of them share the same point in space.

One of the outputs of Wannier90 is the position of the centre of the orbitals, therefore we can see if two orbitals share the same centre or not. So we could simply apply the rule 1 orbital = 1 sublattice (simpler but less elegant) or extract the information regarding the location of the different orbitals.

Then, we would implement setindex!(::Hamiltonian), so that we can do h[(i, j), dn] = -2.7, where i,j are now unitcell sites, and dn = (nx,ny,nz). That is by far the cleanest way to do this, I think, and could even make it possible to interact with Julia's broadcast mechanism for added convenience: h[wannier_indices] .= wannier_values

I think this notation would be really clean! Another thing that could be done is to simply read the wannier90 list of hoppings, from them build the vector of HamiltonianHarmonics and directly build the Hamiltonian with that, not using HoppingTerms or OnsiteTerms. Can it be done in that way? I am still not sure how the Hamiltonian is currently stored.

But in all honesty, I am not a Wannier90 user. I just had a brief experience with using hamiltonians that are outputed by Wannier90. All I know about Wannier90 is thanks to @joselado

This issue comes after the conversation in Slack with @pablosanjose.
Dangling bonds forming in a lattice is, most of the time, an inconvenient occurrence when building a lattice since they are not expected to be formed in an experiment and can lead to spurious wave-function localisation at these sites.
Avoiding dangling bonds formation in a lattice of arbitrary shape and dimension is generally a difficult task following the two currently available approaches, either using a composite Quantica.RegionPresets.Region or the kwarg indices in siteselector.
Both are troublesome since a prior knowledge about the exact indices (possibly many) in which the dangling bonds will be formed is required and, in any case, they cannot be systematically implemented.

As @pablosanjose commented it could be really advantageous in these situations to define a new sort of siteselector dependent on the coordination: something like neighbours = not(1). However, since this should act at the level of the Hamiltonian and not when building the lattice object, I propose to simply add this kwarg to unitcell. How does it sound?

vlplot only seems to work for 1D models of for 1D bandstructure cuts. I would expect it to also work for 2D bandstructures (and believe it did at some point in the past?)

It turns out we can do an elegant generalization of our degeneracy resolution machinery to include parametric Hamiltonians, for which a "velocity" along a parameter direction is not defined in terms of Bloch Harmonics. The trick is replace our whole velocity computation code with calls to bloch! where the phases are dual numbers (from the DualNumbers.jl package). This performs automatic differentiation along specific directions. It is fully general (applicable to Bloch phases and parameters) and very, very efficient (surprisingly there seems to be hardly any overhead!). DualNumbers is moreover an extremely light dependence.

This issue discusses the design of GreenFunction, which holds the essential spectral information of a Hamiltonian or an EffectiveHamiltonian (see pablosanjose/Elsa.jl#50 for discussion on the latter). With a GreenFunction we should be able to compute global/local density of states, transport or linear response coefficients of a Hamiltonian (the latter would be two-body GreenFunctions). In more general terms, one can use the GreenFunction of a given Hamiltonian to compute the self-energy induced on another in the presence of a coupling.

CC: @fernandopenaranda, @BacAmorim

Workflow

The workflow would typically be the following. You want to compute some quantity that can be expressed in terms of the Green function of a Hamiltonian h. You would do g = greenfunction(h; projector = ..., solver =...), where projector is a projector on a unit cell subspace. Then g(ω, celldst, cellsrc) would yield the Green function as a matrix on the projected space between the two specified cells. In the above, one could replace h::Hamiltonian by an h::EffectiveHamiltonian describing a coupled system. In that case, the g would be restricted to the intra-Hamiltonian propagator, which I expect greatly simplifies the implementation.

(If we would really want the inter-Hamiltonian propagator from subsystem 1 to 2 in a coupled system, it is rather simple to write as G_12 = G_1 * W * g_2 where G_1 includes the self-energy from 2, but g_2 is the decoupled Green function of 2, and W is the coupling.)

In the following we discuss how we could compute g and G of a Hamiltonian and an EffectiveHamiltonian, respectively.

Hamiltonian

Zero-dimensional Hamiltonian

Let us start with the g::GreenFunction of a zero-dimensional h::Hamiltonian, with matrix m = bloch(h). The goal is to be able to obtain the g(ω) projected matrix (alias for g(ω, celldst, cellsrc) with zero celldst and cellsrc). This can be done in a number of ways. I can think of the following:

• Compute the full inverse , and then project onto the desired sites with projector P: g(ω) = P * inv(ω - m) * P'
• Compute specific columns of the inverse by doing a linear solve: g(ω)[:,i] = P * ((ω - m) \ δ_i) where δ_i is a vector with zero in all elements except a 1 in the i-th. This can be done with e.g. SparseArrays, Pardiso or Mumps
• Compute the full spectrum of m and then do g(ω)=Σ_i Pψᵢψᵢ⁺P'/(ω-ε_i))
• Compute the KPM momenta μ_n of m (they will be matrices in the projected space) and then do g(ω) = ∑_n μ_n G_n(ω) for the appropriate Chebyshev transforms G_n(ω).

The preferred solution in many codes, as far as I know, seems to be to use Pardiso (equus) or Mumps (kwant) to compute specific elements of the inverse. Definitely the worse is computing the full inverse and then projecting. In either case we need to do the full calculation again for any ω (as far as I understand). There is no way to reuse the factorization of ω - m for different ωs.

In contrast, a full diagonalization yields g(ω) for all ω's, but it is impossibly slow for large systems (does not exploit sparsity of the problem). KPM also yields all ω's in one go, and is very efficient, but its accuracy depends on the expansion order.

The ideal would be to support as many methods as possible. We would need to set up different GreenFunctionSolvers, each with its relevant information, that is passed using the solver keyword. The we define

struct GreenFunction{H<:Union{Hamiltonian,EffectiveHamiltonian},S<:GreenFunctionSolver{H}}
    solver::S
end
(s::GreenFunction)(ω) = s.solver(ω)

For example, we could have the trivial solver written as

struct InversionSolver{H<:Hamiltonian{<:Any,0},P<:AbstractMatrix} <: AbstractGreenFunctionSolver
    hamiltonian::H
    projector::P
end
(s::InversionSolver)(ω) = P' * inv(ω*I - bloch(s.hamiltonian)) * P

(Recall that type parameters are Hamiltonian{E,L,T,...}.)

Similarly for other methods, storing in each case whatever data is required to produce g(ω) by calling solver(ω).

Unbounded L-dimensional Hamiltonian

In this case we have less options to draw from, since we need to somehow integrate over momenta or otherwise incorporate the infinite character of the Hamiltonian. As mentioned in Workflow, the Green function projected matrix could be computed with a syntax g(ω, celldst, cellsrc).

I can think of the following methods

• If the Hamiltonian unit cell is small enough to be exactly diagonalized, build its Bandstructure, which is properly split into smoothly connected subbands. Then sum the contribution of each simplex in each subband to yield g(ω, celldst, cellsrc). This is analogous to the full-spectrum approach above. It is unfeasible in the case of very large unit cells, however.
• Use KPM to compute the momenta μ_n projected on pairs of cells up to a cutoff, and keep it in the form of an unflattened L-dimensional tensor in dcell = celldst - cellsrc space or similar. There is a limit to the range of dcell that can be computed from these momenta (given by the expansion order), but likewise there is a limit (set by the momentum-sampling resolution) for which the Bandstructure approach can give reliable results.
• Can we write down a method analogous to the linear-solve approach for periodic Hamiltonians? Maybe writing the linear equation for g(ω, dcell) using Bloch theorem?

In view of the above, we could perhaps concentrate on only two generic solvers to start with: exact-diagonalization (small unit cells) and KPM (arbitrary unit cells)

Semibounded Hamiltonians

This case involves the method of images relative to the boundary. The g should then be in general g(ω, celldst, cellsrc) = g´(ω, celldst, cellsrc) - g´(ω, celldst, cellsrc´) + ..., where g´ are unbounded GreenFunctions and cellsrc´ is the specular image of cellsrc on a given boundary.

EffectiveHamiltonian

As discussed in pablosanjose/Elsa.jl#50, we'll probably end up implementing coupled systems in terms of EffectiveHamiltonians, i.e. Hamiltonians with additional self-energies, computed from the g's of attached hamiltonians. The question then arises as to how one should compute G (coupled Green function) from the EffectiveHamiltonian.

Unbounded EffectiveHamiltonian

The core of the problem is to solve Dyson's equation. The coupled Green function reads, symbolically,

G(ω) = P * (ω - h - Σ(ω))^{-1} * P'

where h should be thought of as a block matrix in cell space, with intracell matrices in the diagonal, and Bloch harmonics in the off-diagonals. Σ(ω) is the total self energy, restricted to a given set of cells, and to sites within the P subspace. P is the site projector.

As the above is the inverse of an infinite matrix in cell space it is unclear how to do this directly. We write an alternative form, using P'(ω-h)^{-1}P = g(ω) and Σ(ω) = PΣ(ω)P':

G(ω) = (1 - g(ω) * Σ(ω))^{-1} * g(ω)

In general, then, a simplistic algorithm to get G(ω) would require that we build the g(ω) * Σ(ω) within the finite set C_Σ of cells with non-zero self-energy, invert the finite-size matrix (1 - g(ω) * Σ(ω)), and multiply by g that propagates from C_Σ to any other cell.

Bounded EffectiveHamiltonian

The case of a bounded (0D) EffectiveHamiltonian can also be solved as described above. However, it allows a second approach based on

G(ω) = P * (ω - h - Σ(ω))^{-1} * P'

since the inverted matrix is now of finite size. It has the added advantage that the matrix is typically sparse. In this case we can also follow the direct solve approach, using SparseArrays, Pardiso or Mumps: G(ω) = P * ((ω - h - Σ(ω)) \ δ_i).

It is not clear a priori which method would be preferable. It all depends on how expensive is your method to compute g(ω) of the bounded Hamiltonian (which may involve a single calculation for all ω), as opposed to doing a factorization and direct solve of ω - h - Σ(ω) (which must be repeated for each ω). The best would be to allow both approaches.

Due to changes in the way Axis and Scenes work in the new GLMakie. PR #155 fixes plots for Hamiltonians. Bandstructures still fail.

I am trying to define the nearest neighbour TB model for graphene, without using any of the presets.

I define a lattice and a model as:

using Quantica

lat = lattice(
    sublat((0.0, 0.0), name = :A), 
    sublat((0.0, 1/sqrt(3)), name = :B);
    bravais = (
        (1/2, sqrt(3)/2),
        (-1/2, sqrt(2)/2)
    )
)

# I know there are better ways to build this
model = 
onsite(1.0; sublats = :A) + 
onsite(2.0; sublats = :B) +
hopping(3.0; sublats = (:B => :A, ), dn = (0, 0)) + 
hopping(4.0; sublats = (:B => :A, ), dn = (1, 0)) + 
hopping(5.0; sublats = (:B => :A, ), dn = (0, 1)) + 
hopping(6.0; sublats = (:A => :B, ), dn = (0, 0)) + 
hopping(7.0; sublats = (:A => :B, ), dn = (-1, 0)) + 
hopping(8.0; sublats = (:A => :B, ), dn = (0, -1))
# the value of the hoppings make no sense.

graphene = hamiltonian(lat, model; 
    orbitals = ((:pz, ), (:pz, ))
)

I get this:

Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 3 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:pz,), (:pz,))
  Element type     : scalar (ComplexF64)
  Onsites          : 2
  Hoppings         : 2
  Coordination     : 1.0

There are 3 harmonics (I was expecting 5). When I inspect the harmonics, I get

> graphene.harmonics[1]
Quantica.HamiltonianHarmonic{2, ComplexF64, SparseMatrixCSC{ComplexF64, Int64}}([0, 0], 
 1.0+0.0im      ⋅    
     ⋅      2.0+0.0im)

> graphene.harmonics[2]
Quantica.HamiltonianHarmonic{2, ComplexF64, SparseMatrixCSC{ComplexF64, Int64}}([0, 1], 
     ⋅      5.0+0.0im
     ⋅          ⋅    )

> graphene.harmonics[3]
Quantica.HamiltonianHarmonic{2, ComplexF64, SparseMatrixCSC{ComplexF64, Int64}}([0, -1], 
     ⋅          ⋅    
 8.0+0.0im      ⋅    )

The hoppings along the directions dn = (1, 0) and dn = (-1, 0) were not included.
Why is that?

Ref: #64

As noted by @BacAmorim, we should perhaps make sure to sort! vertices in each extracted band in the case where it makes sense, i.e. for 1D bands.

It seems the scaling with system size of the build time of a ParametricHamiltonian is quite bad. While allocations are perfectly under control, build time seems to be quite superlinear with the number of sites. This leads to unacceptable overhead when combining them with O(N) methods such as KPM.

(Thanks @fernandopenaranda for the heads up!)

There is a case to be made for an extension of the current ParametricHamiltonian that knows what parameters it depends on. Currently it does not. If hp::ParametricHamiltonian is called as hp(x=2, y=3) it forwards those kwargs to its ElementModifiers, which will simply error if there is any modifier that does not accept all the input kwargs.

In this issue (which started as an email thread between @BacAmorim, @fernandopenaranda and myself) I'd like to discuss the merits of a ParametricHamiltonian that knows about the accepted kwargs for each of its modifiers. These would need to be declared explicitly. Some possible benefits

1- We could call hp(2, 3) in addition to hp(x=2, y=3), with args matched to the list of declared parameters directly

2- We would then be able to easily pass hp to bandstructure, which would interpret each of its declared parameters as "Bloch axes" for its bands.

An alternative to this is to declare a help function hf(z) = hp(x=z, y = z) and allow bandstructure to accept a function in addition to Hamiltonians. A problem of this approach is that if we are in global scope, we would need to declare const hp = parametric(...) as a const explicitly, otherwise hf will be type-unstable

3- Currently

hp = parametric(h, onsite!((o; x) ->...), hopping!((t; y) ->...))

fails when calling it with hp(x=1, y=2), because y is not a kwarg of the onsite!, and x is not a kwarg of the hopping!. A workaround is currently to do

hp = parametric(h, onsite!((o; x, _...) ->...), hopping!((t; y, _...) ->...))

but it is nasty to impose such rule on the user. With compulsory parameter declaration, such as

hp = parametric(h, onsite!((o; x) ->..., parameters = :x), hopping!((t; y) ->..., parameters = :y))

we could avoid the problem.

Some questions:

1- Is the added complexity of parameter declaration worth it? For point 2, is it not clearer to use the hf approach? I am currently not at all sure... The most compelling reason for me is point 3.

2- Can parameter declaration be kept optional? What is clear is that, as long as parameters rely on the kwarg mechanism, they cannot be extracted automatically. We could have them declared optionally as in (3) and apply filtering only if they are declared, so that absent a declaration things work as currently.

3- Can we think of further use cases of this? Perhaps the possibility of creating meshes from a ParametricHamiltonian in which vertices are NamedTuple's instead of SVectors could be useful somehow...

Without this, we need to allocate a new matrix each time we evaluate a ParametricHamiltonian that we want to flatten

Similarly, we need a Matrix(::ParametricHamiltonian) to make it dense

Besides the model approach, it is also possible to specify hamiltonian entries, using the notation h[dn][i, j] = value. What are the integers (i, j) indexing? Sublattice, sublattice+site, or sublattice+site+orbital? Would it make sense to have a notation of the form

h[dn][(sublat1, site1, orb1), (sublat2, site2, orb2)] = value # a scalar

or

h[dn][(sublat1, site1), (sublat2, site2)] = matrix # matrix in orbital space

?

EDIT by @pablosanjose: the actionable change in this issue would be to support h[dn, siteselector(...)] and h[dn, hopselector(...)]. Not sure this is something we need, however, given that we have siteindices(...). Perhaps this could be a motivation to introduce hopindices(...)?

The documentation of hamiltonian has not been updated with #113 (concerning compatible algebraic operations between compatible hamiltonians)

MWE, a Haldane model (only on AA). The sign of the hopping does not alternate (they are all +λ*im), due to a wrong ptrdatum computation

ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> parametric(@hopping!((t, r, dr; λ) -> λ*im*sK(dr); sublats = :A=>:A))

I wonder whether the capability to redefine the onsiteselectors and hoppingselectors of a given Tightbindingmodel made out of several hopping and onsite terms could be useful or not.
I was thinking about a new function:
modifierfunc!(model, region = ...)
where model = onsite_1 + ... + onsite_N+ hopping_1 + ... +hopping_M and onsite_x (hopping_x) are arbitrary onsite(hopping) terms.
This new functionality would enable us to initially write region-independent TBmodels. This may reduce the code complexity when building: a complex heterostructure or different heterostructures out of the same model, since the model is created only once.
For instance, if we want to build the hamiltonians corresponding to two SNS junctions (1 and 2) with different geometries using the same model for the normal(N) and superconducting (S) regions, we could simply do:

hamSNS1 = hamiltonian(lat1, modifierfunc!(modelN,region=(r,dr)->regionN1) + modifierfunc!(modelS,regionS1))
hamSNS2 = hamiltonian(lat2, modifierfunc!(modelN,region=(r,dr)->regionN2) + modifierfunc!(modelS,regionS2))

This is just to say that this exists:
https://github.com/antoine-levitt/wannier A julia code that performs wannierization.

It would be cool to be able to use this in the future to obtain effective few-band models from tight-binding superlattice calculations (I am thinking of obtaining few band models for tblg for example).

Even bether would be to have a full julia workflow: perform a dft calculation with
https://github.com/JuliaMolSim/DFTK.jl Wannierize with
https://github.com/antoine-levitt/wannier and use the tight-binding Hamiltonian with Quantica.jl

As discussed in pablosanjose/Elsa.jl#22 (comment)

Yeah, we could definitely enable such a functionality further on. But to make it really flexible, I think we would need to make a new type for orbital names (maybe OrbitalType) that can store quantum number information. We could then translate any known orbital indentifier like orbitals = :pz into the correcto OrbitalType, and any unknown identifier into OrbitalType with missing quantum numbers. Do open an issue about this so that we don't forget, and if possible describe a use case for the functionality in some more detail to inform the design (thanks a lot for the feedback, by the way!)

It would be a good ideia to create a OrbitalType struct in order to store the name of the orbital an other aditional information, such as, but not limitted to:

• Atomic species
• Angular momentum quantum numbers
• Spin
• Nambu index to describe superconductors

The imediate reason I have to propose this is that in order to model ARPES within a tight-binding model, besides knowing the band structure, it is also necessary to have some information regarding the structure of the Wannier functions in real space.

For concreteness, assuming well localized Wannier orbitals, which can be well described in terms of a single Spherical Harmonic, the ARPES signal is approximatelly give by (for more technical details see https://arxiv.org/abs/1711.02499):

where the amplitude , is given by

In this expression, , are the entries in orbital space, of the eigenvector of the Tight-Binding Hamiltonian (that is, the amplitudes of the Bloch function in the Wannier basis), and is the Fourier Transform of the Wannier state (assuming that it can be separated in radial and angular parts). While inclusion of the radial part is not essential, inclusion of the real spherical harmonic, , is! Failure to include this term, leads to a violation of the rotational symmetries of the original system.

For example, if we have a model of px and py orbitals in a triangular lattice, but do not include the terms when evaluating ARPES, we will obtain Fermi surface maps that do not respect the 6-fold rotational invariance of the system.

Currently, onsite, hopping, @onsite! and @hopping! support the forcehermitian keyword. When true, the Hamiltonian and ParametricHamiltonian constructors try their best to implement the provided model terms so that the result is Hermitian. This is often not simply a matter of adding the adjoint of a term. For example, for onsite we need to add one half the term and its adjoint. For hopping there is not a one half, but when other keywords such as sublats or dn are unspecified, we need to be careful not to add hoppings twice. The code is quite confusing as a result, and bugs creep in easily.

One such bug arose recently upon merging #36 (it is currently present in master). A @hopping! modifier cumulatively applied its adjoint, resulting in a wrong modification. The reason for the bug is that the pointers stored for fast updates in a ParametricHamiltonian could be repeated if forcehermitian=true (once for an element, another for its adjoint) and the cumulative application of both (in an undefined order) could lead to wrong and non-commutative outcomes.

After thinking about this, I think the best solution is not to try to guess what the user wants when writing a model. onsite(2im) should really add 2im to all onsite energies (currently it adds zero, because it symmetrizes by default). Similarly, hopping(2im, sublats = :A => :B) should not add the adjoint hopping from :B to :A. If the user wants that, she should be explicit about it, and do something like hop = hopping(2im, sublats = :A => :B, addadjoint=true), where addadjoint=false by default. (The addadjoint should not be available in onsite, @onsite! or @hopping!, only in hopping.) hop should then resolve into two hopping terms, hop == h+h', where h=hopping(2im, sublats = :A => :B), and h' = adjoint(h) = hopping(-2im, sublats = :B => :A, addadjoint=true)

Doing things this way (never adding adjoints unless the user asks) can produce a little more verbose code, but it would result in far less surprising behavior and much more explicit interface. It would involve introducing adjoint(::TightBindingModel), and would remove the code complexity associated to forcehermitian.

Unfortunately, the change would be very breaking for existing models.

Example:
A user writing

hop = hopping(im, dn = (1,0))
model = hop + hop'

would get dn = (1,0) and dn=(-1,0) in model (currently, you get that already with hop because of the forcehermitian=true default). However,

hop = hopping(im)
model = hop + hop'

would give a model with all zeros (im + im' = 0), because hop and hop' apply to the same links (they are unconstrained). Similarly

hop = hopping(1, range = 1)
model = hop + hop'

would have all hoppings within range equal to 2 (all hoppings pairs r1 => r2 and r2=>r1 have the same range). I think this outcome is clear and unsurprising using this explicit definition.

It might be necessary to benchmark such a change, as applying two terms instead of a symmetrized one may (or may not) be a bit slower.