Taking Operators seriously

[pablosanjose]

We currently have Operator (an AbstractHamiltonian wrapper, created e.g. with current(h; direction = ...), with scalar eltype), BarebonesOperator (a container of sparse Bloch harmonics of arbitrary eltype) and harmonics::Vector{<:HybridSparseMatrix} (inside Hamiltonian).

This is screaming for unification. We should create a single Operator type that holds a Lattice and any type of AbstractMatrix or HybridSparseMatrix (should the latter be an AbstractSparseMatrix also?). Hamiltonian is then an OrbitalOperator <: Operator, with OrbtialElementType. current(h) and position(w) should be VectorSparseOperator <: Operators, with SVector type.

supercell and many other functions should be simplified to operate on any operator, unless their eltype makes it impossible.

This is a rather major refactor, but would simplify things considerably, and would allow us to place AbstractHamiltonian's and other Operators on the same footing. Perhaps a post-v1.1 undertaking.

[pablosanjose]

What are the features we need out of Operator?

• Conceptually, it should encode a block-diagonal matrix on a periodic lattice of cells. It should perhaps be named LatticeOperator
• It should be generic, storing a general object T (usually some subtype of AbstractArray) for each cell distance dn::SVector{L,Int}. It should then be LatticeOperator{T,L}
• It should allow indexing o[dn] to give the actual T for that dn
• It should also allow view(o, dn)
• It should allow iteration like Dict or Dictionary, yielding pairs dn => view(o, dn)
• It should allow insert!, delete! and set!, just like Dictionary. Should it then probably be implemented internally as a Dictionary?
• It should support doing bloch sums in place, using e.g. sum!(f::Function, m::T, o::LatticeOperator{T}), where f(dn) = cis(-dot(dn, ϕs)). Functionally this should just do for (dn, m´) in o; m .+= f(dn) .* m´; end; m.
• It should support mul!, so that mul!(o, o1, o2, α, β; grow = true) updates o::LatticeOperator, so that o[dn] = α Σ_k o1[k] * o2[dn-k] + βo[dn] . If grow, new entries dn should be generated in o from the ones in o1 and o2. Otherwise, the number of entries in o should not change and the operation should be allocation free, as long as mul! over Ts is allocation free
• It should take advantage of multithreading
• Optional: It may have a fastpath for mul! based on FFT, as it is essentially an L-dimensional convolution?

[pablosanjose]

We note that an equilibrium density matrix will eventually be a LatticeOperator{Matrix{Complex{T}}, L}. This will be used e.g. in real-time equations of motion. For this application it is crucial that mul! is made as efficient as possible. This could inform about the best internal representation. A Dictionary from Dictionaries.jl? With what key order?

Also, we must consider if a state, or ket, should be a special case of operator. Conceptually, they are not, but functionally, they behave very similarly.