`setindex!` for multisite/multiorbital models

[BacAmorim]

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(...)?

[pablosanjose]

The i, j are the site indices, actually, which go from 1 to N = size(h, 1). These are the numbers you see in the tooltips when you plot a Hamiltonian with using VegaLite; vlplot(h). They are also the numbers you use when you specify, e.g. onsite(3, indices = (2,4,9)). When you do h[dn][i,j] = s::SMatrix you are assigning a block between sites i and j, which includes all orbitals in the SMatrix. If there is only one orbital you can also use a scalar instead of an SMatrix.

So, we never index sites and orbitals using an index within a sublattice. We could offer an interface of this type but it is unclear to me how useful that would be, since the actual indices are rather arbitrary. You actually have to look them up by plotting the Hamiltonian. So in practice what we do most often is talk about site positions when addressing them. We could think about a setindex of the form h[dn][siteselector] .= s and h[dn][hopselector] .= s where siteselector/hopselector are any selector we currently use to define hopping or onsite. These include indices, sublats, regions, etc. I would be in favor of such a change.

There is however a caveat: when setting a value in a sparse matrix, which is the internal representation of the Hamiltonian, you induce a reallocation of the whole matrix in memory. So setting a large number of elements this way is not efficient. That's why Kwant has a Builder that gets finalized, i.e. assembled into a sparse matrix when you're done.

[BacAmorim]

So with your proposal, what I was suggesting, e.g. for a hopping between site 1 of sublattice A to site 2 of sublattice B, would be written as:

h[dn][hopselector(sublats = (:A => :B, ), indices = (1 => 2, ))]

Given that dn can actually be passed to ` hopselector it would perhaps make more sense:

h[hopselector(sublats = (:A => :B, ), indices = (1 => 2, ), dn = dn)]

There is however a caveat: when setting a value in a sparse matrix, which is the internal representation of the Hamiltonian, you induce a reallocation of the whole matrix in memory. So setting a large number of elements this way is not efficient. That's why Kwant has a Builder that gets finalized, i.e. assembled into a sparse matrix when you're done.

Yes, I am aware of that. In my own codes I have stored hoppings as a kind of sorted IJV list (more live IJDnV to include dn vectors). In this way adding a new hopping is just push! a new (i, j, dn, v), but I am not sure if building the bloch hamiltonian for that is the most effective.

I think it would be nice to keep this kind of functionality, but discourage its use as it is not very performant.

[pablosanjose]

Not quite, because indices takes unique site indices, not indices within a sublattice. You would still need to find the indices of the sites you want to connect. Then sublats would not be needed.

This form would rather be desirable when you want to, say, change the hopping between all sites at a certain distance, or the onsite of sites in a certain region. Indexing with site indices is possible but, as I said, requires that you first identify your indices via a plot or something.

Storing hoppings in an IJV form is nice for construction, but not so much for computations. If we see that, for future functionality (perhaps Wannier90 imports or other) a Builder type makes sense, it would be very simple to add, with a method hamiltonian(b::Builder) to "finalize" it.

[BacAmorim]

Not quite, because indices takes unique site indices, not indices within a sublattice. You would still need to find the indices of the sites you want to connect. Then sublats would not be needed.

I see.

requires that you first identify your indices via a plot or something.

This raises another question: assuming we have N sublattices, each with Nn sites (n = 1, ..., N) do indices i = N{n-1} + 1, ..., N_{n-1} + N_n, label site i from sublattice n ?

Storing hoppings in an IJV form is nice for construction, but not so much for computations. If we see that, for future functionality (perhaps Wannier90 imports or other) a Builder type makes sense, it would be very simple to add, with a method hamiltonian(b::Builder) to "finalize" it.

That would probably make the most sense I agree.

[pablosanjose]

This raises another question: assuming we have N sublattices, each with Nn sites (n = 1, ..., N) do indices i = N{n-1} + 1, ..., N_{n-1} + N_n, label site i from sublattice n ?

Yes, that's guaranteed. Labels are ordered by sublattice. Actually, doing Quantica.offsets(ham) gives you the offsets for the site indices in each sublattice (the first is 0 and the N+1 is the total number of sites

[BacAmorim]

Thanks for the explanation!

[pablosanjose]

So, the actionable part of this issue is the implementation of siteselectors and hopselectors with setindex!, right? Could you please add an EDIT to the OP with the proposal so we can track that?

[pablosanjose]

This is effectively closed by #179 .