Wannier90 interface

[BacAmorim]

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

[pablosanjose]

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?

The setindex! proposal is essentially that. You start with no HamiltonianHarmonics and start generating and filling them. Hamiltonian is just a vector of HamiltonianHarmonics, each of which is a sparse matrix over the unit cell and an integer vector indicating the bravais distance between unit cells.

All I know about Wannier90 is thanks to @joselado

Hehe, same as me...

[pablosanjose]

With PR pablosanjose/Elsa.jl#25 one could implement a Wannier90 interface as follows:

1- Obtain orbital positions from Wannier (actually optional, we can place all at the origin, unless we need to do something that depends on real space)
2- lat = lattice(sublat(positions))
3 - ham = hamiltonian(lat) to create an empty hamiltonian, with a single orbital per site
4 - foreach(w->(push!(ham, w.dn); ham[w.dn...][w.i, w.j] = w.v), wannier_list), where wannier_list = [(i = 1, j = 2, dn = (0,0), v = -2.7), ...] is a Vector{<:NamedTuple} with the Wannier90 output.

This would be a way to do it one by one. We could also be more efficient using broadcasting if we first collect all entries with the same dn, and do ham[w.dn...][w.is, w.js] .= w.vs, where is, js are vectors of indices and vs is a matrix of values.

[pablosanjose]

Not so fast, I think we're not done here.

The above approach is ok, and is rather transparent, but it has at least one drawback: it does not participate in the ParametricHamiltonian machinery. Namely, it becomes cumbersome to build a Hamiltonian that depends on parameters.

An alternative approach would be the following

hamiltonian(lat, onsite(1) - hopping(2) + matrixmodel(dns, is, js, vs))

where matrixmodel creates a new type of model (in addition to OnsiteTerm and HoppingTerm) that implements hoppings and onsite energies explicitly, like Wannier90. dns is a vector of dn::NTuple{L,Int}, is and js are vectors of integer orbital indices, and vs is a vector of matrix elements. This is similar to how sparse(I,J,V) works. The input to matrixmodel are basically the columns of the Wannier90 output files.

The good thing is that then, we can do parametric hamiltonians like

hamiltonian(lat, (; multiplier = 1) -> onsite(1) - hopping(2) + matrixmodel(dns, is, js, multiplier .* vs))

Of course this approach would be an alternative to the setindex! functionality already in place.

[BacAmorim]

As an alternative, wouldn't it be possible to create a am empty TightBindingModel and then incrementally add new terms, using the hopping algebra.

For example:

model = tbmodel() # creates a model with no terms, I believe this doesn't exist yet
for (i, j, dn, v) in wannier_list
    model += hopping(v, dn = dn, sublats = i => j)
end

Or is this silly?

[pablosanjose]

In its current implementation that's probably going to be problematic, because TighbindingModel stores its terms in a tuple, as opposed to a Vector{Any}, for performance reasons (so that the type of each term is inferable). The downside is that you cannot efficiently have thousands of terms, the Julia compiler just chokes. There is the possibility to switch to a heap-allocated vector of terms, and have function barriers deal with the type uncertainties (Julia has become really good in this regard, lately), or we could alternatively have a single hopping term encode many terms from wannier_list. A HopSelector is general enough to store any number of specifications efficiently (it specializes on the inputs to dn, sublats, indices, etc, hence the tuple issue mentioned above).

The nice thing is that using that approach we wouldn't need to build a whole new builder type, just a method like hopping(wannier_list) or something. The downside is that Hamiltonian construction is much less optimal than it could. The way Hamiltonian building works for hoppings (onsites are simpler) is essentially by looping over sites, finding all potentially connected sites at a given range distance, and then query the hopselector to see if they should be selected (i.e. if they are included in sublats, in indices, etc). So you see that, while the approach is completely general, it is O(N^2) in the number of sites if we don't specify a range (if we have a range we use KDTrees, which are much more efficient). If we know the hopping pairs beforehand (in wannier_list), it is rather silly to do things this way.

[BacAmorim]

because TighbindingModel stores its terms in a tuple, as opposed to a Vector{Any}, for performance reasons (so that the type of each term is inferable)

I see. I thought that the addition of model terms, just just sintatic sugar for push!. It is deeper than that.

The downside is that Hamiltonian construction is much less optimal than it could.

So even if a HopSelector as defined dn and indices, it still performs a range search? Just to make sure:i => j indices and dn completely specify a hopping correct?

So what is the approach for this? Define a WannierHoppingstruct that stores hamiltonians in the w90 format:

struct WannierHopping{Ti, Tv, D}
    scr::Ti # source index
    dst::Ti # destinty index
    dn::SVector{D, Ti}
    val::Tv
end

[Or maybe use a WannierModel that is a struct of vectors (not sure what is better)]
And then add a method hopping(lat, ::Vector{WannierHopping}) ?

There is still the issue of whether to encode or not sublattice and orbital information. In a previous discussion, I think you mentioned is was better to associate one Quantica index with one orbital per site to each Wannier orbital. Do you still think this is the best approach? Surelly, it is the simplest. Perhaps in addition to that behaviour, we could have Quantica pack the hoppings into orbital subspaces, provided a OrbitalStructure is passed as an argument.

Having this Vector{WannierHopping}) / WannierModel would probably make #172 (and even the whole setindex! apparatus) unnecessary? If we want to specify the hoppings by indices and dn use a WannierModel/Hopping, otherwise use the magic of hopping.

What do you think?

[pablosanjose]

So even if a HopSelector as defined dn and indices, it still performs a range search?

It does, but it is a missed optimization opportunity, currently. It would be very easy to specialize a function called Quantica.targets to iterate over pre-specified indices, and thus avoid the O(N^2) overhead

Just to make sure:i => j indices and dn completely specify a hopping correct?

That is correct

So what is the approach for this? Define a WannierHoppingstruct that stores hamiltonians in the w90 format:

The internal applyterm! routine that processes HoppingTerms and OnsiteTerms currently works by adding elements to (i, j, v) arrays, that later, after all terms are processed, are passed to the sparse constructor. So indeed, what we need is just to define a new type of WannierHopping <: TightbindingTerm with its own applyterm! method that converts whatever it carries directly into i,j,v entries, without going thought the whole siteselector/hopselector machinery.

So, that new applyterm!(builder::IJVBuilder, term::WannierTerm) can be very trivial if WannierTerm can produce the i,j,v terms directly, i.e. if we can have, as you propose, something like

struct WannierHoppings{Ti, Tv, D} <: TightbindingTerm
    scrs::Vector{Ti} # source indices
    dsts::Vector{Ti} # destinty indices
    dn::SVector{D, Ti}
    vals::Vector{Tv}
end

The catch is that the indices chosen by Wannier90 should perfectly match the ones in our lattice. So the lattice itself should be built using the same file input. The nice thing is that, as Rui was explaining to me the other day, Wannier90 also gives you the expectation value of position between any two Wannier orbitals with a hopping. We could then define the following generic term, that includes onsites and hoppings, and imports from the file also the position expectation values for each pair

struct WannierTerm{T,E,L} <: TightbindingTerm
    dn::SVector{L, Ti}
    scrs::Vector{Int} # source indices
    dsts::Vector{Int} # destinty indices
    vals::Vector{Tv}
    positions::Vector{SVector{T,E}}
end

struct WannierModel{T,E,L}
    terms::Vector{WannierTerm{T,E,L}}
    bravais::Bravais{T,E,L}
end

where WannierModel is distinct from TighbindingModel (which is the current container of HoppingTerms and OnsiteTerms). We can then define a hamiltonian(model::WannierModel), a special method that does not take a lattice as input, but that builds it using the Wannier model itself.

We could then consider if it would be worthwhile to allow combing TightbindingModels and WannierModels. Neither is really a superset of the other, so that is tricky. We could discuss if it would make sense at all. If we do I think the best generalization would be to make WannierModel a strictly more general model type, and do

struct WannierModel{T,E,L,N,TB<:Tuple{Vararg{TightbindingModelTerm,N}}}
    wterms::Vector{WannierTerm{T,E,L}}
    tbterms::TB
    bravais::Bravais{T,E,L}
end

and then define term addition to push into the relevant field, depending of the type of the added term. The WannierModel type retains the ability to define a lattice at the same time as the Hamiltonian, and then, when we're done adding i,j,v terms from the wterms we continue with the tbterms. What do you think?

[pablosanjose]

Uhm, a question that I just had: is the Wannier90 output always giving you a scalar hopping amplitude per orbital pair? How does it deal with degeneracies, e.g. spin degeneracies? Do we ever get vals to be a vector of inter-orbital blocks, in analogy to Quantica's SMatrix Hamiltonian blocks for multiorbital sites?

EDIT: I expect we don't, so that if there are "degeneracies" (in the sense of several Wannier orbitals with the same spatial position) we should either (1) accept having several single-orbital sites at the same point or (2) detect there are degeneracies (by looking at positions) and then bundle together these wannier orbitals as belonging to the same Quantica site. That would require some preprocessing of the file, but it might be relatively easy to do robustly (though it would not be type-stable - vals type would depend on the file content - but that's not a real issue)

[BacAmorim]

is the Wannier90 output always giving you a scalar hopping amplitude per orbital pair?
Yes, the hopping between two orbitals is always a scalar. The output, is a text document where each line has information of one hopping/onsite energy, in the form:

n1   n2   n3   m   n   tr   ti 

Which encodes the hopping:

with and m, n are orbital indices. tr and ti are the real and imaginary parts of the hopping which are always scalar quantities.

In Wannier90, two orbitals that share the same center are treated in the same way as two orbitals with different centers. In general, the fact that two wannier orbitals share the same center is an accident, as the center of the orbital is also a parameter that can be varied in the minimization of the spread.

EDIT: I expect we don't, so that if there are "degeneracies" (in the sense of several Wannier orbitals with the same spatial position) we should either (1) accept having several single-orbital sites at the same point or (2) detect there are degeneracies (by looking at positions) and then bundle together these wannier orbitals as belonging to the same Quantica site. That would require some preprocessing of the file, but it might be relatively easy to do robustly (though it would not be type-stable - vals type would depend on the file content - but that's not a real issue).

Probably we should support both with an option group_orbitals = true/false.

[lxvm]

I am curious about how this package encodes tight-binding Hamiltonians and I came across this thread. A recent development for Wannier90 in Julia is https://github.com/qiaojunfeng/Wannier.jl. Hope this helps!

[pablosanjose]

That's a great pointer @lxvm, thanks! We will certainly look into doing an extension into Wannier.jl instead of a custom solution when we activate work on this issue.