add OrbitalSliceMatrix

OrbitalSliceMatrix draft

tests pass?

more tests

update docs

no doctest for OrbitalSliceArrays

docs/src/tutorial/greenfunctions.md

@@ -152,7 +152,7 @@ We first define a single lead `Greenfunction` and the central Hamiltonian
 ```julia
 julia> glead = LP.square() |> onsite(4) - hopping(1) |> supercell((1, 0), region = r -> abs(r[2]) <= 5/2) |> greenfunction(GS.Schur(boundary = 0));
 
-julia> hcentral = LP.square() |> onsite(4) - hopping(1) |> supercell(region = RP.circle(10) | RP.rectangle((22, 5)) | RP.rectangle((5, 22)))
+julia> hcentral = LP.square() |> onsite(4) - hopping(1) |> supercell(region = RP.circle(10) | RP.rectangle((22, 5)) | RP.rectangle((5, 22)));
 ```
 The two rectangles overlayed on the circle above create the stubs where the leads will be attached:
 ```@raw html
@@ -199,7 +199,7 @@ Note that since we did not specify the `solver` in `greenfunction`, the `L=0` de
 As explained above, a `g::GreenFunction` represents a Green function of an `OpenHamiltonian` (i.e. `AbstractHamiltonian` with zero or more self-energies), but it does so for any energy `ω` or lattice sites.
     - To specify `ω` (plus any parameters `params` in the underlying `AbstractHamiltonian`) we use the syntax `g(ω; params...)`, which yields an `gω::GreenSolution`
     - To specify source (`sⱼ`) and drain (`sᵢ`) sites we use the syntax `g[sᵢ, sⱼ]` or `g[sᵢ] == g[sᵢ,sᵢ]`, which yields a `gs::GreenSlice`. `sᵢ` and `sⱼ` can be `SiteSelectors(; sites...)`, or an integer denoting a specific contact (i.e. sites with an attached self-energy) or `:` denoting all contacts merged together.
-    - If we specify both of the above we get the Green function between the orbitals of the specified sites at the specified energy, in the form of a `Matrix`
+    - If we specify both of the above we get the Green function between the orbitals of the specified sites at the specified energy, in the form of an `OrbitalSliceMatrix`, which is a special `AbstractMatrix` that knows about the orbitals in the lattice over which its matrix elements are defined.
 
 Let us see this in action using the example from the previous section
 ```julia
@@ -210,17 +210,17 @@ julia> g(0.2)
 GreenSolution{Float64,2,0}: Green function at arbitrary positions, but at a fixed energy
 
 julia> g(0.2)[1, 3]
-5×5 Matrix{ComplexF64}:
- -0.370342-0.0778282im   0.0575525-0.211484im   0.0245456-0.129385im     0.174425-0.155446im       0.100593+0.0134301im
- 0.0575525-0.211484im   -0.0619157+0.0480224im   0.156603+0.256013im    -0.342883+0.0760708im    -0.0414971+0.0510385im
- 0.0245456-0.129385im     0.156603+0.256013im    -0.13008-0.156987im     0.129202-0.139979im       0.155843-0.0597696im
-  0.174425-0.155446im    -0.342883+0.0760708im   0.129202-0.139979im   -0.0515859+0.000612582im   0.0298279+0.109486im
-  0.100593+0.0134301im  -0.0414971+0.0510385im   0.155843-0.0597696im   0.0298279+0.109486im     0.00445114+0.0242172im
-
-  julia> g(0.2)[siteselector(region = RP.circle(1, (0.5, 0))), 3]
-2×5 Matrix{ComplexF64}:
- -0.0051739-0.0122979im  0.258992+0.388052im   0.01413-0.192581im  0.258992+0.388052im   -0.0051739-0.0122979im
-   0.265667+0.296249im   0.171343-0.022414im  0.285251+0.348008im  0.171247+0.0229456im   0.0532086+0.24404im
+5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}:
+ -2.56906+0.000123273im  -4.28767+0.00020578im   -4.88512+0.000234514im  -4.28534+0.00020578im    -2.5664+0.000123273im
+ -4.28767+0.00020578im   -7.15613+0.00034351im   -8.15346+0.000391475im  -7.15257+0.00034351im    -4.2836+0.000205781im
+ -4.88512+0.000234514im  -8.15346+0.000391475im  -9.29002+0.000446138im  -8.14982+0.000391476im  -4.88095+0.000234514im
+ -4.28534+0.00020578im   -7.15257+0.00034351im   -8.14982+0.000391476im  -7.14974+0.000343511im  -4.28211+0.000205781im
+  -2.5664+0.000123273im   -4.2836+0.000205781im  -4.88095+0.000234514im  -4.28211+0.000205781im  -2.56469+0.000123273im
+
+julia> g(0.2)[siteselector(region = RP.circle(1, (0.5, 0))), 3]
+2×5 OrbitalSliceMatrix{Matrix{ComplexF64}}:
+ 0.0749214+3.15744e-8im   0.124325+5.27948e-8im   0.141366+6.01987e-8im   0.124325+5.27948e-8im  0.0749214+3.15744e-8im
+ -0.374862+2.15287e-5im  -0.625946+3.5938e-5im   -0.712983+4.09561e-5im  -0.624747+3.59379e-5im   -0.37348+2.15285e-5im
 ```
 
 ### Diagonal slices
@@ -242,21 +242,25 @@ GreenFunction{Float64,2,2}: Green function of a Hamiltonian{Float64,2,2}
     Coordination     : 3.0
 
 julia> g(0.5)[diagonal(cells = (0, 0))]
-4-element Vector{ComplexF64}:
- -0.34973634684887517 - 0.3118358260293383im
-  -0.3497363468428337 - 0.3118358260293383im
-   -0.349736346839396 - 0.31183582602933824im
- -0.34973634684543714 - 0.3118358260293383im
+4-element OrbitalSliceVector{Vector{ComplexF64}}:
+  -0.3497363468480622 - 0.3118358260294266im
+  -0.3497363468271048 - 0.31183582602942655im
+  -0.3497363468402952 - 0.31183582602942667im
+ -0.34973634686125243 - 0.3118358260294267im
 ```
-Note that we get a vector, which is equal to the diagonal `diag(g(0.5)[cells = (0, 0)])`. Like the `g` Matrix, this vector is resolved in orbitals, of which there are two per site and four per unit cell in this case. Using `diagonal(sᵢ; kernel = K)` we can collect all the orbitals of different sites, and compute `tr(g[site, site] * K)` for a given matrix `K`. This is useful to obtain spectral densities. In the above example, and interpreting the two orbitals per site as the electron spin, we could obtain the spin density along the `x` axis, say, using `σx = SA[0 1; 1 0]` as `kernel`,
+
+Note that we get an `OrbitalSliceVector`, which is equal to the diagonal `diag(g(0.5)[cells = (0, 0)])`. Like the `g` `OrbitalSliceMatrix`, this vector is resolved in orbitals, of which there are two per site and four per unit cell in this case. Using `diagonal(sᵢ; kernel = K)` we can collect all the orbitals of different sites, and compute `tr(g[site, site] * K)` for a given matrix `K`. This is useful to obtain spectral densities. In the above example, and interpreting the two orbitals per site as the electron spin, we could obtain the spin density along the `x` axis, say, using `σx = SA[0 1; 1 0]` as `kernel`,
 ```julia
 julia> g(0.5)[diagonal(cells = (0, 0), kernel = SA[0 1; 1 0])]
-2-element Vector{ComplexF64}:
- -1.1268039540527714e-11 - 2.3843717644870095e-17im
-   1.126802874880133e-11 + 1.9120152589671175e-17im
+2-element OrbitalSliceVector{Vector{ComplexF64}}:
+   4.26186044627701e-12 - 2.2846013280115095e-17im
+ -4.261861877528737e-12 + 1.9177925470610777e-17im
 ```
 which is zero in this spin-degenerate case
 
+!!! tip "Slicing `OrbitalSliceArray`s"
+    An `v::OrbitalSliceVector` and `m::OrbitalSliceMatrix` are both `a::OrbitalSliceArray`, and wrap conventional arrays, with e.g. conventional `axes`. They also provide, however, `orbaxes(a)`, which are a tuple of `OrbitalSliceGrouped`. These are `LatticeSlice`s that represent orbitals grouped by sites. They allow an interesting additional functionality. You can index `v[sitelector(...)]` or `m[rowsiteselector, colsiteselector]` to obtain a new `OrbitalSliceArray` of the selected rows and cols. The full machinery of `siteselector` applies. One can also use a lower-level `v[cellsites(cell_index, site_indices)]` to obtain an unwrapped `AbstractArray`, without building new `orbaxes`. See `OrbitalSliceArray` for further details.
+
 ## Visualizing a Green function
 
 We can use `qplot` to visualize a `GreenSolution` in space. Here we define a bounded square lattice with an interesting shape, and attach a model self-energy to the right. Then we compute the Green function from each orbital in the contact to any other site in the lattice, and compute the norm over contact sites. The resulting vector is used as a shader for the color and radius of sites when plotting the system

docs/src/tutorial/observables.md

@@ -5,6 +5,7 @@ We are almost at our destination now. We have defined a `Lattice`, a `Model` for
 Currently, we have the following observables built into Quantica.jl (with more to come in the future)
 
 - `ldos`: computes the local density of states at specific energy and sites
+- `densitymatrix`: computes the density matrix at thermal equilibrium on specific sites.
 - `current`: computes the local current density along specific directions, and at specific energy and sites
 - `transmission`: computes the total transmission between contacts
 - `conductance`: computes the differential conductance `dIᵢ/dVⱼ` between contacts `i` and `j`
@@ -32,18 +33,16 @@ julia> qplot(h, hide = :hops, sitecolor = d, siteradius = d, minmaxsiteradius =
 
 Note that `d[sites...]` produces a vector with the LDOS at sites defined by `siteselector(; sites...)` (`d[]` is the ldos over all sites). We can also define a `kernel` to be traced over orbitals to obtain the spectral density of site-local observables (see `diagonal` slicing in the preceding section).
 
-We can also compute the convolution of the density of states with the Fermi distribution `f(ω)=1/(exp((ω-μ)/kBT) + 1)`, which yields the density matrix in thermal equilibrium, at a given temperature `kBT` and chemical potential `μ`. This is computed with `ρ = densitymatrix(gs, (ωmin, ωmax); kBT = kBT, mu = μ)`. Here `gs = g[sites...]` is a `GreenSlice`, and `(ωmin, ωmax)` are integration bounds (they should span the full bandwidth of the system). Then, `ρ(; params...)` will yield a matrix over the selected `sites` for a set of model `params`.
+## Density matrix
+
+We can also compute the convolution of the density of states with the Fermi distribution `f(ω)=1/(exp((ω-μ)/kBT) + 1)`, which yields the density matrix in thermal equilibrium, at a given temperature `kBT` and chemical potential `μ`. This is computed with `ρ = densitymatrix(gs, (ωmin, ωmax))`. Here `gs = g[sites...]` is a `GreenSlice`, and `(ωmin, ωmax)` are integration bounds (they should span the full bandwidth of the system). Then, `ρ(µ, kBT = 0; params...)` will yield a matrix over the selected `sites` for a set of model `params`.
 ```julia
-julia> ρ = densitymatrix(g[region = RP.circle(1)], (-0.1, 8.1), mu = 4)
-Integrator: Complex-plane integrator
-  Integration path    : (-0.1 + 1.4901161193847656e-8im, 1.95 + 2.050000014901161im, 4.0 + 1.4901161193847656e-8im)
-  Integration options : (atol = 1.0e-7,)
-  Integrand           : GFermi{Float64}
-    kBT               : 0.0
+julia> ρ = densitymatrix(g[region = RP.circle(1)], (-0.1, 8.1))
+DensityMatrix: density matrix on specified sites using solver of type DensityMatrixIntegratorSolver
 
 julia> @time ρ(4)
-  5.436825 seconds (57.99 k allocations: 5.670 GiB, 1.54% gc time)
-5×5 Matrix{ComplexF64}:
+  6.150548 seconds (57.84 k allocations: 5.670 GiB, 1.12% gc time)
+5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}:
           0.5+0.0im          -7.34893e-10-3.94035e-15im  0.204478+1.9366e-14im   -7.34889e-10-1.44892e-15im  -5.70089e-10+5.48867e-15im
  -7.34893e-10+3.94035e-15im           0.5+0.0im          0.200693-2.6646e-14im   -5.70089e-10-1.95251e-15im  -7.34891e-10-2.13804e-15im
      0.204478-1.9366e-14im       0.200693+2.6646e-14im        0.5+0.0im              0.200693+3.55692e-14im      0.204779-4.27255e-14im
@@ -56,14 +55,14 @@ Note that the diagonal is `0.5`, indicating half-filling.
 The default algorithm used here is slow, as it relies on numerical integration in the complex plane. Some GreenSolvers have more efficient implementations. If they exist, they can be accessed by omitting the `(ωmin, ωmax)` argument. For example, using `GS.Spectrum`:
 ```julia
 julia> @time g = h |> greenfunction(GS.Spectrum());
- 38.024189 seconds (2.81 M allocations: 2.929 GiB, 0.33% gc time, 1.86% compilation time)
+ 37.638522 seconds (105 allocations: 2.744 GiB, 0.79% gc time)
 
 julia> ρ = densitymatrix(g[region = RP.circle(1)])
 DensityMatrix: density matrix on specified sites with solver of type DensityMatrixSpectrumSolver
 
 julia> @time ρ(4)
-  0.000723 seconds (8 allocations: 430.531 KiB)
-5×5 Matrix{ComplexF64}:
+  0.001659 seconds (9 allocations: 430.906 KiB)
+5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}:
           0.5+0.0im  -2.21437e-15+0.0im  0.204478+0.0im   2.67668e-15+0.0im   3.49438e-16+0.0im
  -2.21437e-15+0.0im           0.5+0.0im  0.200693+0.0im  -1.40057e-15+0.0im  -2.92995e-15+0.0im
      0.204478+0.0im      0.200693+0.0im       0.5+0.0im      0.200693+0.0im      0.204779+0.0im
@@ -75,14 +74,14 @@ Note, however, that the computation of `g` is much slower in this case, due to t
 ```julia
 julia> @time g = h |> attach(nothing, region = RP.circle(1)) |> greenfunction(GS.KPM(order = 10000, bandrange = (0,8)));
 Computing moments: 100%|█████████████████████████████████████████████████████████████████████████████████| Time: 0:00:01
-  1.704793 seconds (31.04 k allocations: 11.737 MiB)
+  2.065083 seconds (31.29 k allocations: 11.763 MiB)
 
 julia> ρ = densitymatrix(g[1])
 DensityMatrix: density matrix on specified sites with solver of type DensityMatrixKPMSolver
 
 julia> @time ρ(4)
-  0.006521 seconds (2 allocations: 992 bytes)
-5×5 Matrix{ComplexF64}:
+  0.006580 seconds (3 allocations: 1.156 KiB)
+5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}:
          0.5+0.0im  2.15097e-17+0.0im   0.20456+0.0im  2.15097e-17+0.0im   3.9251e-17+0.0im
  2.15097e-17+0.0im          0.5+0.0im  0.200631+0.0im  1.05873e-16+0.0im  1.70531e-18+0.0im
      0.20456+0.0im     0.200631+0.0im       0.5+0.0im     0.200631+0.0im      0.20482+0.0im

ext/QuanticaMakieExt/plotlattice.jl

@@ -131,7 +131,7 @@ function _hoppingprimitives(ls::LatticeSlice{<:Any,E}, h, opts, siteradii) where
     hp = HoppingPrimitives{E}()
     hp´ = HoppingPrimitives{E}()
     lat = parent(ls)
-    sdict = sitedict(ls)
+    sdict = Quantica.siteindexdict(ls)
     for (j´, cs) in enumerate(Quantica.cellsites(ls))
         nj = Quantica.cell(cs)
         j = Quantica.siteindex(cs)
@@ -144,7 +144,7 @@ function _hoppingprimitives(ls::LatticeSlice{<:Any,E}, h, opts, siteradii) where
             for ptr in nzrange(mat, j)
                 i = rows[ptr]
                 Quantica.isonsite((i, j), dn) && continue
-                i´ = get(sdict, (i, ni), nothing)
+                i´ = get(sdict, Quantica.CellSite(ni, i), nothing)
                 if i´ === nothing
                     opts´ = maybe_getindex.(opts, j´)
                     push_hopprimitive!(hp´, opts´, lat, (i, j), (ni, nj), siteradius, mat[i, j], false)
@@ -158,9 +158,6 @@ function _hoppingprimitives(ls::LatticeSlice{<:Any,E}, h, opts, siteradii) where
     return hp, hp´
 end
 
-sitedict(ls::LatticeSlice) =
-    Dict([(Quantica.siteindex(cs), Quantica.cell(cs)) => j for (j, cs) in enumerate(Quantica.cellsites(ls))])
-
 maybe_evaluate_observable(o::Quantica.IndexableObservable, ls) = o[ls]
 maybe_evaluate_observable(x, ls) = x
 

src/Quantica.jl

@@ -29,7 +29,8 @@ export sublat, bravais_matrix, lattice, sites, supercell,
        spectrum, energies, states, bands, subdiv,
        greenfunction, selfenergy, attach, contact, cellsites,
        plotlattice, plotlattice!, plotbands, plotbands!, qplot, qplot!, qplotdefaults,
-       conductance, josephson, ldos, current, transmission, densitymatrix
+       conductance, josephson, ldos, current, transmission, densitymatrix,
+       OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix, orbaxes
 
 export LatticePresets, LP, RegionPresets, RP, HamiltonianPresets, HP
 export EigenSolvers, ES, GreenSolvers, GS

src/apply.jl

@@ -298,10 +298,15 @@ end
 
 apply(c::CellSites{L,Vector{Int}}, ::Lattice{<:Any,<:Any,L}) where {L} = c
 
+apply(c::CellSites{L,Int}, ::Lattice{<:Any,<:Any,L}) where {L} = c
+
 apply(c::CellSites{L,Colon}, l::Lattice{<:Any,<:Any,L}) where {L} =
     CellSites(cell(c), collect(siterange(l)))
 
 apply(c::CellSites{L}, l::Lattice{<:Any,<:Any,L}) where {L} =
     CellSites(cell(c), [i for i in siteindices(c) if i in siterange(l)])
 
+apply(::CellSites{L}, l::Lattice{<:Any,<:Any,L´}) where {L,L´} =
+    argerror("Cell sites must have $(L´)-dimensional cell indices")
+
 #endregion

src/convert.jl

@@ -23,7 +23,10 @@ Sublat{T,E}(s::Sublat, name = s.name) where {T<:AbstractFloat,E} =
     Sublat{T,E}([sanitize_SVector(SVector{E,T}, site) for site in sites(s)], name)
 
 CellSites{L,V}(c::CellSites) where {L,V} =
-    CellSites{L,V}(convert(SVector{L,Int}, cell(c)), convert(V, siteindices(c)))
+    CellIndices(convert(SVector{L,Int}, cell(c)), convert(V, siteindices(c)), SiteLike())
+
+CellSites{L,V}(c::CellSite) where {L,V<:Vector} =
+    CellIndices(convert(SVector{L,Int}, cell(c)), convert(V, [siteindices(c)]), SiteLike())
 
 function Hamiltonian{E}(h::Hamiltonian) where {E}
     lat = lattice(h)