Provide wrappers to matrices / vectors to indicate parent Hamiltonian an allow use of siteselector

[BacAmorim]

It would be really useful if one could use the siteselector and LatticeSlice machinery of Quantica with user provided matrices and vectors.

Example:

Assume we have a periodic system with Quantica Hamiltonian h. We want to evaluate the Green's function in momentum space for a given k point as:

gk = inv(omega*I - h(k))

Here gk will be a simple matrix. If would be useful if we could say that gk exists in the same space as h, such that we could use the indexing machinery that works for h. Maybe this could de achieved by doing something like:

struct Operator{H, O}
    parent::H
    operator::O
end

function operator(h, matrix)
    
    op = Hamiltonian(
        h.lattice, 
        h.blockstruct, 
        [Harmonic((0, ), matrix)],
        matrix
        )
        
    return Operator(h, op)
    
end

This would require some promotion of matrix to a HybridSparseMatrix, such that then it can be internatally interpreted as a Quantica Hamiltonian, thus allowing for the use of indexing.

I would also be nice if the same logic could be extended to states or selected columns of an operator. For example if psi is an eigenstate of h, then by defining someting like

struct StateVector{H, V}
    parent::H
    vector::V
end

then calling state(h, psi) would create a wrapper of psi, promising that psi lives in the same Hilbert space as h, allowing for the use use of sitelector via psi[site...]. The same should also work for collections of states / operator columns. I have no idea how this could be implemented.

[pablosanjose]

This issue points to a core problem in current master: we have no concept for an AbstractMatrix defined on the orbitals of a unit cell or, more generally, on a LatticeSlice. That is the subject of #230 too.

We do have AbstractHamiltonian that wraps a number of Harmonics. Each of these are special types of sparse matrices that know about the unit cell and orbital structure. That allows us to update parameters and compute Bloch matrices. But this information is completely lost when we do h(ϕs; params...) or h[dn] (we get a simple SparseMatrixCSC). Similarly, when we get a GreenFunction matrix, we "forget" what sites or orbitals it spans.

As a result, as you point out @BacAmorim, we cannot do siteselector indexing of the matrices produced by Quantica.

So, to summarize the problem:

1. We need a way to make use of matrices produced by Quantica using siteselectors instead of indices, e.g. the output of h(φs; params...) and h[dn] should be indexable with siteselectors, but also of g(ω; params...)[siteselector´, siteselector] (g is a GreenFunction). These objects are currently AbstractMatrices over orbitals living on a certain LatticeSlice, but the matrices know nothing of the underlying orbital structure or LatticeSlice.
2. In addition to siteselector indexing, it would be desirable to also be able to index these matrices with single-site identifiers (not just site selectors). By single-site identifiers I mean cellsite(dn::SVector{L,Int}, i::Int) that identify a single site in a lattice.

To address these two needs we need a new concept: a matrix defined over the orbitals in a LatticeSlice that knows about its domain, and the OrbitalStructure defined on it. The Operator struct you propose wraps a Hamiltonian, but we only need a LatticeSlice and an OrbitalBlockStructure. We could then define

struct LatticeSliceMatrix{C,M<:AbstractMatrix{C},L<:LatticeSlice,B} <: AbstractMatrix{C}
    parent::M
    rows::L
    cols::L                                 # can be cols === rows
    obstruct::OrbitalBlockStructure{B}      # B here is the block type enconding orbitals
end

This type could implement the full AbtractMatrix interface, and therefore act exactly as the parent matrix for all usual matrix operations, including sparse. But then it would also implement two new types of indexing, using SiteSelector indices and CellSite indices. Indexing a LatticeSliceMatrix with one or two SiteSelectors would again return a LatticeSliceMatrix. In contrast, indexing with two CellSites would return an element or type B.

This indexing behavior is analogous to that of a conventional m::AbstractMatrix, where m[1:3, 2:4] and m[3, 5:6] returns a new matrix, while m[2,3] returns an element. So SiteSelectors act as index containers, and CellSites as scalar indices.

How would that type be produced by different Quantica objects?

• When we do hφ = h(φs; params...), then, we would get a LatticeSliceMatrix{<:SparseMatrixCSC}, and we could probably just use that as if it where a SparseMatrixCSC. We could also unwrap it with parent(hφ) if we wanted the actual sparse matrix. If we do hφ[siteselector] we would get another LatticeSliceMatrix with the selected slice. If we fo view(hφ, siteselector) we would get a SubArray{LatticeSliceMatrix}.
• When we do g(ω; params...)[selector1, selector2] we would get a LatticeSliceMatrix{<:Matrix}. If we indexed this with a selector3 we would then get another LatticeSliceMatrix{<:Matrix}, where any site not included in selector1 and selector2 would be dropped from the wrapped LatticeSlices.
• When we do ρ(; params...) we would get a ρs::DensityMatrixSolution, analogous to g(ω; params...)::GreenSolution. This we could index using siteselectors once more, or with cellsites, so we could do ρs[cellsites((0,0), 1), cellsites((1,1), 2)]. This could then be passed to a HartreeFock model using the approach in Self-consistent mean fields #229

Note: the HybridSparseMatrix in Harmonic probably cannot trivially be made a LatticeSliceMatrix, since it doesn't conform to the AbstractMatrix interface. It has different and more complicated requirements. In particular it doesn't have one well-defined eltype, but rather two. It also has two different indexing and iteration behaviors (one for flat and another for unflat)

You mention also the concept of a AbstractVector, which could be implemented by generalizing LatticeSliceMatrix to a more general LatticeSliceArray, or we could simply define an additional LatticeSliceVector type with a single rows field. But I'm not sure where we would end up needing such an object. What exactly do you have in mind?

[BacAmorim]

OK, your proposal looks like what I had in mind.

Regarding the LatticeSliceArray. My motivation is for density matrices and Green's functions. Some algorithms (like traditional KPM) naturally compute some columns of a hamiltonian function. For each column, each of the entries correspond to a lattice sites, but the each column might not be associated to a particular lattice site, but to some superposition of sites (examplle: in traditional KPM using stotastic trace evaluation the columns correspond to random vectors, or KPM can also be used to compute k-resolved density of states, where the initial vectors correspond to plane waves).

Does this make sense?

[pablosanjose]

Note: this issue strongly overlaps with the long and winding #230 discussion, and it is quite a bit simpler. So if @BacAmorim finds this agreeable, I wouldn't be opposed to closing #230 and abandoning for the moment the caching machinery discussed there, which would no longer be needed to implement #229.

[pablosanjose]

Does this make sense?

I guessed this is what you had in mind. Yes, it makes sense, but wouldn't a one- or few-column LatticeSliceMatrix work too? Do we really need a new Vector type for that?

[pablosanjose]

So what would be a sketch for a general LatticeSliceArray? Something like:

struct LatticeSliceArray{C,N,M<:AbstractArray{C,N},L<:LatticeSlice,B} <: AbstractArray{C,N}
    parent::M
    latslices::NTuple{N,L}
    obstruct::OrbitalBlockStructure{B}      # B here is the block type enconding orbitals
end

So not that bad, we just have to add one more parameter N. Since the implementation of SiteSelector indexing of this would probably need to go through the computation of orbital indices along each axis, there is no cost to doing a general implementation for arbitrary N, I think, so yes, let's do this. It might come in handy eventually as you say.

[BacAmorim]

one- or few-column LatticeSliceMatrix work too? Do we really need a new Vector type for that?

From what I understood in your proposal of LatticeSliceMatrix both the rows and column indices correspond to LatticeSlices. In what I am thinking the row indices correspond to a LatticeSlice but the column indices don´t.

This is ilustrated by computing the eigenvalues of an Hamiltonian. We get a Matrix U where each column, U[:, n] is an eigenvector. So the row indices correspond to a LatticeSlice, but the column indices don´t (they are just the eigenvector/value index). So in this case, it would make sense to call U[siteselector, n] for the n-th eigenvector, but is doesn´t make sense to call `U[siteslector1, siteselector2]'.

Question: could all of this issue be solved by just using U[getindex(siteselector), n] ? I have tried this in the past, but faced the problem that getindex(siteselector) does not actually imediatly return indices that can be used to index a Matrix. The actual indices are hidden somewhere inside the return LatticeSlice. (in that way the name getindex is perhaps a bit misleading).

[pablosanjose]

Ah, that's an excellent point! Fortunately, Julia extremely nice AbstractArray API allows us to do this very cleanly. We need

struct OrbitalAxis{T,E,L,B}
    latslice::LatticeSlice{T,E,L}
    orbstruct::OrbitalBlockStructure{B}
end
struct LatticeSliceArray{C,N,M<:AbstractArray{C,N},A<:NTuple{N,Union{OrbitalAxis,UnitRange}}} <: AbstractArray{C,N}
    parent::M
    axes::A
end

In your StateVector case the axes field would be (l::OrbitalAxis, 1:n). For a Matrix it would be (::OrbitalAxis, ::OrbitalAxis). (EDIT: correction)

To make this work we just need to add a method Base.to_indices(::LatticeSliceMatrix, ::SiteSelector) which yields the actual indices you want, and a Base.to_indices(::LatticeSliceMatrix, ::CellSite) that returns a scalar. Then Julia takes care of the rest, and we have access to consistent indexing rules and to the actual indices of the parent array.

[pablosanjose]

I've been iterating the above design of a LatticeSliceArray a bit more on a piece of paper. I think this is a good opportunity to clean up the internals of our representations of Lattice slices, which is somewhat confusing. I write it here so I don't loose it.

We currently have CellSites{L,V} and CellOrbitals{L,V,R} which are very similar parametric types representing a subset of sites and a subset of orbitals (of generic type V) in a unit cell of indices cell::SVector{L,Int}. When the subset is a Vector{Int} we call them AppliedCellSites and AppliedCellOrbitals (awful names, I know). These variants then get wrapped in LatticeSlice and OrbitalSlice. But internally, each of these pairs are quite different, since the orbital variants need to keep information about the site groupings of orbitals. The R type for example corresponds to a list of orbital groupings for each site in the subset, that can be Vector{UnitRange} or Missing. Additionally, CellOrbitals is part of the general GreenSolver API: they are the basic index that every GreenSolver need to support. The result is a rather messy type structure that I've never been really satisfied with.

EDIT: tweaks of the original post follow

The refactor would be as follows: we will have a single parametric type CellIndices{L,I,G} that can represent any subset of any type in a given cell.

struct CellIndices{L,I,G<:Union{Missing,OrbitalGroups, OrbitalGroupsSkipped}}
    cell::SVector{L,Int}
    inds::I    # can be anything, Int, Colon, Vector{Int}, etc. Represents one or more entities in cell.
    groups::G  # if G<:OrbitalGroups, `inds` are orbitals, and `groups` are their grouping in sites.
end

Then we define a number of aliases:

• CellSites{L,I} = CellIndices{L,I,Missing}
• CellOrbitals{L,I} = CellIndices{L,I,OrbitalGroups}
• CellOrbitalsNoGroups{L,I} = CellIndices{L,I,OrbitalGroupsSkipped}
• Site{T,E,L} = CellSites{L,Pair{Int,SVector{T,E}}}. This replaces the current CellSite. It stores i=>r, i.e. both the site index i and its position r. For use in mean field models and also in the few places where we currently use single CellSites.

In CellOrbitals above we have used orbital groups defined as

const OrbitalGroups = Dictionary{Int,UnitRange{Int}}
struct OrbitalGroupsSkipped end; # singleton type to avoid allocating a Dictionary when we don't need groups

OrbitalGroups are thus a site => inds_range Dictionary, that we can skip allocating if it is not needed.
We also define another Dictionary alias that maps cell => CellIndices

const Subcells{L,C<:CellIndices} = Dictionary{SVector{L,Int},C}

Then, we will define slices with a unique type

struct LatticeSlice{T,E,L,C<:CellIndices}
    lat::Lattice{T,E,L}
    subcells::Subcells{L,C}
end

What we currently call OrbitalSlice{L} is now renamed to

const OrbitalSubcells{L} = Subcells{L,CellOrbitals{L,Vector{Int}}}

which, note, is just a Dictionary{SVector{L,Int},CellOrbitals{L,Vector{Int}}}.
The name OrbitalSlice is now reserved for a special type of LatticeSlice with C<:CellOrbitals encoding a vector of orbital indices,

const OrbitalSlice{T,E,L} = LatticeSlice{T,E,L,CellOrbitals{L,Vector{Int}}}

Note the use of Dictionary (from Dictionary.jl) instead of Vector. Dictionarys are like Dict in that they can be indexed with hashed keys, but can be iterated like a Vector (they are actually just Vectors with Dict-like indexing). This is a long due change that should pay off when searching for unit cells in a LatticeSlice or OrbitalSlice.

After such a refactor, the OrbitalAxis for LatticeSliceArray type becomes unnecessary. We can just define

struct OrbitalSliceArray{C,N,M<:AbstractArray{C,N},A<:NTuple{N,Union{OrbitalSlice,UnitRange}}} <: AbstractArray{C,N}
    parent::M
    axes::A
end

which is conceptually much clearer: the type for a matrix defined over a given OrbitalSlice should be called OrbitalSliceArray, and its axes are simply OrbitalSlices (or UnitRanges for non-slice axes as proposed by @BacAmorim). These OrbitalSlice axes contain enough information to convert a SiteSelector into a bunch of orbital indices for parent, due to the definition of GroupedOrbitals.

The above should be implemented in two steps: one PR with an invisible internal refactor, and a second PR introducing OrbitalSliceArray.