Inconsistent API between orbitals and sublats

[pablosanjose]

As discussed in #21, that PR creates an inconsistency. It allows us to do

hamiltonian(lat, onsite(@SMatrix[1 0; 0 -1], sublats = 1), orbitals = (:A => Val(2), :B => Val(1)))

This creates a slight problem: what is sublat = 1? Sublattice names are assigned when defining lat, but are referenced only later, when calling hamiltonian. The old API had a milder form of this "nonlocality" issue, as onsite referred to sublattices by index.... With the new API the issue is more severe, as one needs also to remember that sublat 1 in lat is called :A (if one uses the named form of orbital as above).

Also, the following currently errors, because sublat expects an integer or a tuple of integers (sublat indices)

hamiltonian(lat, onsite(@SMatrix[1 0; 0 -1], sublats = :A), orbitals = (:A => Val(2), :B => Val(1)))

A possible solution:
If no names are given to sublattices, they are currently given unique names in lexicographic order as :A, :B, :C... This makes it easy to remember which is which. We would then need to force the sublats to take only names, as onsite(..., sublats = (:A, :D)), thus imposing consistency with the orbital keyword.
If sublattices are given other (distinct) names, we can assume that the user remembers which is which!

[BacAmorim]

I am not completelly sure if the sublats option should always expect a symbol. I think it would be very desirable to make ELSA friendly to Wannier90. The tight-binding Hamiltonian obtained from Wannier90 is outputed as a table, with the values in each line givining us:

nx ny nz i j real imag

with nx, ny and nz being integers that describe the direction of the hopping; i an j being natural numbers that index orbitals. This would represent a hopping with value real + i imag from orbital i at the origin unit cell to an orbital j located in the unit cell at nx * a1 + ny * a2 + nz * a3.

For example, in this notation the nearest-neighbour pz hamiltonian for graphene would be given by:

0 0 0 1 2 -2.7 0.0
0 0 0 2 1 -2.7 0.0
1 0 0 2 1 -2.7 0.0
0 1 0 2 1 -2.7 0.0
-1 0 0 1 2 -2.7 0.0
0 -1 0 1 2 -2.7 0.0

where orbital 1 would be an A site and orbital 2 would be a B site. Wannier90 uses integers to label orbitals (actually sublattice+orbital is labelled as a single flavour). Would it be easy to use the current API to deal with this kind of input?

[BacAmorim]

Also it would be interesting to give not only names to the orbitals, but to also give information about its angular momentum character (l, m) numbers. This information is necessary in order to apply rotations to the Hamiltonian or to model ARPES.

Some presets could be defined, such that :pz = (1, 0) for example.

[pablosanjose]

A Wannier90 interface is of course extremely desirable, so it should definitely be kept in mind when refining our API

Wannier90 uses integers to label orbitals (actually sublattice+orbital is labelled as a single flavour). Would it be easy to use the current API to deal with this kind of input?

So two orbitals (say :px and :py) would be labeled as different i's, even if they belong to the same atom? 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.

To define hoppings using the current API becomes awkward. One possibility would be to create a different sublat per orbital in the unit cell, and define one hopping term per each line in the Wannier90 output. That is not elegant.

I think the best approach in this case would be to define setindex! methods, as alternatives onsite(...) and hopping(...). So, to create such a model we would first create an empty hamiltonian, like for example
h = hamiltonian(lat)
This lat would have one single sublat, with sites in place of Wannier orbitals, even if some share position in space.

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
(For performance one would probably first make h dense. Currently you do hd = Matrix(h), which converts it from sparse to dense representations. Probably, setindex! will be faster in that case (no reallocations). Then you can always go back to sparse with sparse(h).)

If you agree with this solution, I can go ahead and do both things: force sublat to accept only Symbols, and implement setindex!

EDIT: one would need to have a mechanism to ensure that the i,j from Wannier90 correspond to the same i,j from ELSA (and possibly a way of reordering sites to make them match). In principle the ordering shouldn't matter... until you begin distorting your system in a way that depends on site positions, of course.

[pablosanjose]

Also it would be interesting to give not only names to the orbitals, but to also give information about its angular momentum character (l, m) numbers. This information is necessary in order to apply rotations to the Hamiltonian or to model ARPES.

Some presets could be defined, such that :pz = (1, 0) for example.

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!)

[pablosanjose]

Closed in #21 (force sublats to take only lattice names)

If needed, reopen the issue