Introducing the current Operator

src/docstrings.jl

@@ -577,7 +577,7 @@ Alternative and less general curried form equivalent to `hamiltonian(h, modifier
 
 Return the Bloch harmonic of an `h::AbstractHamiltonian` in the form of a `SparseMatrixCSC`
 with complex scalar `eltype`. This matrix is "flat", in the sense that it contains matrix
-elements between indivisual orbitals, not sites. This distinction is only relevant for
+elements between individual orbitals, not sites. This distinction is only relevant for
 multiorbital Hamiltonians. To access the non-flattened matrix use `h[unflat(dn)]` (see
 also `unflat`).
 
@@ -1636,6 +1636,13 @@ true
 ldos
 
 """
+    current(h::AbstractHamiltonian; charge = -I, direction = 1)
+
+Build an `Operator` object that behaves like a `ParametricHamiltonian` in regards to calls
+and getindex, but whose matrix elements are hoppings ``im*(rⱼ-rᵢ)[direction]*charge*tⱼᵢ``,
+where `tᵢⱼ` are the hoppings in `h`. This operator is equal to ``∂h/∂Aᵢ``, where `Aᵢ`
+is a gauge field along `direction = i`.
+
     current(gs::GreenSlice; charge = -I, direction = missing)
 
 Build `Js::CurrentDensitySlice`, a partially evaluated object representing the equilibrium

src/hamiltonian.jl

@@ -440,7 +440,12 @@ Base.getindex(h::AbstractHamiltonian, dn::HybridInds{<:Union{Integer,Tuple}}) =
 Base.getindex(h::AbstractHamiltonian{<:Any,<:Any,L}, ::HybridInds{Tuple{}}) where {L} =
     h[hybrid(zero(SVector{L,Int}))]
 
-function Base.getindex(h::AbstractHamiltonian{<:Any,<:Any,L}, dn::HybridInds{SVector{L,Int}}) where {L}
+function Base.getindex(h::ParametricHamiltonian{<:Any,<:Any,L}, dn::HybridInds{SVector{L,Int}}; params...) where {L}
+    h´ = call!(h; params...)
+    return getindex(h´, dn)
+end
+
+function Base.getindex(h::Hamiltonian{<:Any,<:Any,L}, dn::HybridInds{SVector{L,Int}}) where {L}
     for har in harmonics(h)
         parent(dn) == dcell(har) && return matrix(har)
     end

src/observables.jl

@@ -50,6 +50,19 @@ check_different_contact_slice(gs) = (slicerows(gs) isa Integer && slicecols(gs)
 
 #endregion
 
+############################################################################################
+# Operators
+#region
+
+function current(h::AbstractHamiltonian; charge = -I, direction = 1)
+    current = parametric(h,
+        @onsite!(o -> zero(o)),
+        @hopping!((t, r, dr) -> im*dr[direction]*charge*t))
+    return Operator(current)
+end
+
+#endregion
+
 ############################################################################################
 # Integrator - integrates a function f along a complex path ωcomplex(ω::Real), connecting ωi
 #   The path is piecewise linear in the form of a sawtooth with a given ± slope

src/show.jl

@@ -152,25 +152,32 @@ Base.summary(::HoppingModifier{N}) where {N} = "HoppingModifier{ParametricFuncti
 # Hamiltonian
 #region
 
-function Base.show(io::IO, h::Union{Hamiltonian,ParametricHamiltonian})
+function Base.show(io::IO, h::Union{Hamiltonian,ParametricHamiltonian,Operator})
     i = get(io, :indent, "")
-    print(io, i, summary(h), "\n",
+    print(io, i, summary(h), "\n", showstring(h, i))
+    showextrainfo(io, i, h)
+end
+
+showstring(h::Union{Hamiltonian,ParametricHamiltonian}, i) =
 "$i  Bloch harmonics  : $(length(harmonics(h)))
 $i  Harmonic size    : $((n -> "$n × $n")(size(h, 1)))
 $i  Orbitals         : $(norbitals(h))
 $i  Element type     : $(displaytype(blocktype(h)))
 $i  Onsites          : $(nonsites(h))
 $i  Hoppings         : $(nhoppings(h))
-$i  Coordination     : $(round(coordination(h), digits = 5))")
-    showextrainfo(io, i, h)
-end
+$i  Coordination     : $(round(coordination(h), digits = 5))"
+
+showstring(o::Operator, i) = showstring(hamiltonian(o), i)
 
 Base.summary(h::Hamiltonian{T,E,L}) where {T,E,L} =
     "Hamiltonian{$T,$E,$L}: Hamiltonian on a $(L)D Lattice in $(E)D space"
 
 Base.summary(h::ParametricHamiltonian{T,E,L}) where {T,E,L} =
     "ParametricHamiltonian{$T,$E,$L}: Parametric Hamiltonian on a $(L)D Lattice in $(E)D space"
 
+Base.summary(h::Operator{H}) where {T,E,L,H<:AbstractHamiltonian{T,E,L}} =
+    "Operator{$T,$E,$L}: Operator on a $(L)D Lattice in $(E)D space"
+
 displaytype(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = "$N × $N blocks ($T)"
 displaytype(::Type{S}) where {N,T,S<:SMatrixView{N,N,T}} = "At most $N × $N blocks ($T)"
 displaytype(::Type{T}) where {T} = "scalar ($T)"
@@ -181,6 +188,7 @@ showextrainfo(io, i, h) = nothing
 showextrainfo(io, i, h::ParametricHamiltonian) = print(io, i, "\n",
 "$i  Parameters       : $(parameters(h))")
 
+showextrainfo(io, i, o::Operator) = showextrainfo(io, i, hamiltonian(o))
 #endregion
 
 ############################################################################################

src/types.jl

@@ -1990,3 +1990,26 @@ Base.:(==)(g::GreenSlice, g´::GreenSlice) = function_not_defined("==")
 
 #endregion
 #endregion
+
+############################################################################################
+# Operator - Hamiltonian-like operator representing observables other than a Hamiltonian
+#   see observables.jl for specific constructors
+#region
+
+struct Operator{H<:AbstractHamiltonian}
+    h::H
+end
+
+#region ## API ##
+
+hamiltonian(o::Operator) = o.h
+
+(c::Operator)(φ...; kw...) = c.h(φ...; kw...)
+
+call!(c::Operator, φ...; kw...) = call!(c.h, φ...; kw...)
+
+Base.getindex(c::Operator, i...) = getindex(c.h, i...)
+
+#endregion
+
+#endregion

test/test_hamiltonian.jl

@@ -1,4 +1,4 @@
-using Quantica: Hamiltonian, ParametricHamiltonian, sites, nsites, nonsites, nhoppings, coordination, flat, hybrid, transform!
+using Quantica: Hamiltonian, ParametricHamiltonian, sites, nsites, nonsites, nhoppings, coordination, flat, hybrid, transform!, nnz, nonzeros
 
 @testset "basic hamiltonians" begin
     presets = (LatticePresets.linear, LatticePresets.square, LatticePresets.triangular, LatticePresets.honeycomb,
@@ -376,3 +376,19 @@ end
     h = combine(hb, h0, ht; coupling = hopping((r,dr) -> exp(-norm(dr)), range = 2))
     @test !iszero(h((0,0))[1:2, 5:6])
 end
+
+
+@testset "current operator" begin
+    h = LP.honeycomb() |> hamiltonian(@onsite((; μ = 0) -> (2-μ)*I) - hopping(SA[0 1; 1 0]), orbitals = 2) |> supercell(2)
+    co = current(h, direction = 2)
+    c = co(SA[0,0])
+    @test c ≈ c'
+    @test iszero(diag(c))
+    @test all(x -> real(x) ≈ 0, c)
+    cp = co[unflat(SA[1,0])]
+    cm = co[unflat(SA[-1,0])]
+    @test nnz(cp) == nnz(cm) == 2
+    @test cp ≈ cm'
+    @test all(x -> x[1] ≈ x[2]', zip(nonzeros(cp), nonzeros(cm)))
+    @test all(x -> iszero(real(x)), nonzeros(cp))
+end

test/test_show.jl

@@ -1,33 +1,38 @@
 @testset "show methods" begin
-    h = HP.graphene(orbitals = 2)
-    b = bands(h, subdiv(0,2pi,10), subdiv(0,2pi,10))
+    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))
+        @test nothing === show(stdout, sublat((0,0)))
+        @test nothing === show(stdout, LP.honeycomb())
+        @test nothing === show(stdout, LP.honeycomb()[cells = (0,0)])
+        @test nothing === show(stdout, siteselector(; cells = (0,0)))
+        @test nothing === show(stdout, hopselector(; range = 1))
+        @test nothing === show(stdout, onsite(1) + hopping(1))
+        @test nothing === show(stdout, @onsite(()->1) + @hopping(()->1))
+        @test nothing === show(stdout, @onsite!(o->1))
+        @test nothing === show(stdout, @hopping!(t->1))
+        @test nothing === show(stdout, h)
+        @test nothing === show(stdout, h |> hamiltonian(@onsite!(o->2o)))
+        @test nothing === show(stdout, h |> attach(nothing, cells = (0,0)))
+        @test nothing === show(stdout, current(h))
+        @test nothing === show(stdout, ES.LinearAlgebra())
+        @test nothing === show(stdout, spectrum(h, (0,0)))
+        @test nothing === show(stdout, b)
+        @test nothing === show(stdout, b[(0,0)])
+        @test nothing === show(stdout, Quantica.slice(b, (0,0)))
+        @test nothing === show(stdout, Quantica.slice(b, (0,0)))
+        @test nothing === show(stdout, g)
+        @test nothing === show(stdout, g[cells = (0,0)])
+        @test nothing === show(stdout, g(0.1))
+        @test nothing === show(stdout, ldos(g[1]))
+        @test nothing === show(stdout, ldos(g(0.1)))
+        @test nothing === show(stdout, current(g[1]))
+        @test nothing === show(stdout, current(g(0.1)))
+        @test nothing === show(stdout, conductance(g[1]))
+        @test nothing === show(stdout, transmission(g[1,2]))
+    end
+    h = first(hs)
     g = greenfunction(supercell(h) |> attach(@onsite(ω -> im*I)) |> attach(nothing))
-    @test nothing === show(stdout, sublat((0,0)))
-    @test nothing === show(stdout, LP.honeycomb())
-    @test nothing === show(stdout, LP.honeycomb()[cells = (0,0)])
-    @test nothing === show(stdout, siteselector(; cells = (0,0)))
-    @test nothing === show(stdout, hopselector(; range = 1))
-    @test nothing === show(stdout, onsite(1) + hopping(1))
-    @test nothing === show(stdout, @onsite(()->1) + @hopping(()->1))
-    @test nothing === show(stdout, @onsite!(o->1))
-    @test nothing === show(stdout, @hopping!(t->1))
-    @test nothing === show(stdout, h)
-    @test nothing === show(stdout, h |> hamiltonian(@onsite!(o->2o)))
-    @test nothing === show(stdout, h |> attach(nothing, cells = (0,0)))
-    @test nothing === show(stdout, ES.LinearAlgebra())
-    @test nothing === show(stdout, spectrum(h, (0,0)))
-    @test nothing === show(stdout, b)
-    @test nothing === show(stdout, b[(0,0)])
-    @test nothing === show(stdout, Quantica.slice(b, (0,0)))
-    @test nothing === show(stdout, Quantica.slice(b, (0,0)))
-    @test nothing === show(stdout, g)
-    @test nothing === show(stdout, g[cells = (0,0)])
-    @test nothing === show(stdout, g(0.1))
-    @test nothing === show(stdout, ldos(g[1]))
-    @test nothing === show(stdout, ldos(g(0.1)))
-    @test nothing === show(stdout, current(g[1]))
-    @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, josephson(g[1], 2))
-end
+end