Add sublattice arguments to onsite and hopping functions

[fernandopenaranda]

I think that it could be a good idea to generalise the onsite and hopping functions to allow definitions that optionally take sublattices as arguments too. For example, something like: hopping((r, dr, s1, s2) -> ...)

[pablosanjose]

Nice suggestion @fernandopenaranda

I've tried out a quick implementation in #43. I'm not fully convinced yet whether it can be made cleaner (internally), but at least there is very little perf regression.

[pablosanjose]

That said, I am not sure we want to encourage writing models in this way (á la MathQ). It tends to litter your models with conditionals (on the sublattice name) instead of splitting your model into terms acting on sublattices, which I think is cleaner. The current API doesn't need to be more verbose if one builds models by pieces. For example a Haldane model could be written as

hop(s) = hopping((r, dr) -> im * cos(3*atan(dr[2], dr[1])), sublats = s)
model = hopping(1, sublats = (:A,:B)) + 0.1 * (hop(:A) - hop(:B))
h = LatticePresets.honeycomb() |> hamiltonian(model)

This is actually clearer (to me) and less verbose than using conditionals on the sublattices:

hop = hopping((r, dr, s1, s2) -> ifelse(s1 == s2, ifelse(s1 == :A, 1, -1), 0) * im * cos(3*atan(dr[2], dr[1])))
model = hopping(1, sublats = (:A,:B)) + 0.1 * hop)
h = LatticePresets.honeycomb() |> hamiltonian(model)

[fernandopenaranda]

Yes, I agree. As you say, verbose conditionals are something to avoid. The first code is not redundant so it makes perfect sense not to merge this PR.

[BacAmorim]

One question, and this might be off-topic but, is it currently possible to define hoppings with orbital arguments?

The use case I have in mind is the definition of Slater-Koster models.

Suppose we have a s, px, py model in a honeycomb lattice. Is there a way to define a hopping from sublattice :A and orbital group (:px, :py) into sublattice :B orbital (:s)? Or would this just be written as something like:

hopping(@SMatrix [0 1 0; 0 0 0; 0 0 0], sublats=(:A=>:B))

?

[pablosanjose]

@BacAmorim , yes, it's always been possible (EDIT: current API established in #21). You just need to pass a rectangular matrix that matches the orbitals you declare

julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(@SMatrix[1; 0], sublats = (:A, :B)), orbitals = (:A => (:px, :py), :B => :s))
Hamiltonian{<:Lattice} : Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 5 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Orbitals         : ((:px, :py), (:s,))
  Element type     : 2 × 2 blocks (Complex{Float64})
  Onsites          : 0
  Hoppings         : 6
  Coordination     : 3.0

Internally, however, sparse Hamiltonians have a well defined eltype for performance. In this case it is a square 2x2 @SMatrix, padded with zeros. In general it is a square SMatrix that is big enough to hold the hoppings between any pair of sublattices. That is an implementation detail, however, that is meant to be mostly invisible to the user. If you do flatten(h), the eltype will be ComplexF64 scalars, with no padding zeros.

EDIT: You can also use the sublat = :A => B syntax, but note that that corresponds to sublat = (:B,:A), as the latter refers to row-column blocks.

julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(@SMatrix[1 0], sublats = :A =>:B), orbitals = (:A => (:px, :py), :B => :s))