diff --git a/docs/make.jl b/docs/make.jl
index 487ba253..977cbf72 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -23,7 +23,8 @@ makedocs(;
             "Hamiltonians" => "tutorial/hamiltonians.md",
             "Bandstructures" => "tutorial/bandstructures.md",
             "GreenFunctions" => "tutorial/greenfunctions.md",
-            "Observables" => "tutorial/observables.md"
+            "Observables" => "tutorial/observables.md",
+            "Advanced topics" => "tutorial/advanced.md"
             ],
         "Examples" => "examples.md",
         "API" => "api.md",
diff --git a/docs/src/tutorial/advanced.md b/docs/src/tutorial/advanced.md
new file mode 100644
index 00000000..b69caa7c
--- /dev/null
+++ b/docs/src/tutorial/advanced.md
@@ -0,0 +1,94 @@
+# Non-spatial models and self-consistent mean-field problems
+
+As briefly mentioned when discussing parametric models and modifiers, we have a special syntax that allows models to depend on sites directly, instead of on their spatial location. We call these non-spatial models. A simple example, with an onsite energy proportional to the site **index**
+```julia
+julia> model = @onsite((i; p = 1) --> ind(i) * p)
+ParametricModel: model with 1 term
+  ParametricOnsiteTerm{ParametricFunction{1}}
+    Region            : any
+    Sublattices       : any
+    Cells             : any
+    Coefficient       : 1
+    Argument type     : non-spatial
+    Parameters        : [:p]
+```
+or a modifier that changes a hopping between different cells
+```julia
+julia> modifier = @hopping!((t, i, j; dir = 1) --> (cell(i) - cell(j))[dir])
+HoppingModifier{ParametricFunction{3}}:
+  Region            : any
+  Sublattice pairs  : any
+  Cell distances    : any
+  Hopping range     : Inf
+  Reverse hops      : false
+  Argument type     : non-spatial
+  Parameters        : [:dir]
+```
+
+Note that we use the special syntax `-->` instead of `->`. This indicates that the positional arguments of the function, here called `i` and `j`, are no longer site positions as up to now. These `i, j` are non-spatial arguments, as noted by the `Argument type` property shown above. Instead of a position, a non-spatial argument `i` represent an individual site, whose index is `ind(i)`, its position is `pos(i)` and the cell it occupies on the lattice is `cell(i)`.
+
+Technically `i` is of type `CellSitePos`, which is an internal type not meant for the end user to instantiate. One special property of this type, however, is that it can efficiently index `OrbitalSliceArray`s. We can use this to build a Hamiltonian that depends on an observable, such as a `densitymatrix`. A simple example of a four-site chain with onsite energies shifted by a potential proportional to the local density on each site:
+```julia
+julia> h = LP.linear() |> onsite(2) - hopping(1) |> supercell(4) |> supercell;
+
+julia> h(SA[])
+4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+  2.0+0.0im  -1.0+0.0im       ⋅           ⋅
+ -1.0+0.0im   2.0+0.0im  -1.0+0.0im       ⋅
+      ⋅      -1.0+0.0im   2.0+0.0im  -1.0+0.0im
+      ⋅           ⋅      -1.0+0.0im   2.0+0.0im
+
+julia> g = greenfunction(h, GS.Spectrum());
+
+julia> ρ = densitymatrix(g[])(0.5, 0) ## density matrix at chemical potential `µ=0.5` and temperature `kBT = 0`  on all sites
+4×4 OrbitalSliceMatrix{Matrix{ComplexF64}}:
+ 0.138197+0.0im  0.223607+0.0im  0.223607+0.0im  0.138197+0.0im
+ 0.223607+0.0im  0.361803+0.0im  0.361803+0.0im  0.223607+0.0im
+ 0.223607+0.0im  0.361803+0.0im  0.361803+0.0im  0.223607+0.0im
+ 0.138197+0.0im  0.223607+0.0im  0.223607+0.0im  0.138197+0.0im
+
+julia> hρ = h |> @onsite!((o, i; U = 0, ρ) --> o + U * ρ[i])
+ParametricHamiltonian{Float64,1,0}: Parametric Hamiltonian on a 0D Lattice in 1D space
+  Bloch harmonics  : 1
+  Harmonic size    : 4 × 4
+  Orbitals         : [1]
+  Element type     : scalar (ComplexF64)
+  Onsites          : 4
+  Hoppings         : 6
+  Coordination     : 1.5
+  Parameters       : [:U, :ρ]
+
+julia> hρ(SA[]; U = 1, ρ = ρ0)
+4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+ 2.1382+0.0im    -1.0+0.0im         ⋅             ⋅
+   -1.0+0.0im  2.3618+0.0im    -1.0+0.0im         ⋅
+        ⋅        -1.0+0.0im  2.3618+0.0im    -1.0+0.0im
+        ⋅             ⋅        -1.0+0.0im  2.1382+0.0im
+```
+
+Note the `ρ[i]` above. This indexes `ρ` at site `i`. For a multiorbital hamiltonian, this will be a matrix (the local density matrix on each site `i`). Here it is just a number, either ` 0.138197` (sites 1 and 4) or `0.361803` (sites 2 and 3).
+
+The above provides the tools to implement self-consistent mean field. We just need to find a fixed point `ρf(ρ) = ρ` of the function `ρf` that produces the density matrix of the system.
+
+In the following example we use the FixedPoint.jl package. It provides a simple utility `afps(f, x0)` that computes the solution of `f(x) = x` with initial condition `x0`. The package requires `x0` to be a (real) `AbstractArray`. Note that any other fixed-point search routine that work with `AbstractArray`s should also work.
+```julia
+julia> using FixedPoint
+
+julia> ρf(hρ; µ = 0.5, kBT = 0, U = 0.1) = ρ -> densitymatrix(greenfunction(hρ(; U, ρ), GS.Spectrum())[])(µ, kBT)
+ρf (generic function with 1 method)
+
+julia> (ρsol, ρerror, iters) = @time afps(ρf(hρ; U = 0.4), real(ρ0), tol = 1e-8, ep = 0.95, vel = 0.0); @show iters, ρerror; ρsol
+  0.000627 seconds (1.91 k allocations: 255.664 KiB)
+(iters, ρerror) = (8, 2.0632459629688071e-10)
+4×4 OrbitalSliceMatrix{Matrix{ComplexF64}}:
+ 0.145836+0.0im  0.227266+0.0im  0.227266+0.0im  0.145836+0.0im
+ 0.227266+0.0im  0.354164+0.0im  0.354164+0.0im  0.227266+0.0im
+ 0.227266+0.0im  0.354164+0.0im  0.354164+0.0im  0.227266+0.0im
+ 0.145836+0.0im  0.227266+0.0im  0.227266+0.0im  0.145836+0.0im
+```
+
+!!! note "Bring your own fixed-point solution!"
+    Note that fixed-point calculations can be tricky, and the search algorithm can have a huge impact in convergence (if the problem converges at all!). For this reason, Quantica.jl does not provide built-in fixed-point routines, only the functionality to write functions such as `ρf` above.
+
+!!! tip "Sparse mean fields"
+    The method explained above to build a Hamiltonian supports all the `SiteSelector` and `HopSelector` functionality of conventional models. Therefore, although the density matrix computed above is dense, its application to the Hamiltonian is sparse: it only touches the onsite matrix elements. Likewise, we could for example use `@hopping!` with a finite `range` to apply a Fock mean field within a finite range. In the future we will support built-in Hartree-Fock model presets with adjustable sparsity.
diff --git a/docs/src/tutorial/glossary.md b/docs/src/tutorial/glossary.md
index 1f71aa74..1d5d6c1e 100644
--- a/docs/src/tutorial/glossary.md
+++ b/docs/src/tutorial/glossary.md
@@ -29,6 +29,5 @@ This is a summary of the type of objects you will be studying.
 - **`GreenFunction`**: an `OpenHamiltonian` combined with a `GreenSolver`, which is an algorithm that can in general compute the retarded or advanced Green function at any energy between any subset of sites of the underlying lattice.
   - **`GreenSlice`**: a `GreenFunction` evaluated on a specific set of sites, but at an unspecified energy
   - **`GreenSolution`**: a `GreenFunction` evaluated at a specific energy, but on an unspecified set of sites
-- **`Observable`**: a physical observable that can be expressed in terms of a `GreenFunction`.
-
-  Examples of supported observables include local density of states, current density, transmission probability, conductance and Josephson current
+- **`OrbitalSliceArray`**: an `AbstractArray` that can be indexed with a `SiteSelector`, in addition to the usual scalar indexing. Particular cases are `OrbitalSliceMatrix` and `OrbitalSliceVector`. This is the most common type obtained from `GreenFunction`s and observables obtained from them.
+- **Observables**: Supported observables, obtained from Green functions using various algorithms, include **local density of states**, **density matrices**, **current densities**, **transmission probabilities**, **conductance** and **Josephson currents**
diff --git a/docs/src/tutorial/models.md b/docs/src/tutorial/models.md
index ca67048b..38e5727a 100644
--- a/docs/src/tutorial/models.md
+++ b/docs/src/tutorial/models.md
@@ -120,7 +120,8 @@ ParametricModel: model with 1 term
     Hopping range     : Neighbors(1)
     Reverse hops      : false
     Coefficient       : 1
-    Parameters        : [:B, :t]
+    Argument type     : spatial
+    Parameters        : [:Bz, :t]
 ```
 
 One can linearly combine parametric and non-parametric models freely, omit parameter default values, and use any of the functional argument forms described for `onsite` and `hopping` (although not the constant argument form):
@@ -134,6 +135,7 @@ ParametricModel: model with 2 terms
     Hopping range     : Neighbors(1)
     Reverse hops      : false
     Coefficient       : -4
+    Argument type     : spatial
     Parameters        : [:t]
   OnsiteTerm{Int64}:
     Region            : any
@@ -142,6 +144,10 @@ ParametricModel: model with 2 terms
     Coefficient       : 2
 ```
 
+!!! tip "Non-spatial parametric models with -->"
+    The `->` in the above parametric models `@onsite` and `@hopping`, but also in the modifiers below, can be changed to `-->`. This indicates that the function arguments are no longer treated as site or link positions `r` and `dr`, but as objects `i, j` representing destination and source sites. This allows to address sites directly instead of through their spatial location. See the Mean Field section for further details.
+
+
 ## Modifiers
 
 There is a third model-related functionality known as `OnsiteModifier`s and `HoppingModifier`s. Given a model that defines a set of onsite and hopping amplitudes on a subset of sites and hops, one can define a parameter-dependent modification of a subset of said amplitudes. This is a useful way to introduce a new parameter dependence on an already defined model. Modifiers are built with
@@ -160,6 +166,7 @@ HoppingModifier{ParametricFunction{3}}:
   Cell distances    : any
   Hopping range     : Inf
   Reverse hops      : false
+  Argument type     : spatial
   Parameters        : [:B]
 ```
 The difference with `model_perierls` is that `model_perierls!` will never add any new hoppings. It will only modify previously existing hoppings in a model. Modifiers are not models themselves, and cannot be summed to other models. They are instead meant to be applied sequentially after applying a model.
diff --git a/src/Quantica.jl b/src/Quantica.jl
index 82727652..37aaaa6e 100644
--- a/src/Quantica.jl
+++ b/src/Quantica.jl
@@ -21,10 +21,9 @@ using Statistics: mean
 using QuadGK
 using SpecialFunctions
 
-export sublat, bravais_matrix, lattice, sites, supercell,
-       hopping, onsite, @onsite, @hopping, @onsite!, @hopping!, plusadjoint, neighbors,
-       siteselector, hopselector, diagonal,
-       hamiltonian,
+export sublat, bravais_matrix, lattice, sites, supercell, hamiltonian,
+       hopping, onsite, @onsite, @hopping, @onsite!, @hopping!, pos, ind, cell,
+       plusadjoint, neighbors, siteselector, hopselector, diagonal,
        unflat, torus, transform, translate, combine,
        spectrum, energies, states, bands, subdiv,
        greenfunction, selfenergy, attach, contact, cellsites,
diff --git a/src/apply.jl b/src/apply.jl
index dbe1ab8c..5ebcedc5 100644
--- a/src/apply.jl
+++ b/src/apply.jl
@@ -130,10 +130,10 @@ end
 apply(m::TightbindingModel, latos) = TightbindingModel(apply.(terms(m), Ref(latos)))
 
 apply(t::ParametricOnsiteTerm, lat::Lattice) =
-    ParametricOnsiteTerm(functor(t), apply(selector(t), lat), coefficient(t))
+    ParametricOnsiteTerm(functor(t), apply(selector(t), lat), coefficient(t), is_spatial(t))
 
 apply(t::ParametricHoppingTerm, lat::Lattice) =
-    ParametricHoppingTerm(functor(t), apply(selector(t), lat), coefficient(t))
+    ParametricHoppingTerm(functor(t), apply(selector(t), lat), coefficient(t), is_spatial(t))
 
 apply(m::ParametricModel, lat) = ParametricModel(apply.(terms(m), Ref(lat)))
 
@@ -152,7 +152,8 @@ function apply(m::OnsiteModifier, h::Hamiltonian, shifts = missing)
     asel = apply(sel, lattice(h))
     ptrs = pointers(h, asel, shifts)
     B = blocktype(h)
-    return AppliedOnsiteModifier(sel, B, f, ptrs)
+    spatial = is_spatial(m)
+    return AppliedOnsiteModifier(sel, B, f, ptrs, spatial)
 end
 
 function apply(m::HoppingModifier, h::Hamiltonian, shifts = missing)
@@ -161,14 +162,16 @@ function apply(m::HoppingModifier, h::Hamiltonian, shifts = missing)
     asel = apply(sel, lattice(h))
     ptrs = pointers(h, asel, shifts)
     B = blocktype(h)
-    return AppliedHoppingModifier(sel, B, f, ptrs)
+    spatial = is_spatial(m)
+    return AppliedHoppingModifier(sel, B, f, ptrs, spatial)
 end
 
-function pointers(h::Hamiltonian{T,E}, s::AppliedSiteSelector{T,E}, shifts) where {T,E}
+function pointers(h::Hamiltonian{T,E,L,B}, s::AppliedSiteSelector{T,E,L}, shifts) where {T,E,L,B}
     isempty(cells(s)) || argerror("Cannot constrain cells in an onsite modifier, cell periodicity is assumed.")
-    ptr_r = Tuple{Int,SVector{E,T},Int}[]
+    ptrs = Tuple{Int,SVector{E,T},CellSitePos{T,E,L,B},Int}[]
     lat = lattice(h)
     har0 = first(harmonics(h))
+    dn0 = zerocell(lat)
     umat = unflat(har0)
     rows = rowvals(umat)
     norbs = norbitals(h)
@@ -179,19 +182,20 @@ function pointers(h::Hamiltonian{T,E}, s::AppliedSiteSelector{T,E}, shifts) wher
         r = apply_shift(shifts, r, col)
         if (scol, r) in s
             n = norbs[scol]
-            push!(ptr_r, (p, r, n))
+            sp = CellSitePos(dn0, col, r, B)
+            push!(ptrs, (p, r, sp, n))
         end
     end
-    return ptr_r
+    return ptrs
 end
 
-function pointers(h::Hamiltonian{T,E}, s::AppliedHopSelector{T,E}, shifts) where {T,E}
+function pointers(h::Hamiltonian{T,E,L,B}, s::AppliedHopSelector{T,E,L}, shifts) where {T,E,L,B}
     hars = harmonics(h)
-    ptr_r_dr = [Tuple{Int,SVector{E,T},SVector{E,T},Tuple{Int,Int}}[] for _ in hars]
+    ptrs = [Tuple{Int,SVector{E,T},SVector{E,T},CellSitePos{T,E,L,B},CellSitePos{T,E,L,B},Tuple{Int,Int}}[] for _ in hars]
     lat = lattice(h)
     dn0 = zerocell(lat)
     norbs = norbitals(h)
-    for (har, ptr_r_dr) in zip(hars, ptr_r_dr)
+    for (har, ptrs) in zip(hars, ptrs)
         mh = unflat(har)
         rows = rowvals(mh)
         for scol in sublats(lat), col in siterange(lat, scol), p in nzrange(mh, col)
@@ -205,11 +209,12 @@ function pointers(h::Hamiltonian{T,E}, s::AppliedHopSelector{T,E}, shifts) where
             if (scol => srow, (r, dr), dn) in s
                 ncol = norbs[scol]
                 nrow = norbs[srow]
-                push!(ptr_r_dr, (p, r, dr, (nrow, ncol)))
+                srow, scol = CellSitePos(dn, row, rrow, B), CellSitePos(dn0, col, rcol, B)
+                push!(ptrs, (p, r, dr, srow, scol, (nrow, ncol)))
             end
         end
     end
-    return ptr_r_dr
+    return ptrs
 end
 
 apply_shift(::Missing, r, _) = r
@@ -288,17 +293,15 @@ function apply_map(mapping, hf::Function, ::Type{S}) where {T,S<:SVector{<:Any,T
     return FunctionWrapper{SparseMatrixCSC{Complex{T},Int},Tuple{S}}(sfunc)
 end
 
-
 #endregion
 
-
 ############################################################################################
 # apply CellSites
 #region
 
-apply(c::CellSites{L,Vector{Int}}, ::Lattice{<:Any,<:Any,L}) where {L} = c
+apply(c::AnyCellSite, ::Lattice{<:Any,<:Any,L}) where {L} = c
 
-apply(c::CellSites{L,Int}, ::Lattice{<:Any,<:Any,L}) where {L} = c
+apply(c::CellSites{L,Vector{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)))
diff --git a/src/docstrings.jl b/src/docstrings.jl
index dc9a332f..64399528 100644
--- a/src/docstrings.jl
+++ b/src/docstrings.jl
@@ -779,6 +779,16 @@ params...)`, to obtain a regular `Hamiltonian` without reconstructing it from sc
 
 Special form of a parametric onsite potential meant to model a self-energy (see `attach`).
 
+    @onsite((i; params...) --> ...; sites...)
+    @onsite((ω, i; params...) --> ...; sites...)
+
+The `-->` syntax allows to treat the argument `i` as a site index, instead of a position. In
+fact, the type of `i` is `CellSitePos`, so they can be used to index `OrbitalSliceArray`s
+(see doctrings for details). The functions `pos(i)`, `cell(i)` and `ind(i)` yield the
+position, cell and site index of the site. This syntax is useful to implement models that
+depend on observables (in the form of `OrbitalSliceArray`s), like in self-consistent mean
+field calculations.
+
 ## Model algebra
 
 Parametric models can be combined using `+`, `-` and `*`, or conjugated with `'`, e.g.
@@ -794,12 +804,14 @@ ParametricModel: model with 2 terms
     Sublattices       : A
     Cells             : any
     Coefficient       : 1
+    Argument type     : spatial
     Parameters        : [:dμ]
   ParametricOnsiteTerm{ParametricFunction{0}}
     Region            : any
     Sublattices       : B
     Cells             : any
     Coefficient       : 1
+    Argument type     : spatial
     Parameters        : [:dμ]
 
 julia> LP.honeycomb() |> supercell(2) |> hamiltonian(model, orbitals = 2)
@@ -815,7 +827,7 @@ ParametricHamiltonian{Float64,2,2}: Parametric Hamiltonian on a 2D Lattice in 2D
 ```
 
 # See also
-    `onsite`, `hopping`, `@hopping`, `@onsite!`, `@hopping!`, `attach`, `hamiltonian`
+    `onsite`, `hopping`, `@hopping`, `@onsite!`, `@hopping!`, `attach`, `hamiltonian`, `OrbitalSliceArray`
 """
 macro onsite end
 
@@ -845,6 +857,17 @@ params...)`, to obtain a regular `Hamiltonian` without reconstructing it from sc
 
 Special form of a parametric hopping amplitude meant to model a self-energy (see `attach`).
 
+    @hopping((i, j; params...) --> ...; sites...)
+    @hopping((ω, i, j; params...) --> ...; sites...)
+
+The `-->` syntax allows to treat the arguments `i, j` as a site indices, instead of a
+positions. Here `i` is the destination (row) and `j` the source (column) site. In fact, the
+type of `i` and `j` is `CellSitePos`, so they can be used to index `OrbitalSliceArray`s (see
+doctrings for details). The functions `pos(i)`, `cell(i)` and `ind(i)` yield the position,
+cell and site index of the site. This syntax is useful to implement models that depend on
+observables (in the form of `OrbitalSliceArray`s), like in self-consistent mean field
+calculations.
+
 ## Model algebra
 
 Parametric models can be combined using `+`, `-` and `*`, or conjugated with `'`, e.g.
@@ -862,6 +885,7 @@ ParametricModel: model with 1 term
     Hopping range     : Neighbors(1)
     Reverse hops      : false
     Coefficient       : 1
+    Argument type     : spatial
     Parameters        : [:t, :A]
 
 julia> LP.honeycomb() |> supercell(2) |> hamiltonian(model)
@@ -877,7 +901,7 @@ ParametricHamiltonian{Float64,2,2}: Parametric Hamiltonian on a 2D Lattice in 2D
 ```
 
 # See also
-    `onsite`, `hopping`, `@onsite`, `@onsite!`, `@hopping!`, `attach`, `hamiltonian`
+    `onsite`, `hopping`, `@onsite`, `@onsite!`, `@hopping!`, `attach`, `hamiltonian`, `OrbitalSliceArray`
 """
 macro hopping end
 
@@ -900,6 +924,15 @@ will be zero, and will not be modified by any `@onsite!` modifier. Conversely, i
 model has been applied, `@onsite!` may modify the onsite potential even if it is zero. The
 same applies to `@hopping!`.
 
+    @onsite((o, i; params...) --> ...; sites...)
+
+The `-->` syntax allows to treat the argument `i` as a site index, instead of a position. In
+fact, the type of `i` is `CellSitePos`, so they can be used to index `OrbitalSliceArray`s
+(see doctrings for details). The functions `pos(i)`, `cell(i)` and `ind(i)` yield the
+position, cell and site index of the site. This syntax is useful to implement models that
+depend on observables (in the form of `OrbitalSliceArray`s), like in self-consistent mean
+field calculations.
+
 # Examples
 ```jldoctest
 julia> model = onsite(0); disorder = @onsite!((o; W = 0) -> o + W * rand())
@@ -907,6 +940,7 @@ OnsiteModifier{ParametricFunction{1}}:
   Region            : any
   Sublattices       : any
   Cells             : any
+  Argument type     : spatial
   Parameters        : [:W]
 
 julia> LP.honeycomb() |> hamiltonian(model) |> supercell(10) |> hamiltonian(disorder)
@@ -922,7 +956,7 @@ ParametricHamiltonian{Float64,2,2}: Parametric Hamiltonian on a 2D Lattice in 2D
 ```
 
 # See also
-    `onsite`, `hopping`, `@onsite`, `@hopping`, `@hopping!`, `hamiltonian`
+    `onsite`, `hopping`, `@onsite`, `@hopping`, `@hopping!`, `hamiltonian`, `OrbitalSliceArray`
 """
 macro onsite! end
 
@@ -946,6 +980,16 @@ will be zero, and will not be modified by any `@onsite!` modifier. Conversely, i
 model has been applied, `@onsite!` may modify the onsite potential even if it is zero. The
 same applies to `@hopping!`.
 
+    @hopping!((t, i, j; params...) --> ...; sites...)
+
+The `-->` syntax allows to treat the arguments `i, j` as a site indices, instead of a
+positions. Here `i` is the destination (row) and `j` the source (column) site. In fact, the
+type of `i` and `j` is `CellSitePos`, so they can be used to index `OrbitalSliceArray`s (see
+doctrings for details). The functions `pos(i)`, `cell(i)` and `ind(i)` yield the position,
+cell and site index of the site. This syntax is useful to implement models that depend on
+observables (in the form of `OrbitalSliceArray`s), like in self-consistent mean field
+calculations.
+
 # Examples
 ```jldoctest
 julia> model = hopping(1); peierls = @hopping!((t, r, dr; A = r -> SA[0,0]) -> t * cis(-dr' * A(r)))
@@ -955,6 +999,7 @@ HoppingModifier{ParametricFunction{3}}:
   Cell distances    : any
   Hopping range     : Inf
   Reverse hops      : false
+  Argument type     : spatial
   Parameters        : [:A]
 
 julia> LP.honeycomb() |> hamiltonian(model) |> supercell(10) |> hamiltonian(peierls)
@@ -970,7 +1015,7 @@ ParametricHamiltonian{Float64,2,2}: Parametric Hamiltonian on a 2D Lattice in 2D
 ```
 
 # See also
-    `onsite`, `hopping`, `@onsite`, `@hopping`, `@onsite!`, `hamiltonian`
+    `onsite`, `hopping`, `@onsite`, `@hopping`, `@onsite!`, `hamiltonian`, `OrbitalSliceArray`
 """
 macro hopping! end
 
diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl
index b3f8e9a1..314bc9ea 100644
--- a/src/hamiltonian.jl
+++ b/src/hamiltonian.jl
@@ -92,6 +92,8 @@ lattice(b::AbstractHamiltonianBuilder) = b.lat
 
 blockstructure(b::AbstractHamiltonianBuilder) = b.blockstruct
 
+blocktype(::AbstractHamiltonianBuilder{<:Any,<:Any,<:Any,B}) where {B} = B
+
 harmonics(b::AbstractHamiltonianBuilder) = b.harmonics
 
 function Base.getindex(b::AbstractHamiltonianBuilder{<:Any,<:Any,L}, dn::SVector{L,Int}) where {L}
@@ -166,11 +168,13 @@ function applyterm!(builder, block, term::AppliedOnsiteTerm)
     ijv = builder[dn0]
     bs = blockstructure(builder)
     bsizes = blocksizes(bs)
+    B = blocktype(builder)
     foreach_site(sel, dn0) do s, i, r
         isinblock(i, block) || return nothing
         n = bsizes[s]
-        v = term(r, n)
-        push!(ijv, (i, i, v))
+        # conventional terms are never non-spatial, only modifiers can be
+        vr = term(r, n)
+        push!(ijv, (i, i, vr))
     end
     return nothing
 end
@@ -180,14 +184,16 @@ function applyterm!(builder, block, term::AppliedHoppingTerm)
     sel = selector(term)
     bs = blockstructure(builder)
     bsizes = blocksizes(bs)
+    B = blocktype(builder)
     foreach_cell(sel) do dn
         ijv = builder[dn]
         found = foreach_hop(sel, trees, dn) do (si, sj), (i, j), (r, dr)
             isinblock(i, j, block) || return nothing
             ni = bsizes[si]
             nj = bsizes[sj]
-            v = term(r, dr, (ni, nj))
-            push!(ijv, (i, j, v))
+            # conventional terms are never non-spatial, only modifiers can be
+            vr = term(r, dr, (ni, nj))
+            push!(ijv, (i, j, vr))
         end
         return found
     end
@@ -381,19 +387,34 @@ applymodifiers!(h, m::Modifier; kw...) = applymodifiers!(h, apply(m, h); kw...)
 
 function applymodifiers!(h, m::AppliedOnsiteModifier; kw...)
     nz = nonzeros(unflat(first(harmonics(h))))
-    @simd for p in pointers(m)
-        (ptr, r, norbs) = p
-        @inbounds nz[ptr] = m(nz[ptr], r, norbs; kw...)   # @inbounds too risky?
+    if is_spatial(m)    # Branch outside loop
+        @simd for p in pointers(m)
+            (ptr, r, s, norbs) = p
+            @inbounds nz[ptr] = m(nz[ptr], r, norbs; kw...)   # @inbounds too risky?
+        end
+    else
+        @simd for p in pointers(m)
+            (ptr, r, s, norbs) = p
+            @inbounds nz[ptr] = m(nz[ptr], s, norbs; kw...)
+        end
     end
     return h
 end
 
 function applymodifiers!(h, m::AppliedOnsiteModifier{B}; kw...) where {B<:SMatrixView}
     nz = nonzeros(unflat(first(harmonics(h))))
-    @simd for p in pointers(m)
-        (ptr, r, norbs) = p
-        val = view(nz[ptr], 1:norbs, 1:norbs)  # this might be suboptimal - do we need view?
-        @inbounds nz[ptr] = m(val, r, norbs; kw...)      # this allocates, currently unavoidable
+    if is_spatial(m)    # Branch outside loop
+        @simd for p in pointers(m)
+            (ptr, r, s, norbs) = p
+            val = view(nz[ptr], 1:norbs, 1:norbs)  # this might be suboptimal - do we need view?
+            @inbounds nz[ptr] = m(val, r, norbs; kw...)   # @inbounds too risky?
+        end
+    else
+        @simd for p in pointers(m)
+            (ptr, r, s, norbs) = p
+            val = view(nz[ptr], 1:norbs, 1:norbs)
+            @inbounds nz[ptr] = m(val, s, norbs; kw...)
+        end
     end
     return h
 end
@@ -401,9 +422,16 @@ end
 function applymodifiers!(h, m::AppliedHoppingModifier; kw...)
     for (har, ptrs) in zip(harmonics(h), pointers(m))
         nz = nonzeros(unflat(har))
-        @simd for p in ptrs
-            (ptr, r, dr, orborb) = p
-            @inbounds nz[ptr] = m(nz[ptr], r, dr, orborb; kw...)
+        if is_spatial(m)    # Branch outside loop
+            @simd for p in ptrs
+                (ptr, r, dr, si, sj, orborb) = p
+                @inbounds nz[ptr] = m(nz[ptr], r, dr, orborb; kw...)
+            end
+        else
+            @simd for p in ptrs
+                (ptr, r, dr, si, sj, orborb) = p
+                @inbounds nz[ptr] = m(nz[ptr], si, sj, orborb; kw...)
+            end
         end
     end
     return h
@@ -412,10 +440,18 @@ end
 function applymodifiers!(h, m::AppliedHoppingModifier{B}; kw...) where {B<:SMatrixView}
     for (har, ptrs) in zip(harmonics(h), pointers(m))
         nz = nonzeros(unflat(har))
-        @simd for p in ptrs
-            (ptr, r, dr, (norbs, norbs´)) = p
-            val = view(nz[ptr], 1:norbs, 1:norbs´)    # this might be suboptimal - do we need view?
-            @inbounds nz[ptr] = m(val, r, dr, (norbs, norbs´); kw...) # this allocates, unavoidable
+        if is_spatial(m)    # Branch outside loop
+            @simd for p in ptrs
+                (ptr, r, dr, si, sj, (oi, oj)) = p
+                val = view(nz[ptr], 1:oi, 1:oj)  # this might be suboptimal - do we need view?
+                @inbounds nz[ptr] = m(val, r, dr, (oi, oj); kw...)
+            end
+        else
+            @simd for p in ptrs
+                (ptr, r, dr, si, sj, (oi, oj)) = p
+                val = view(nz[ptr], 1:oi, 1:oj)
+                @inbounds nz[ptr] = m(val, si, sj, (oi, oj); kw...)
+            end
         end
     end
     return h
@@ -625,17 +661,17 @@ harmonics_map(h, h´, S) = [last(harmonic_index(h´, S*dcell(har))) for har in h
 
 function stitch_modifier(m::AppliedOnsiteModifier, ptrmap, _)
     ptrs´ = first(ptrmap)
-    p´ = [(ptrs´[ptr], r, orbs) for (ptr, r, orbs) in pointers(m)]
+    p´ = [(ptrs´[ptr], r, s, orbs) for (ptr, r, s, orbs) in pointers(m)]
     return AppliedOnsiteModifier(m, p´)
 end
 
 function stitch_modifier(m::AppliedHoppingModifier, ptrmap, harmap)
     ps = pointers(m)
     ps´ = [similar(first(ps), 0) for _ in 1:maximum(harmap)]
-    for (i, p) in enumerate(ps), (ptr, r, dr, orborbs) in p
+    for (i, p) in enumerate(ps), (ptr, r, dr, si, sj, orborbs) in p
         i´ = harmap[i]
         ptrs´ = ptrmap[i]
-        push!(ps´[i´], (ptrs´[ptr], r, dr, orborbs))
+        push!(ps´[i´], (ptrs´[ptr], r, dr, si, sj, orborbs))
     end
     sort!.(ps´)
     check_ptr_duplicates(first(ps´))
diff --git a/src/models.jl b/src/models.jl
index 14c498a3..70e5328b 100644
--- a/src/models.jl
+++ b/src/models.jl
@@ -22,9 +22,9 @@ _onsite(o::OnsiteTerm; kw...) =
 _hopping(t::HoppingTerm; kw...) =
     TightbindingModel(HoppingTerm(functor(t), hopselector(selector(t); kw...), coefficient(t)))
 _onsite(o::ParametricOnsiteTerm; kw...) =
-    ParametricModel(ParametricOnsiteTerm(functor(o), siteselector(selector(o); kw...), coefficient(o)))
+    ParametricModel(ParametricOnsiteTerm(functor(o), siteselector(selector(o); kw...), coefficient(o), is_spatial(o)))
 _hopping(t::ParametricHoppingTerm; kw...) =
-    ParametricModel(ParametricHoppingTerm(functor(t), hopselector(selector(t); kw...), coefficient(t)))
+    ParametricModel(ParametricHoppingTerm(functor(t), hopselector(selector(t); kw...), coefficient(t), is_spatial(t)))
 
 #endregion
 
@@ -42,59 +42,63 @@ _hopping(t::ParametricHoppingTerm; kw...) =
 
 # version with site selector kwargs
 macro onsite(kw, f)
-    f, N, params = get_f_N_params(f, "Only @onsite(args -> body; kw...) syntax supported. Mind the `;`.")
+    f, N, params, spatial = parse_term(f, "Only @onsite(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
     return esc(:(Quantica.ParametricModel(Quantica.ParametricOnsiteTerm(
-        Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector($kw), 1))))
+        Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector($kw), 1, $(spatial)))))
 end
 
 # version without site selector kwargs
 macro onsite(f)
-    f, N, params = get_f_N_params(f, "Only @onsite(args -> body; kw...) syntax supported.  Mind the `;`.")
+    f, N, params, spatial = parse_term(f, "Only @onsite(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
     return esc(:(Quantica.ParametricModel(Quantica.ParametricOnsiteTerm(
-            Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector(), 1))))
+            Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector(), 1, $(spatial)))))
 end
 
 # version with hop selector kwargs
 macro hopping(kw, f)
-    f, N, params = get_f_N_params(f, "Only @hopping(args -> body; kw...) syntax supported. Mind the `;`.")
+    f, N, params, spatial = parse_term(f, "Only @hopping(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
     return esc(:(Quantica.ParametricModel(Quantica.ParametricHoppingTerm(
-        Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector($kw), 1))))
+        Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector($kw), 1, $(spatial)))))
 end
 
 # version without hop selector kwargs
 macro hopping(f)
-    f, N, params = get_f_N_params(f, "Only @hopping(args -> body; kw...) syntax supported. Mind the `;`.")
+    f, N, params, spatial = parse_term(f, "Only @hopping(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
     return esc(:(Quantica.ParametricModel(Quantica.ParametricHoppingTerm(
-        Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector(), 1))))
+        Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector(), 1, $(spatial)))))
 end
 
 ## Model modifiers ##
 
 macro onsite!(kw, f)
-    f, N, params = get_f_N_params(f, "Only @onsite!(args -> body; kw...) syntax supported. Mind the `;`.")
-    return esc(:(Quantica.OnsiteModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector($kw))))
+    f, N, params, spatial = parse_term(f, "Only @onsite!(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
+    return esc(:(Quantica.OnsiteModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector($kw), $(spatial))))
 end
 
 macro onsite!(f)
-    f, N, params = get_f_N_params(f, "Only @onsite!(args -> body; kw...) syntax supported.  Mind the `;`.")
-    return esc(:(Quantica.OnsiteModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector())))
+    f, N, params, spatial = parse_term(f, "Only @onsite!(args -> body; kw...) syntax supported (or with -->).  Mind the `;`.")
+    return esc(:(Quantica.OnsiteModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.siteselector(), $(spatial))))
 end
 
 # Since the default hopping range is neighbors(1), we need change the default to Inf for @hopping!
 macro hopping!(kw, f)
-    f, N, params = get_f_N_params(f, "Only @hopping!(args -> body; kw...) syntax supported. Mind the `;`.")
-    return esc(:(Quantica.HoppingModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector_infrange($kw))))
+    f, N, params, spatial = parse_term(f, "Only @hopping!(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
+    return esc(:(Quantica.HoppingModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector_infrange($kw), $(spatial))))
 end
 
 macro hopping!(f)
-    f, N, params = get_f_N_params(f, "Only @hopping!(args -> body; kw...) syntax supported. Mind the `;`.")
-    return esc(:(Quantica.HoppingModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector_infrange())))
+    f, N, params, spatial = parse_term(f, "Only @hopping!(args -> body; kw...) syntax supported (or with -->). Mind the `;`.")
+    return esc(:(Quantica.HoppingModifier(Quantica.ParametricFunction{$N}($f, $(params)), Quantica.hopselector_infrange(), $(spatial))))
 end
 
 # Extracts normalized f, number of arguments and kwarg names from an anonymous function f
-function get_f_N_params(f, msg)
-    (f isa Expr && f.head == :->) || throw(ArgumentError(msg))
+function parse_term(f, msg)
+    (f isa Expr && (f.head == :-> || f.head == :-->)) || throw(ArgumentError(msg))
+    # change --> to -> and record change in spatial
+    spatial = f.head == :->
+    !spatial && (f.head = :->)
     d = ExprTools.splitdef(f)
+    # process keyword arguments, add splat
     kwargs = convert(Vector{Any}, get!(d, :kwargs, []))
     d[:kwargs] = kwargs  # in case it wasn't Vector{Any} originally
     if isempty(kwargs)
@@ -110,7 +114,7 @@ function get_f_N_params(f, msg)
     end
     N = haskey(d, :args) ? length(d[:args]) : 0
     f´ = ExprTools.combinedef(d)
-    return f´, N, params
+    return f´, N, params, spatial
 end
 
 get_kwname(x::Symbol) = x
@@ -122,6 +126,7 @@ hopselector_infrange(; kw...) = hopselector(; range = Inf, kw...)
 
 ############################################################################################
 # @onsite, @hopping conversions
+#   each ParametricTerm gets converted, using `basemodel` and `modifier`, into two pieces
 #region
 
 zero_model(term::ParametricOnsiteTerm) =
@@ -132,13 +137,13 @@ zero_model(term::ParametricHoppingTerm) =
 function modifier(term::ParametricOnsiteTerm{N}) where {N}
     f = (o, args...; kw...) -> o + term(args...; kw...)
     pf = ParametricFunction{N+1}(f, parameters(term))
-    return OnsiteModifier(pf, selector(term))
+    return OnsiteModifier(pf, selector(term), is_spatial(term))
 end
 
 function modifier(term::ParametricHoppingTerm{N}) where {N}
     f = (t, args...; kw...) -> t + term(args...; kw...)
     pf = ParametricFunction{N+1}(f, parameters(term))
-    return HoppingModifier(pf, selector(term))
+    return HoppingModifier(pf, selector(term), is_spatial(term))
 end
 
 basemodel(m::ParametricModel) = nonparametric(m) + TightbindingModel(zero_model.(terms(m)))
@@ -153,14 +158,14 @@ function model_ω_to_param(term::ParametricOnsiteTerm{N}, default = 0) where {N}
     # parameters(term) only needed for reporting, we omit adding :ω_internal
     f = (args...; ω_internal = default, kw...) -> term(ω_internal, args...; kw...)
     pf = ParametricFunction{N-1}(f, parameters(term))
-    return ParametricOnsiteTerm(pf, selector(term), coefficient(term))
+    return ParametricOnsiteTerm(pf, selector(term), coefficient(term), is_spatial(term))
 end
 
 function model_ω_to_param(term::ParametricHoppingTerm{N}, default = 0) where {N}
     # parameters(term) only needed for reporting, we omit adding :ω_internal
     f = (args...; ω_internal = default, kw...) -> term(ω_internal, args...; kw...)
     pf = ParametricFunction{N-1}(f, parameters(term))
-    return ParametricHoppingTerm(pf, selector(term), coefficient(term))
+    return ParametricHoppingTerm(pf, selector(term), coefficient(term), is_spatial(term))
 end
 
 #endregion
@@ -189,4 +194,4 @@ function blockindices(hams::NTuple{N,<:Any}) where {N}
     return inds
 end
 
-#endregion
\ No newline at end of file
+#endregion
diff --git a/src/observables.jl b/src/observables.jl
index 8170ce2d..eb7d80ad 100644
--- a/src/observables.jl
+++ b/src/observables.jl
@@ -44,6 +44,10 @@ check_same_contact_slice(gs) = (rows(gs) isa Integer && cols(gs) === rows(gs)) |
 check_different_contact_slice(gs) = (rows(gs) isa Integer && cols(gs) != rows(gs)) ||
     argerror("Please use a Green slice of the form `g[i::Integer, j::Integer] with `i ≠ j`")
 
+check_nodiag_axes(gs::GreenSlice) = check_nodiag_axes(rows(gs)), check_nodiag_axes(rows(gs))
+check_nodiag_axes(::DiagIndices) = argerror("Diagonal indexing not yet supported for this function")
+check_nodiag_axes(_) = nothing
+
 #endregion
 
 ############################################################################################
@@ -429,7 +433,8 @@ end
 (s::DensityMatrixIntegratorSolver)(mu, kBT; params...) = s.ifunc(mu, kBT; params...)
 
 # redirects to specialized method
-densitymatrix(gs::GreenSlice; kw...) = densitymatrix(solver(parent(gs)), gs::GreenSlice; kw...)
+densitymatrix(gs::GreenSlice; kw...) =
+    densitymatrix(solver(parent(gs)), gs::GreenSlice; kw...)
 
 # generic fallback
 densitymatrix(s::AppliedGreenSolver, gs::GreenSlice; kw...) =
@@ -439,6 +444,7 @@ densitymatrix(s::AppliedGreenSolver, gs::GreenSlice; kw...) =
 densitymatrix(gs::GreenSlice, ωmax::Number; opts...) = densitymatrix(gs, (-ωmax, ωmax); opts...)
 
 function densitymatrix(gs::GreenSlice{T}, (ωmin, ωmax)::Tuple; omegamap = Returns((;)), imshift = missing, atol = 1e-7, opts...) where {T}
+    check_nodiag_axes(gs)
     result = similar_Matrix(gs)
     opts´ = (; omegamap, imshift, slope = 1, post = gf_to_rho!, atol, opts...)
     integratorfunc(mu, kBT; params...) = iszero(kBT) ?
@@ -516,6 +522,7 @@ end
 (s::JosephsonIntegratorSolver)(kBT; params...) = s.ifunc(kBT; params...)
 
 function josephson(gs::GreenSlice{T}, ωmax; omegamap = Returns((;)), phases = missing, imshift = missing, atol = 1e-7, opts...) where {T}
+    check_nodiag_axes(gs)
     check_same_contact_slice(gs)
     contact = rows(gs)
     g = parent(gs)
diff --git a/src/sanitizers.jl b/src/sanitizers.jl
index bab08b41..3160e9f2 100644
--- a/src/sanitizers.jl
+++ b/src/sanitizers.jl
@@ -101,6 +101,19 @@ end
 
 #endregion
 
+############################################################################################
+# block sanitizers
+#   if a is the result of indexing an OrbitalSliceArray with an CellSitePos, ensure it can
+#   be the result of a model term function. Required for e.g. mean-field models.
+#region
+
+sanitize_block(::Type{C}, a) where {C<:Number} = complex(C)(only(a))
+sanitize_block(::Type{C}, a) where {C<:SMatrix} = C(a)
+# here we assume a is already of the correct size and let if fail later otherwise
+sanitize_block(::Type{C}, a) where {C<:SMatrixView} = a
+
+#endregion
+
 ############################################################################################
 # Supercell sanitizers
 #region
diff --git a/src/show.jl b/src/show.jl
index 910ffd1b..647cb86d 100644
--- a/src/show.jl
+++ b/src/show.jl
@@ -140,6 +140,7 @@ function Base.show(io::IO, t::Union{AbstractModelTerm,Modifier})
         print(io, "\n", "$(i)  Coefficient       : $(t.coefficient)")
     end
     if t isa AbstractParametricTerm || t isa Modifier
+        print(io, "\n", "$(i)  Argument type     : $(display_argument_type(t))")
         print(io, "\n", "$(i)  Parameters        : $(parameters(t))")
     end
 end
@@ -153,6 +154,8 @@ Base.summary(::ParametricHoppingTerm{N}) where {N} = "ParametricHoppingTerm{Para
 Base.summary(::OnsiteModifier{N}) where {N} = "OnsiteModifier{ParametricFunction{$N}}:"
 Base.summary(::HoppingModifier{N}) where {N} = "HoppingModifier{ParametricFunction{$N}}:"
 
+display_argument_type(t) = is_spatial(t) ? "spatial" : "non-spatial"
+
 #endregion
 
 ############################################################################################
@@ -475,3 +478,21 @@ Base.summary(::CurrentDensitySlice{T}) where {T} =
     "CurrentDensitySlice{$T} : current density at a fixed location and arbitrary energy"
 
 #endregion
+
+############################################################################################
+# OrbitalSliceMatrix
+#region
+
+# For simplified printing of the array typename
+
+function Base.showarg(io::IO, ::OrbitalSliceMatrix{<:Any,M}, toplevel) where {M}
+    toplevel || print(io, "::")
+    print(io,  "OrbitalSliceMatrix{$M}")
+end
+
+function Base.showarg(io::IO, ::OrbitalSliceVector{<:Any,M}, toplevel) where {M}
+    toplevel || print(io, "::")
+    print(io,  "OrbitalSliceVector{$M}")
+end
+
+#endregion
diff --git a/src/solvers/green/kpm.jl b/src/solvers/green/kpm.jl
index e027edc7..8571b5ea 100644
--- a/src/solvers/green/kpm.jl
+++ b/src/solvers/green/kpm.jl
@@ -233,6 +233,7 @@ end
 ## Constructor
 
 function densitymatrix(s::AppliedKPMGreenSolver, gs::GreenSlice{T}) where {T}
+    check_nodiag_axes(gs)
     has_selfenergy(gs) && argerror("The KPM densitymatrix solver currently support only `nothing` contacts")
     momenta = slice_momenta(s.momenta, gs)
     solver = DensityMatrixKPMSolver(momenta, s.bandCH)
diff --git a/src/solvers/green/spectrum.jl b/src/solvers/green/spectrum.jl
index dbc23277..cf454b94 100644
--- a/src/solvers/green/spectrum.jl
+++ b/src/solvers/green/spectrum.jl
@@ -82,6 +82,7 @@ end
 function densitymatrix(s::AppliedSpectrumGreenSolver, gs::GreenSlice)
     # SpectrumGreenSlicer is 0D, so there is a single cellorbs in dict.
     # If rows/cols are contacts, we need their orbrows/orbcols (unlike for gs(ω; params...))
+    check_nodiag_axes(gs)
     i, j = only(cellsdict(orbrows(gs))), only(cellsdict(orbcols(gs)))
     oi, oj = orbindices(i), orbindices(j)
     es, psis = spectrum(s)
diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl
index 2dc9a58b..5919f456 100644
--- a/src/specialmatrices.jl
+++ b/src/specialmatrices.jl
@@ -480,3 +480,106 @@ checkcontactindices(allcontactinds, hdim) = maximum(allcontactinds) <= hdim ||
     internalerror("InverseGreenBlockSparse: unexpected contact indices beyond Hamiltonian dimension")
 
 #endregion
+
+############################################################################################
+## OrbitalSliceArray
+#region
+
+# AbstractArray interface
+Base.size(a::OrbitalSliceArray) = size(parent(a))
+Base.iterate(a::OrbitalSliceArray, i...) = iterate(parent(a), i...)
+Base.length(a::OrbitalSliceArray) = length(parent(a))
+Base.IndexStyle(::Type{T}) where {M,T<:OrbitalSliceArray{<:Any,<:Any,M}} = IndexStyle(M)
+Base.similar(a::OrbitalSliceArray) = OrbitalSliceArray(similar(parent(a)), orbaxes(a))
+Base.similar(a::OrbitalSliceArray, t::Type) = OrbitalSliceArray(similar(parent(a), t), orbaxes(a))
+# doesn't make sense to keep orbaxes in similar with different dimensions.
+Base.similar(a::OrbitalSliceArray, dims::Tuple) = similar(parent(a), dims)
+Base.copy(a::OrbitalSliceArray) = OrbitalSliceArray(copy(parent(a)), orbaxes(a))
+Base.@propagate_inbounds Base.getindex(a::OrbitalSliceArray, i::Int) =
+    getindex(parent(a), i)
+Base.@propagate_inbounds Base.getindex(a::OrbitalSliceArray, I::Vararg{Int, N}) where {N} =
+    getindex(parent(a), I...)
+Base.@propagate_inbounds Base.setindex!(a::OrbitalSliceArray, v, i::Int) = setindex!(parent(a), v, i)
+Base.@propagate_inbounds Base.setindex!(a::OrbitalSliceArray, v, I::Vararg{Int, N}) where {N} = setindex!(parent(a), v, I...)
+
+# Additional indexing over sites
+Base.getindex(a::OrbitalSliceMatrix; sites...) = getindex(a, siteselector(; sites...))
+Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::NamedTuple = i) =
+    getindex(a, siteselector(i), siteselector(j))
+Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::SiteSelector) =
+    getindex(a, siteselector(i), j)
+Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::NamedTuple) =
+    getindex(a, i, siteselector(j))
+
+# SiteSelector: return a new OrbitalSliceMatrix
+function Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::SiteSelector = i)
+    rowslice, colslice = orbaxes(a)
+    rowslice´, colslice´ = rowslice[i], colslice[j]
+    rows = collect(indexcollection(rowslice, rowslice´))
+    cols = i === j && rowslice === colslice ? rows : indexcollection(colslice, colslice´)
+    m = parent(a)[rows, cols]
+    return OrbitalSliceMatrix(m, (rowslice´, colslice´))
+end
+
+# CellSites: return an unwrapped Matrix or a view thereof (non-allocating)
+Base.getindex(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) =
+    copy(view(a, i, j))
+
+Base.getindex(a::OrbitalSliceMatrix, i::C, j::C = i) where {B,C<:CellSitePos{<:Any,<:Any,<:Any,B}} =
+    sanitize_block(B, view(a, i, j))
+
+function Base.view(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i)
+    rowslice, colslice = orbaxes(a)
+    i´, j´ = apply(i, lattice(rowslice)), apply(j, lattice(colslice))
+    rows = indexcollection(rowslice, i´)
+    cols = j === i && rowslice === colslice ? rows : indexcollection(colslice, j´)
+    return view(parent(a), rows, cols)
+end
+
+## broadcasting
+
+# following the manual: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array
+
+Broadcast.BroadcastStyle(::Type{<:OrbitalSliceArray}) = Broadcast.ArrayStyle{OrbitalSliceArray}()
+
+Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{OrbitalSliceArray}}, ::Type{ElType}) where {ElType} =
+    OrbitalSliceArray(similar(Array{ElType}, axes(bc)), orbaxes(find_osa(bc)))
+
+find_osa(bc::Base.Broadcast.Broadcasted) = find_osa(bc.args)
+find_osa(args::Tuple) = find_osa(find_osa(args[1]), Base.tail(args))
+find_osa(x) = x
+find_osa(::Tuple{}) = nothing
+find_osa(a::OrbitalSliceArray, rest) = a
+find_osa(::Any, rest) = find_osa(rest)
+
+# taken from https://github.com/JuliaArrays/OffsetArrays.jl/blob/756e839563c88faa4ebe4ff971286747863aaff0/src/OffsetArrays.jl#L469
+
+Base.dataids(A::OrbitalSliceArray) = Base.dataids(parent(A))
+Broadcast.broadcast_unalias(dest::OrbitalSliceArray, src::OrbitalSliceArray) =
+    parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src)
+
+## conversion
+
+maybe_OrbitalSliceArray(i) = x -> maybe_OrbitalSliceArray(x, i)
+
+maybe_OrbitalSliceArray(x::AbstractVector, i) = maybe_OrbitalSliceVector(x, i)
+maybe_OrbitalSliceArray(x::AbstractMatrix, i) = maybe_OrbitalSliceMatrix(x, i)
+
+maybe_OrbitalSliceVector(x, i::DiagIndices{Missing,<:OrbitalSliceGrouped}) =
+    maybe_OrbitalSliceVector(x, parent(i))
+maybe_OrbitalSliceVector(x, i::DiagIndices{<:Any,<:OrbitalSliceGrouped}) =
+    maybe_OrbitalSliceVector(x, scalarize(parent(i)))
+maybe_OrbitalSliceVector(x, i::OrbitalSliceGrouped) = OrbitalSliceVector(x, (i,))
+
+# fallback
+maybe_OrbitalSliceVector(x, i) = x
+
+maybe_OrbitalSliceMatrix(x, i::OrbitalSliceGrouped) =
+    maybe_OrbitalSliceMatrix(x, (i, i))
+maybe_OrbitalSliceMatrix(x, (i, j)::Tuple{OrbitalSliceGrouped,OrbitalSliceGrouped}) =
+    OrbitalSliceMatrix(x, (i, j))
+
+# fallback
+maybe_OrbitalSliceMatrix(x, i) = x
+
+#endregion
diff --git a/src/types.jl b/src/types.jl
index 9bc881d3..40ad9cdf 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -262,6 +262,66 @@ Base.NamedTuple(s::HopSelector) =
 #endregion
 #endregion
 
+############################################################################################
+# MatrixElementTypes
+#region
+
+struct SMatrixView{N,M,T,NM}
+    s::SMatrix{N,M,T,NM}
+    SMatrixView{N,M,T,NM}(mat) where {N,M,T,NM} = new(sanitize_SMatrix(SMatrix{N,M,T}, mat))
+end
+
+const MatrixElementType{T} = Union{
+    Complex{T},
+    SMatrix{N,N,Complex{T}} where {N},
+    SMatrixView{N,N,Complex{T}} where {N}}
+
+const MatrixElementUniformType{T} = Union{
+    Complex{T},
+    SMatrix{N,N,Complex{T}} where {N}}
+
+const MatrixElementNonscalarType{T,N} = Union{
+    SMatrix{N,N,Complex{T}},
+    SMatrixView{N,N,Complex{T}}}
+
+#region ## Constructors ##
+
+SMatrixView(mat::SMatrix{N,M,T,NM}) where {N,M,T,NM} = SMatrixView{N,M,T,NM}(mat)
+
+SMatrixView{N,M}(mat::AbstractMatrix{T}) where {N,M,T} = SMatrixView{N,M,T}(mat)
+
+SMatrixView{N,M,T}(mat) where {N,M,T} = SMatrixView{N,M,T,N*M}(mat)
+
+SMatrixView(::Type{SMatrix{N,M,T,NM}}) where {N,M,T,NM} = SMatrixView{N,M,T,NM}
+
+#endregion
+
+#region ## API ##
+
+unblock(s::SMatrixView) = s.s
+unblock(s) = s
+
+Base.parent(s::SMatrixView) = s.s
+
+Base.view(s::SMatrixView, i...) = view(s.s, i...)
+
+Base.getindex(s::SMatrixView, i...) = getindex(s.s, i...)
+
+Base.zero(::Type{SMatrixView{N,M,T,NM}}) where {N,M,T,NM} = SMatrixView(zero(SMatrix{N,M,T,NM}))
+Base.zero(::S) where {S<:SMatrixView} = zero(S)
+
+Base.adjoint(s::SMatrixView) = SMatrixView(s.s')
+
+Base.:+(s::SMatrixView...) = SMatrixView(+(parent.(s)...))
+Base.:-(s::SMatrixView) = SMatrixView(-parent.(s))
+Base.:*(s::SMatrixView, x::Number) = SMatrixView(parent(s) * x)
+Base.:*(x::Number, s::SMatrixView) = SMatrixView(x * parent(s))
+
+Base.:(==)(s::SMatrixView, s´::SMatrixView) = parent(s) == parent(s´)
+
+#endregion
+#endregion
+
 ############################################################################################
 # LatticeSlice - see slices.jl for methods
 #   Encodes subsets of sites (or orbitals) of a lattice in different cells. Produced e.g. by
@@ -270,6 +330,11 @@ Base.NamedTuple(s::HopSelector) =
 
 struct SiteLike end
 
+struct SiteLikePos{T,E,B<:MatrixElementType{T}}
+    r::SVector{E,T}
+    blocktype::Type{B}
+end
+
 struct OrbitalLike end
 
 struct OrbitalLikeGrouped
@@ -277,16 +342,20 @@ struct OrbitalLikeGrouped
                                             # (non-selected sites are not present)
 end
 
-struct CellIndices{L,I,G<:Union{SiteLike,OrbitalLike,OrbitalLikeGrouped}}
+struct CellIndices{L,I,G<:Union{SiteLike,SiteLikePos,OrbitalLike,OrbitalLikeGrouped}}
     cell::SVector{L,Int}
     inds::I    # can be anything: Int, Colon, Vector{Int}, etc.
     type::G    # can be SiteLike, OrbitalLike or OrbitalLikeGrouped
 end
 
 const CellSites{L,I} = CellIndices{L,I,SiteLike}
+const CellSite{L} = CellIndices{L,Int,SiteLike}
+const CellSitePos{T,E,L,B} = CellIndices{L,Int,SiteLikePos{T,E,B}} # for non-spatial models
+const AnyCellSite = Union{CellSite,CellSitePos}
+const AnyCellSites = Union{CellSites,CellSitePos}
+
 const CellOrbitals{L,I} = CellIndices{L,I,OrbitalLike}
 const CellOrbitalsGrouped{L,I} = CellIndices{L,I,OrbitalLikeGrouped}
-const CellSite{L} = CellSites{L,Int}
 const CellOrbital{L} = CellIndices{L,Int,OrbitalLike}
 const AnyCellOrbitals = Union{CellOrbital,CellOrbitals,CellOrbitalsGrouped}
 
@@ -319,6 +388,7 @@ const AnyOrbitalSlice = Union{OrbitalSlice,OrbitalSliceGrouped}
 #region ## Constructors ##
 
 CellSite(cell, ind::Int) = CellIndices(sanitize_SVector(Int, cell), ind, SiteLike())
+CellSite(c::CellSitePos) = CellSite(c.cell, c.inds)
 CellSites(cell, inds = Int[]) = CellIndices(sanitize_SVector(Int, cell), sanitize_cellindices(inds), SiteLike())
 # exported lowercase constructor for general inds
 cellsites(cell, inds) = CellSites(cell, inds)
@@ -331,6 +401,9 @@ CellOrbital(cell, ind::Int) =
 CellOrbitalsGrouped(cell, inds, groups::Dictionary) =
     CellIndices(sanitize_SVector(Int, cell), sanitize_cellindices(inds), OrbitalLikeGrouped(groups))
 
+# B <: MatrixElementType{T}
+CellSitePos(cell, ind, r, B) = CellIndices(sanitize_SVector(Int, cell), ind, SiteLikePos(r, B))
+
 # LatticeSlice from an AbstractVector of CellIndices
 LatticeSlice(lat::Lattice, cs::AbstractVector{<:CellIndices}) =
     LatticeSlice(lat, cellinds_to_dict(cs))
@@ -376,9 +449,10 @@ lattice(ls::LatticeSlice) = ls.lat
 siteindices(s::CellSites) = s.inds
 siteindices(s::CellOrbitalsGrouped) = keys(s.type.groups)
 orbindices(s::AnyCellOrbitals) = s.inds
-siteindex(s::CellSite) = s.inds
+siteindex(s::AnyCellSite) = s.inds
 orbindex(s::CellOrbital) = s.inds
 
+indexcollection(l::LatticeSlice, c::CellSitePos) = siteindexdict(l)[CellSite(c)]
 indexcollection(l::LatticeSlice, c::CellSite) = siteindexdict(l)[c]
 indexcollection(l::LatticeSlice, cs::CellSites) =
     [j for i in siteindices(cs) for j in indexcollection(l, CellSite(cell(cs), i))]
@@ -414,6 +488,11 @@ findsubcell(cell::SVector, l::LatticeSlice) = findsubcell(cell, cellsdict(l))
 
 boundingbox(l::LatticeSlice) = boundingbox(keys(cellsdict(l)))
 
+# interface for non-spatial models
+pos(s::CellSitePos) = s.type.r
+ind(s::CellSitePos) = s.inds
+cell(s::CellSitePos) = s.cell
+
 # iterators
 
 sites(l::LatticeSlice) =
@@ -504,24 +583,25 @@ struct ParametricOnsiteTerm{N,S<:Union{SiteSelector,AppliedSiteSelector},F<:Para
     f::F
     selector::S
     coefficient::T
+    spatial::Bool   # If true, f is a function of position r. Otherwise it takes a single CellSite
 end
 
 struct ParametricHoppingTerm{N,S<:Union{HopSelector,AppliedHopSelector},F<:ParametricFunction{N},T<:Number} <: AbstractModelTerm
     f::F
     selector::S
     coefficient::T
+    spatial::Bool   # If true, f is a function of positions r, dr. Otherwise it takes two CellSite's
 end
 
 const AbstractParametricTerm{N} = Union{ParametricOnsiteTerm{N},ParametricHoppingTerm{N}}
+const AppliedParametricTerm{N} = Union{ParametricOnsiteTerm{N,<:AppliedSiteSelector},
+                                       ParametricHoppingTerm{N,<:AppliedHopSelector}}
 
 struct ParametricModel{T<:NTuple{<:Any,AbstractParametricTerm},M<:TightbindingModel} <: AbstractModel
     npmodel::M  # non-parametric model to use as base
     terms::T    # Collection of `AbstractParametricTerm`s
 end
 
-const AppliedParametricTerm{N} = Union{ParametricOnsiteTerm{N,<:AppliedSiteSelector},
-                                       ParametricHoppingTerm{N,<:AppliedHopSelector}}
-
 ## BlockModels ##
 
 struct InterblockModel{M<:AbstractModel,N}
@@ -570,10 +650,20 @@ coefficient(t::OnsiteTerm) = t.coefficient
 coefficient(t::HoppingTerm) = t.coefficient
 coefficient(t::AbstractParametricTerm) = t.coefficient
 
+narguments(t::OnsiteTerm{<:Function}) = 1
+narguments(t::HoppingTerm{<:Function}) = 2
+narguments(t::OnsiteTerm) = 0
+narguments(t::HoppingTerm) = 0
+narguments(t::AbstractParametricTerm) = narguments(t.f)
+narguments(::ParametricFunction{N}) where {N} = N
+
 Base.parent(m::InterblockModel) = m.model
 
 block(m::InterblockModel) = m.block
 
+is_spatial(t::AbstractParametricTerm) = t.spatial
+is_spatial(t) = true
+
 ## call API##
 
 (term::OnsiteTerm{<:Function})(r) = term.coefficient * term.f(r)
@@ -596,12 +686,12 @@ block(m::InterblockModel) = m.block
 # orbital structure, not only to a site selection
 function (t::ParametricOnsiteTerm{N})(; kw...) where {N}
     f = ParametricFunction{N}((args...) -> t.f(args...; kw...)) # no params
-    return ParametricOnsiteTerm(f, t.selector, t.coefficient)
+    return ParametricOnsiteTerm(f, t.selector, t.coefficient, t.spatial)
 end
 
 function (t::ParametricHoppingTerm{N})(; kw...) where {N}
     f = ParametricFunction{N}((args...) -> t.f(args...; kw...)) # no params
-    return ParametricHoppingTerm(f, t.selector, t.coefficient)
+    return ParametricHoppingTerm(f, t.selector, t.coefficient, t.spatial)
 end
 
 ## Model term algebra
@@ -623,9 +713,9 @@ Base.:-(m::AbstractModel, m´::AbstractModel) = m + (-m´)
 Base.:*(x::Number, o::OnsiteTerm) = OnsiteTerm(o.f, o.selector, x * o.coefficient)
 Base.:*(x::Number, t::HoppingTerm) = HoppingTerm(t.f, t.selector, x * t.coefficient)
 Base.:*(x::Number, o::ParametricOnsiteTerm) =
-    ParametricOnsiteTerm(o.f, o.selector, x * o.coefficient)
+    ParametricOnsiteTerm(o.f, o.selector, x * o.coefficient, o.spatial)
 Base.:*(x::Number, t::ParametricHoppingTerm) =
-    ParametricHoppingTerm(t.f, t.selector, x * t.coefficient)
+    ParametricHoppingTerm(t.f, t.selector, x * t.coefficient, t.spatial)
 
 Base.adjoint(m::TightbindingModel) = TightbindingModel(adjoint.(terms(m))...)
 Base.adjoint(m::ParametricModel) = ParametricModel(adjoint.(terms(m))...)
@@ -636,12 +726,12 @@ Base.adjoint(t::HoppingTerm) = HoppingTerm(t.f', t.selector', t.coefficient')
 
 function Base.adjoint(o::ParametricOnsiteTerm{N}) where {N}
     f = ParametricFunction{N}((args...; kw...) -> o.f(args...; kw...)', o.f.params)
-    return ParametricOnsiteTerm(f, o.selector, o.coefficient')
+    return ParametricOnsiteTerm(f, o.selector, o.coefficient', o.spatial)
 end
 
 function Base.adjoint(t::ParametricHoppingTerm{N}) where {N}
     f = ParametricFunction{N}((args...; kw...) -> t.f(args...; kw...)', t.f.params)
-    return ParametricHoppingTerm(f, t.selector, t.coefficient')
+    return ParametricHoppingTerm(f, t.selector, t.coefficient', t.spatial)
 end
 
 #endregion
@@ -656,27 +746,31 @@ abstract type AbstractModifier end
 struct OnsiteModifier{N,S<:SiteSelector,F<:ParametricFunction{N}} <: AbstractModifier
     f::F
     selector::S
+    spatial::Bool
 end
 
-struct AppliedOnsiteModifier{B,N,R<:SVector,F<:ParametricFunction{N},S<:SiteSelector} <: AbstractModifier
+struct AppliedOnsiteModifier{B,N,R<:SVector,F<:ParametricFunction{N},S<:SiteSelector,P<:CellSitePos} <: AbstractModifier
     parentselector::S        # unapplied selector, needed to grow a ParametricHamiltonian
     blocktype::Type{B}  # These are needed to cast the modification to the sublat block type
     f::F
-    ptrs::Vector{Tuple{Int,R,Int}}
-    # [(ptr, r, norbs)...] for each selected site, dn = 0 harmonic
+    ptrs::Vector{Tuple{Int,R,P,Int}}
+    # [(ptr, r, si, norbs)...] for each selected site, dn = 0 harmonic
+    spatial::Bool   # If true, f is a function of position r. Otherwise it takes a single CellSite
 end
 
 struct HoppingModifier{N,S<:HopSelector,F<:ParametricFunction{N}} <: AbstractModifier
     f::F
     selector::S
+    spatial::Bool  # If true, f is a function of positions r, dr. Otherwise it takes two CellSite's
 end
 
-struct AppliedHoppingModifier{B,N,R<:SVector,F<:ParametricFunction{N},S<:HopSelector} <: AbstractModifier
+struct AppliedHoppingModifier{B,N,R<:SVector,F<:ParametricFunction{N},S<:HopSelector,P<:CellSitePos} <: AbstractModifier
     parentselector::S        # unapplied selector, needed to grow a ParametricHamiltonian
     blocktype::Type{B}  # These are needed to cast the modification to the sublat block type
     f::F
-    ptrs::Vector{Vector{Tuple{Int,R,R,Tuple{Int,Int}}}}
-    # [[(ptr, r, dr, (norbs, norbs´)), ...], ...] for each selected hop on each harmonic
+    ptrs::Vector{Vector{Tuple{Int,R,R,P,P,Tuple{Int,Int}}}}
+    # [[(ptr, r, dr, si, sj, (norbs, norbs´)), ...], ...] for each selected hop on each harmonic
+    spatial::Bool  # If true, f is a function of positions r, dr. Otherwise it takes two CellSite's
 end
 
 const Modifier = Union{OnsiteModifier,HoppingModifier}
@@ -685,10 +779,10 @@ const AppliedModifier = Union{AppliedOnsiteModifier,AppliedHoppingModifier}
 #region ## Constructors ##
 
 AppliedOnsiteModifier(m::AppliedOnsiteModifier, ptrs) =
-    AppliedOnsiteModifier(m.parentselector, m.blocktype, m.f, ptrs)
+    AppliedOnsiteModifier(m.parentselector, m.blocktype, m.f, ptrs, m.spatial)
 
 AppliedHoppingModifier(m::AppliedHoppingModifier, ptrs) =
-    AppliedHoppingModifier(m.parentselector, m.blocktype, m.f, ptrs)
+    AppliedHoppingModifier(m.parentselector, m.blocktype, m.f, ptrs, m.spatial)
 
 #endregion
 
@@ -705,6 +799,10 @@ pointers(m::AppliedModifier) = m.ptrs
 
 blocktype(m::AppliedModifier) = m.blocktype
 
+is_spatial(m::AbstractModifier) = m.spatial
+
+narguments(m::AbstractModifier) = narguments(m.f)
+
 @inline (m::AppliedOnsiteModifier{B,1})(o, r, orbs; kw...) where {B} =
     mask_block(B, m.f.f(o; kw...), (orbs, orbs))
 @inline (m::AppliedOnsiteModifier{B,2})(o, r, orbs; kw...) where {B} =
@@ -715,10 +813,10 @@ blocktype(m::AppliedModifier) = m.blocktype
 @inline (m::AppliedHoppingModifier{B,3})(t, r, dr, orborb; kw...) where {B} =
     mask_block(B, m.f.f(t, r, dr; kw...), orborb)
 
-Base.similar(m::A) where {A <: AppliedModifier} = A(m.blocktype, m.f, similar(m.ptrs, 0))
+Base.similar(m::A) where {A <: AppliedModifier} = A(m.blocktype, m.f, similar(m.ptrs, 0), m.spatial)
 
-Base.parent(m::AppliedOnsiteModifier) = OnsiteModifier(m.f, m.parentselector)
-Base.parent(m::AppliedHoppingModifier) = HoppingModifier(m.f, m.parentselector)
+Base.parent(m::AppliedOnsiteModifier) = OnsiteModifier(m.f, m.parentselector, m.spatial)
+Base.parent(m::AppliedHoppingModifier) = HoppingModifier(m.f, m.parentselector, m.spatial)
 
 #endregion
 #endregion
@@ -728,11 +826,6 @@ Base.parent(m::AppliedHoppingModifier) = HoppingModifier(m.f, m.parentselector)
 #    Block structure for Hamiltonians, sorted by sublattices
 #region
 
-struct SMatrixView{N,M,T,NM}
-    s::SMatrix{N,M,T,NM}
-    SMatrixView{N,M,T,NM}(mat) where {N,M,T,NM} = new(sanitize_SMatrix(SMatrix{N,M,T}, mat))
-end
-
 struct OrbitalBlockStructure{B}
     blocksizes::Vector{Int}       # number of orbitals per site in each sublattice
     subsizes::Vector{Int}         # number of blocks (sites) in each sublattice
@@ -744,29 +837,8 @@ struct OrbitalBlockStructure{B}
     end
 end
 
-const MatrixElementType{T} = Union{
-    Complex{T},
-    SMatrix{N,N,Complex{T}} where {N},
-    SMatrixView{N,N,Complex{T}} where {N}}
-
-const MatrixElementUniformType{T} = Union{
-    Complex{T},
-    SMatrix{N,N,Complex{T}} where {N}}
-
-const MatrixElementNonscalarType{T,N} = Union{
-    SMatrix{N,N,Complex{T}},
-    SMatrixView{N,N,Complex{T}}}
-
 #region ## Constructors ##
 
-SMatrixView(mat::SMatrix{N,M,T,NM}) where {N,M,T,NM} = SMatrixView{N,M,T,NM}(mat)
-
-SMatrixView{N,M}(mat::AbstractMatrix{T}) where {N,M,T} = SMatrixView{N,M,T}(mat)
-
-SMatrixView{N,M,T}(mat) where {N,M,T} = SMatrixView{N,M,T,N*M}(mat)
-
-SMatrixView(::Type{SMatrix{N,M,T,NM}}) where {N,M,T,NM} = SMatrixView{N,M,T,NM}
-
 @inline function OrbitalBlockStructure(T, orbitals, subsizes)
     orbitals´ = sanitize_orbitals(orbitals) # <:Union{Val,distinct_collection}
     B = blocktype(T, orbitals´)
@@ -804,7 +876,7 @@ flatsize(b::OrbitalBlockStructure) = blocksizes(b)' * subsizes(b)
 flatsize(b::OrbitalBlockStructure{B}, ls::SiteSlice) where {B<:Union{Complex,SMatrix}} =
     length(ls) * blocksize(b)
 flatsize(b::OrbitalBlockStructure, ls::SiteSlice) = sum(cs -> flatsize(b, cs), cellsites(ls); init = 0)
-flatsize(b::OrbitalBlockStructure, cs::CellSite) = blocksize(b, siteindex(cs))
+flatsize(b::OrbitalBlockStructure, cs::AnyCellSite) = blocksize(b, siteindex(cs))
 
 unflatsize(b::OrbitalBlockStructure) = sum(subsizes(b))
 
@@ -891,27 +963,6 @@ unflatindex_and_blocksize(::OrbitalBlockStructure{<:Number}, iflat::Integer) = (
 Base.copy(b::OrbitalBlockStructure{B}) where {B} =
     OrbitalBlockStructure{B}(copy(blocksizes(b)), copy(subsizes(b)))
 
-unblock(s::SMatrixView) = s.s
-unblock(s) = s
-
-Base.parent(s::SMatrixView) = s.s
-
-Base.view(s::SMatrixView, i...) = view(s.s, i...)
-
-Base.getindex(s::SMatrixView, i...) = getindex(s.s, i...)
-
-Base.zero(::Type{SMatrixView{N,M,T,NM}}) where {N,M,T,NM} = SMatrixView(zero(SMatrix{N,M,T,NM}))
-Base.zero(::S) where {S<:SMatrixView} = zero(S)
-
-Base.adjoint(s::SMatrixView) = SMatrixView(s.s')
-
-Base.:+(s::SMatrixView...) = SMatrixView(+(parent.(s)...))
-Base.:-(s::SMatrixView) = SMatrixView(-parent.(s))
-Base.:*(s::SMatrixView, x::Number) = SMatrixView(parent(s) * x)
-Base.:*(x::Number, s::SMatrixView) = SMatrixView(x * parent(s))
-
-Base.:(==)(s::SMatrixView, s´::SMatrixView) = parent(s) == parent(s´)
-
 Base.:(==)(b::OrbitalBlockStructure{B}, b´::OrbitalBlockStructure{B}) where {B} =
     b.blocksizes == b´.blocksizes && b.subsizes == b´.subsizes
 
@@ -2128,7 +2179,7 @@ Base.getindex(c::Operator, i...) = getindex(c.h, i...)
 #endregion
 
 ############################################################################################
-# OrbitalSliceArray - array type over orbital slice
+# OrbitalSliceArray: array type over orbital slice - see specialmatrices.jl
 #   orbaxes is a tuple of OrbitalSlice's or Missing's, one per dimension
 #   if missing, the dimension does not span orbitals but something else
 #region
@@ -2146,84 +2197,6 @@ OrbitalSliceMatrix(m::AbstractMatrix, axes) = OrbitalSliceArray(m, axes)
 
 orbaxes(a::OrbitalSliceArray) = a.orbaxes
 
-# AbstractArray interface
 Base.parent(a::OrbitalSliceArray) = a.parent
-Base.size(a::OrbitalSliceArray) = size(a.parent)
-Base.getindex(a::OrbitalSliceArray, i::Int) = getindex(a.parent, i)
-Base.getindex(a::OrbitalSliceArray, I::Vararg{Int, N}) where {N} = getindex(a.parent, I...)
-Base.IndexStyle(::Type{T}) where {M,T<:OrbitalSliceArray{<:Any,<:Any,M}} = IndexStyle(M)
-Base.setindex!(a::OrbitalSliceArray, v, i::Int) = setindex!(a.parent, v, i)
-Base.setindex!(a::OrbitalSliceArray, v, I::Vararg{Int, N}) where {N} = setindex!(a.parent, v, I...)
-Base.iterate(a::OrbitalSliceArray) = iterate(a.parent)
-Base.iterate(a::OrbitalSliceArray, state) = iterate(a.parent, state)
-Base.length(a::OrbitalSliceArray) = length(a.parent)
-# doesn't make sense to keep orbaxes in similar (dimensions could change)
-Base.similar(a::OrbitalSliceArray, t::Tuple) = similar(a.parent, t)
-Base.copy(a::OrbitalSliceArray) = OrbitalSliceArray(copy(a.parent), a.orbaxes)
-
-# Additional indexing over sites
-Base.getindex(a::OrbitalSliceMatrix; sites...) = getindex(a, siteselector(; sites...))
-Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::NamedTuple = i) =
-    getindex(a, siteselector(i), siteselector(j))
-Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::SiteSelector) =
-    getindex(a, siteselector(i), j)
-Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::NamedTuple) =
-    getindex(a, i, siteselector(j))
-
-# SiteSelector: return a new OrbitalSliceMatrix
-function Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::SiteSelector = i)
-    rowslice, colslice = orbaxes(a)
-    rowslice´, colslice´ = rowslice[i], colslice[j]
-    rows = collect(indexcollection(rowslice, rowslice´))
-    cols = i === j && rowslice === colslice ? rows : indexcollection(colslice, colslice´)
-    m = parent(a)[rows, cols]
-    return OrbitalSliceMatrix(m, (rowslice´, colslice´))
-end
-
-# CellSites: return an unwrapped Matrix or a view thereof (non-allocating)
-Base.getindex(a::OrbitalSliceMatrix, i::CellSites, j::CellSites = i) = copy(view(a, i, j))
-
-function Base.view(a::OrbitalSliceMatrix, i::CellSites, j::CellSites = i)
-    rowslice, colslice = orbaxes(a)
-    i´, j´ = apply(i, lattice(rowslice)), apply(j, lattice(colslice))
-    rows = indexcollection(rowslice, i´)
-    cols = j === i && rowslice === colslice ? rows : indexcollection(colslice, j´)
-    return view(a.parent, rows, cols)
-end
-
-# For simplified printing
-function Base.showarg(io::IO, ::OrbitalSliceMatrix{<:Any,M}, toplevel) where {M}
-    toplevel || print(io, "::")
-    print(io,  "OrbitalSliceMatrix{$M}")
-end
-
-function Base.showarg(io::IO, ::OrbitalSliceVector{<:Any,M}, toplevel) where {M}
-    toplevel || print(io, "::")
-    print(io,  "OrbitalSliceVector{$M}")
-end
-
-## conversion
-
-maybe_OrbitalSliceArray(i) = x -> maybe_OrbitalSliceArray(x, i)
-
-maybe_OrbitalSliceArray(x::AbstractVector, i) = maybe_OrbitalSliceVector(x, i)
-maybe_OrbitalSliceArray(x::AbstractMatrix, i) = maybe_OrbitalSliceMatrix(x, i)
-
-maybe_OrbitalSliceVector(x, i::DiagIndices{Missing,<:OrbitalSliceGrouped}) =
-    maybe_OrbitalSliceVector(x, parent(i))
-maybe_OrbitalSliceVector(x, i::DiagIndices{<:Any,<:OrbitalSliceGrouped}) =
-    maybe_OrbitalSliceVector(x, scalarize(parent(i)))
-maybe_OrbitalSliceVector(x, i::OrbitalSliceGrouped) = OrbitalSliceVector(x, (i,))
-
-# fallback
-maybe_OrbitalSliceVector(x, i) = x
-
-maybe_OrbitalSliceMatrix(x, i::OrbitalSliceGrouped) =
-    maybe_OrbitalSliceMatrix(x, (i, i))
-maybe_OrbitalSliceMatrix(x, (i, j)::Tuple{OrbitalSliceGrouped,OrbitalSliceGrouped}) =
-    OrbitalSliceMatrix(x, (i, j))
-
-# fallback
-maybe_OrbitalSliceMatrix(x, i) = x
 
 #endregion
diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl
index c40236d9..2d330ac9 100644
--- a/test/test_greenfunction.jl
+++ b/test/test_greenfunction.jl
@@ -55,8 +55,8 @@ end
     s1´ = GS.Schur(boundary = -1)
     sites´ = (; region = r -> abs(r[2]) < 2 && r[1] == 0)
     # non-hermitian Σ model
-    mod = @onsite((ω, r; o = 1) -> (o - im*ω)*I) +
-          plusadjoint(@onsite((ω; p=1)-> p*I) +  @hopping((ω, r, dr; t = 1) -> im*dr[1]*t*I; range = 1))
+    mod = @onsite((ω, r; o = 1) -> (o - im*ω)*I) + @onsite((ω, s; o = 1, b = 2) --> o*b*pos(s)[1]*I) +
+          plusadjoint(@onsite((ω; p=1)-> p*I) + @hopping((ω, r, dr; t = 1) -> im*dr[1]*t*I; range = 1))
     g0, g0´, g1´ = greenfunction(h0, s0), greenfunction(h0, s0´), greenfunction(h1, s1´)
     for (h, s) in zip((h0, h0, h1, h1), (s0, s0´, s1, s1´))
         oh = h |> attach(nothing; sites´...)
@@ -329,6 +329,11 @@ end
     @test ρ0[cellsites(SA[0], 1), cellsites(SA[1], 1)] isa Matrix
     @test size(view(ρ0, cellsites(SA[0], 1), cellsites(SA[1], 1))) == (2, 2)
 
+    # Diagonal slicing not yet supported
+    @test_broken densitymatrix(g1[diagonal(cells = SA[1])], 5)
+    @test_broken densitymatrix(gSpectrum[diagonal(cells = SA[])])
+    @test_broken densitymatrix(gKPM[diagonal(1)])
+
     glead = LP.square() |> hamiltonian(hopping(1)) |> supercell((0,1), region = r -> -1 <= r[1] <= 1) |> attach(nothing; cells = SA[10]) |> greenfunction(GS.Schur(boundary = 0));
     contact1 = r -> r[1] ≈ 5 && -1 <= r[2] <= 1
     contact2 = r -> r[2] ≈ 5 && -1 <= r[1] <= 1
@@ -346,3 +351,18 @@ end
     testcond(g0; nambu = true)
     testjosephson(g0)
 end
+
+@testset "mean-field models" begin
+    h1 = LP.honeycomb() |> supercell(2) |> supercell |> hamiltonian(onsite(0I) + hopping(I), orbitals = 1)
+    h2 = LP.honeycomb() |> supercell(2) |> supercell |> hamiltonian(onsite(0I) + hopping(I), orbitals = 2)
+    h3 = LP.honeycomb() |> supercell(2) |> supercell |> hamiltonian(onsite(0I) + hopping(I), orbitals = (1,2))
+    for h0 in (h1, h2, h3)
+        ρ0 = densitymatrix(greenfunction(h0, GS.Spectrum())[cells = SA[]])();
+        h = h0 |> @onsite!((o, s; ρ = ρ0, t) --> o + t*ρ[s])
+        @test diag(h(t = 2)[()]) ≈ 2*diag(ρ0) atol = 0.0000001
+        h = h0 |> @hopping!((t, si, sj; ρ = ρ0, α = 2) --> α*ρ[si, sj])
+        @test h() isa Quantica.Hamiltonian
+        diff = (h()[()] - 2ρ0) .* h()[()]
+        @test iszero(diff)
+    end
+end
diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl
index bd1e2853..9d441ec0 100644
--- a/test/test_hamiltonian.jl
+++ b/test/test_hamiltonian.jl
@@ -117,8 +117,17 @@ end
 end
 
 @testset "models" begin
-    mo = (onsite(1), onsite(r-> r[1]), @onsite((; o) -> o), @onsite((r; o) -> r[1]*o))
-    mh = (hopping(1), hopping((r, dr)-> im*dr[1]), @hopping((; t) -> t), @hopping((r, dr; t) -> r[1]*t))
+    mo = (onsite(1), onsite(r-> r[1]), @onsite((; o) -> o), @onsite((r; o=2) -> r[1]*o),
+         @onsite((s; o, p) --> pos(s)[1]*o))
+    mh = (hopping(1), hopping((r, dr)-> im*dr[1]), @hopping((; t) -> t), @hopping((r, dr; t) -> r[1]*t),
+        @hopping((si, sj) --> im*ind(si)), @hopping((si, sj; t, p = 2) --> pos(sj)[1]*t))
+    argso, argsh = (0, 1, 0, 1, 1), (0, 2, 0, 2, 2, 2)
+    for (o, no) in zip(mo, argso)
+        @test Quantica.narguments(only(Quantica.terms(o))) == no
+    end
+    for (h, nh) in zip(mh, argsh)
+        @test Quantica.narguments(only(Quantica.terms(h))) == nh
+    end
     for o in mo, h in mh
         @test length(Quantica.allterms(-o - 2*h)) == 2
         @test Quantica.ParametricModel(o+h) isa Quantica.ParametricModel
@@ -329,8 +338,13 @@ end
     h0 = LP.square() |> hopping(1) |> supercell(3) |> @hopping!((t, r, dr; A = SA[1,2]) -> t*cis(A'dr))
     h = torus(h0, (0.2,:))
     @test h0((0.2, 0.3)) ≈ h((0.3,))
+    # non-spatial models
+    h = LP.linear() |> @hopping((i,j) --> ind(i) + ind(j)) + @onsite((i; k = 1) --> pos(i)[k])
+    @test ishermitian(h())
+
 end
 
+
 @testset "hamiltonian nrange" begin
     lat = LatticePresets.honeycomb(a0 = 2)
     @test Quantica.nrange(1, lat) ≈ 2/√3
diff --git a/test/test_show.jl b/test/test_show.jl
index a9886274..ddabf567 100644
--- a/test/test_show.jl
+++ b/test/test_show.jl
@@ -2,7 +2,7 @@
     hs = HP.graphene(orbitals = 2), HP.graphene(orbitals = (2,1))
     for h in hs
         b = bands(h, subdiv(0,2pi,10), subdiv(0,2pi,10))
-        g = greenfunction(supercell(h) |> attach(@onsite(ω -> im*I)) |> attach(nothing))
+        g = greenfunction(supercell(h) |> attach(@onsite(ω -> im*I)) |> attach(nothing), GS.Spectrum())
         @test nothing === show(stdout, sublat((0,0)))
         @test nothing === show(stdout, LP.honeycomb())
         @test nothing === show(stdout, LP.honeycomb()[cells = (0,0)])
@@ -24,6 +24,7 @@
         @test nothing === show(stdout, Quantica.slice(b, (0,0)))
         @test nothing === show(stdout, g)
         @test nothing === show(stdout, g[cells = ()])
+        @test nothing === show(stdout, MIME("text/plain"), g[diagonal(cells = ())](0.1))
         @test nothing === show(stdout, g(0.1))
         @test nothing === show(stdout, ldos(g[1]))
         @test nothing === show(stdout, ldos(g(0.1)))
@@ -31,6 +32,8 @@
         @test nothing === show(stdout, current(g(0.1)))
         @test nothing === show(stdout, conductance(g[1]))
         @test nothing === show(stdout, transmission(g[1,2]))
+        @test nothing === show(stdout, densitymatrix(g[1]))
+        @test nothing === show(stdout, MIME("text/plain"), densitymatrix(g[1])())
     end
     h = first(hs)
     g = greenfunction(supercell(h) |> attach(@onsite(ω -> im*I)) |> attach(nothing))