Merge pull request #273 from pablosanjose/qplot_openham

Add `qplot(::OpenHamiltonian)`

ext/QuanticaMakieExt/QuanticaMakieExt.jl

@@ -4,10 +4,10 @@ using Makie
 using Quantica
 using Makie.GeometryBasics
 using Makie.GeometryBasics: Ngon
-using Quantica: Lattice, LatticeSlice, AbstractHamiltonian, Hamiltonian,
+using Quantica: Lattice, LatticeSlice, AbstractHamiltonian, Hamiltonian, OpenHamiltonian,
       ParametricHamiltonian, Harmonic, Bravais, SVector, GreenFunction, GreenSolution,
       argerror, harmonics, sublats, siterange, site, norm,
-      normalize, nsites, nzrange, rowvals, nonzeros, sanitize_SVector
+      normalize, nsites, nzrange, rowvals, nonzeros, sanitize_SVector, default_hamiltonian
 
 import Quantica: plotlattice, plotlattice!, plotbands, plotbands!, qplot, qplot!, qplotdefaults
 
@@ -17,7 +17,8 @@ const PlotLatticeArgumentType{E} = Union{
     Lattice{<:Any,E},
     LatticeSlice{<:Any,E},
     AbstractHamiltonian{<:Any,E},
-    GreenFunction{<:Any,E}}
+    GreenFunction{<:Any,E},
+    OpenHamiltonian{<:Any,E}}
 
 const PlotBandsArgumentType{E} =
     Union{Quantica.Bandstructure{<:Any,E},

ext/QuanticaMakieExt/docstrings.jl

@@ -6,14 +6,14 @@ specialized plot recipes:
 
 - `object::Lattice`             -> `plotlattice` (supports also `LatticeSlice`s)
 - `object::Hamiltonian`         -> `plotlattice` (supports also `ParametricHamiltonian`)
-- `object::GreenFunction`       -> `plotlattice` (see below)
+- `object::GreenFunction`       -> `plotlattice` (supports also `OpenHamiltonian`, see below)
 - `object::Bandstructure`       -> `plotbands`   (supports also slices of `Bandstructure`)
 - `object::Subband`             -> `plotbands`   (supports also collections of `Subbands`)
 
 Supported `Makie` backends include `GLMakie`, `CairoMakie`, `WGLMakie`, etc. Instead of
 `using Makie`, load a specific backend directly, e.g. `using GLMakie`.
 
-    qplot(g::GreenFunction; children = (plotkw₁::NamedTuple, plotkw₂::NamedTuple, ...), kw...)
+    qplot(g::Union{GreenFunction,OpenHamiltonian}; children = (plotkw₁::NamedTuple, plotkw₂::NamedTuple, ...), kw...)
 
 Render the parent Hamiltonian of `g`, along with a representation of each of `g`'s
 self-energies. Specific `plotkwᵢ` keywords for each self-energy can be given in `children`,

ext/QuanticaMakieExt/plotlattice.jl

@@ -494,11 +494,11 @@ function Makie.plot!(plot::PlotLattice{Tuple{H,S,S´}}) where {H<:Hamiltonian,S<
     return plot
 end
 
-function Makie.plot!(plot::PlotLattice{Tuple{G}}) where {G<:GreenFunction}
+function Makie.plot!(plot::PlotLattice{Tuple{G}}) where {G<:Union{GreenFunction,OpenHamiltonian}}
     g = to_value(plot[1])
     gsel = haskey(plot, :selector) && plot[:selector][] !== missing ?
         plot[:selector][] : green_selector(g)
-    h = hamiltonian(g)
+    h = default_hamiltonian(g)
     latslice = lattice(h)[gsel]
     latslice´ = Quantica.growdiff(latslice, h)
     L = Quantica.latdim(h)
@@ -528,7 +528,7 @@ function Makie.plot!(plot::PlotLattice{Tuple{G}}) where {G<:GreenFunction}
     # plot contacts
     if !ishidden((:contact, :contacts), plot)
         Σkws = Iterators.cycle(parse_children(plot[:children]))
-        Σs = Quantica.selfenergies(Quantica.contacts(g))
+        Σs = Quantica.selfenergies(g)
         hideΣ = Quantica.tupleflatten(:bravais, plot[:hide][])
         for (Σ, Σkw) in zip(Σs, Σkws)
             Σplottables = Quantica.selfenergy_plottables(Σ)
@@ -575,7 +575,7 @@ boundary_selector(g) = siteselector(cells = n -> isboundarycell(n, g))
 
 function green_bounding_box(g)
     L = Quantica.latdim(hamiltonian(g))
-    return isempty(Quantica.contacts(g)) ?
+    return isempty(Quantica.selfenergies(g)) ?
         boundary_bounding_box(Val(L), Quantica.boundaries(g)...) :
         broaden_bounding_box(Quantica.boundingbox(g), Quantica.boundaries(g)...)
 end
@@ -748,9 +748,9 @@ Makie.convert_arguments(::PointBased, lat::Lattice, sublat = missing) =
     (Point.(Quantica.sites(lat, sublat)),)
 Makie.convert_arguments(p::PointBased, lat::Lattice{<:Any,1}, sublat = missing) =
     Makie.convert_arguments(p, Quantica.lattice(lat, dim = Val(2)), sublat)
-Makie.convert_arguments(p::PointBased, h::Union{AbstractHamiltonian,GreenFunction,GreenSolution}, sublat = missing) =
+Makie.convert_arguments(p::PointBased, h::Union{AbstractHamiltonian,GreenFunction,GreenSolution,OpenHamiltonian}, sublat = missing) =
     Makie.convert_arguments(p, Quantica.lattice(h), sublat)
-Makie.convert_arguments(p::Type{<:LineSegments}, g::Union{GreenFunction,GreenSolution}, sublat = missing) =
+Makie.convert_arguments(p::Type{<:LineSegments}, g::Union{GreenFunction,GreenSolution,OpenHamiltonian}, sublat = missing) =
     Makie.convert_arguments(p, Quantica.hamiltonian(g), sublat)
 
 function Makie.convert_arguments(::Type{<:LineSegments}, h::AbstractHamiltonian{<:Any,E}, sublat = missing) where {E}

src/greenfunction.jl

@@ -21,14 +21,25 @@ default_green_solver(::AbstractHamiltonian) = GS.Bands()
 
 #endregion
 
+############################################################################################
+# Contacts call! API
+#region
+
+function call!(c::Contacts, ω; params...)
+    Σblocks = selfenergyblocks(c)
+    call!.(c.selfenergies, Ref(ω); params...) # updates matrices in Σblocks
+    return Σblocks
+end
+
+call!_output(c::Contacts) = selfenergyblocks(c)
+
+#endregion
+
 ############################################################################################
 # GreenFuntion call! API
 #region
 
-## TODO: test copy(g) for potential aliasing problems
-(g::GreenFunction)(; params...) = minimal_callsafe_copy(call!(g; params...))
 (g::GreenFunction)(ω; params...) = minimal_callsafe_copy(call!(g, ω; params...))
-(g::GreenSlice)(; params...) = minimal_callsafe_copy(call!(g; params...))
 (g::GreenSlice)(ω; params...) = copy(call!(g, ω; params...))
 
 call!(g::G, ω; params...) where {T,G<:Union{GreenFunction{T},GreenSlice{T}}} =
@@ -49,9 +60,6 @@ function call!(g::GreenFunction{T}, ω::Complex{T}; params...) where {T}
     return GreenSolution(g, slicer, Σblocks, corbs)
 end
 
-call!(g::GreenSlice; params...) =
-    GreenSlice(call!(greenfunction(g); params...), greenindices(g)...)
-
 call!(g::GreenSlice{T}, ω::T; params...) where {T} =
     call!(g, retarded_omega(ω, solver(parent(g))); params...)
 
@@ -70,18 +78,43 @@ needs_omega_shift(s::AppliedGreenSolver) = true
 #endregion
 
 ############################################################################################
-# Contacts call! API
+# DeparametrizedGreenSolver
+#   support for g(; params...) --> GreenFunction (not a wrapper, completely independent)
+#   DeparametrizedGreenSolver doesn't need to implement the AppliedGreenSolver API, since it
+#   forwards to its parent
 #region
 
-call!(c::Contacts; params...) = Contacts(call!.(c.selfenergies; params...), c.orbitals, c.orbslice)
+struct DeparametrizedGreenSolver{P,G<:GreenFunction} <: AppliedGreenSolver
+    gparent::G
+    params::P
+end
 
-function call!(c::Contacts, ω; params...)
-    Σblocks = selfenergyblocks(c)
-    call!.(c.selfenergies, Ref(ω); params...) # updates matrices in Σblocks
-    return Σblocks
+Base.parent(s::DeparametrizedGreenSolver) = s.gparent
+parameters(s::DeparametrizedGreenSolver) = s.params
+
+function (g::GreenFunction)(; params...)
+    h´ = minimal_callsafe_copy(parent(g))
+    c´ = minimal_callsafe_copy(contacts(g))
+    s´ = minimal_callsafe_copy(solver(g), h´, c´)
+    gparent = GreenFunction(h´, s´, c´)
+    return GreenFunction(h´, DeparametrizedGreenSolver(gparent, params), c´)
 end
 
-call!_output(c::Contacts) = selfenergyblocks(c)
+# params are ignored, solver.params are used instead. T required to disambiguate.
+function call!(g::GreenFunction{T,<:Any,<:Any,<:DeparametrizedGreenSolver}, ω::Complex{T}; params...) where {T}
+    s = solver(g)
+    return call!(s.gparent, ω; s.params...)
+end
+
+function minimal_callsafe_copy(s::DeparametrizedGreenSolver, parentham, parentcontacts)
+    solver´ = minimal_callsafe_copy(solver(s.gparent), parentham, parentcontacts)
+    gparent = GreenFunction(parentham, solver´, parentcontacts)
+    s´ = DeparametrizedGreenSolver(gparent, s.params)
+    return s´
+end
+
+default_hamiltonian(g::GreenFunction{<:Any,<:Any,<:Any,<:Quantica.DeparametrizedGreenSolver}) =
+    default_hamiltonian(parent(g); parameters(solver(g))...)
 
 #endregion
 

src/types.jl

@@ -1404,6 +1404,9 @@ function harmonic_index(h::AbstractHamiltonian, dn)
     return first(harmonics(h)), 1  # unreachable
 end
 
+# Unless params are given, it returns the Hamiltonian with defaults parameters
+default_hamiltonian(h::AbstractHamiltonian; params...) = h(; params...)
+
 ## Hamiltonian
 
 Hamiltonian(l::Lattice{T,E,L}, b::OrbitalBlockStructure{B}, h::Vector{Harmonic{T,L,B}}, bl) where {T,E,L,B} =
@@ -1769,34 +1772,13 @@ Base.length(b::Bandstructure) = length(b.subbands)
 #endregion
 
 ############################################################################################
-# SelfEnergy solvers - see solvers/selfenergy.jl for self-energy solvers
+# SelfEnergySolvers - see solvers/selfenergy.jl for self-energy solvers
 #region
 
 abstract type AbstractSelfEnergySolver end
 abstract type RegularSelfEnergySolver <: AbstractSelfEnergySolver end
 abstract type ExtendedSelfEnergySolver <: AbstractSelfEnergySolver end
 
-# Support for call!(s::AbstractSelfEnergySolver; params...) -> AbstractSelfEnergySolver
-## TODO: revisit to see if there is a better/simpler approach
-
-struct WrappedRegularSelfEnergySolver{F} <: RegularSelfEnergySolver
-    f::F
-end
-
-struct WrappedExtendedSelfEnergySolver{F} <: ExtendedSelfEnergySolver
-    f::F
-end
-
-call!(s::WrappedRegularSelfEnergySolver, ω; params...) = s.f(ω)
-call!(s::WrappedExtendedSelfEnergySolver, ω; params...) = s.f(ω)
-
-call!(s::RegularSelfEnergySolver; params...) =
-    WrappedRegularSelfEnergySolver(ω -> call!(s, ω; params...))
-call!(s::ExtendedSelfEnergySolver; params...) =
-    WrappedExtendedSelfEnergySolver(ω -> call!(s, ω; params...))
-call!(s::WrappedRegularSelfEnergySolver; params...) = s
-call!(s::WrappedExtendedSelfEnergySolver; params...) = s
-
 #endregion
 
 ############################################################################################
@@ -1852,11 +1834,12 @@ has_selfenergy(s::SelfEnergy) = has_selfenergy(solver(s))
 has_selfenergy(s::AbstractSelfEnergySolver) = true
 # see nothing.jl for override for the case of SelfEnergyEmptySolver
 
-call!(Σ::SelfEnergy; params...) = SelfEnergy(call!(Σ.solver; params...), Σ.orbslice)
 call!(Σ::SelfEnergy, ω; params...) = call!(Σ.solver, ω; params...)
 
 call!_output(Σ::SelfEnergy) = call!_output(solver(Σ))
 
+(Σ::SelfEnergy)(; params...) = SelfEnergy(Σ.solver(; params...), Σ.orbslice)
+
 minimal_callsafe_copy(Σ::SelfEnergy) =
     SelfEnergy(minimal_callsafe_copy(Σ.solver), Σ.orbslice)
 
@@ -1885,6 +1868,8 @@ selfenergies(oh::OpenHamiltonian) = oh.selfenergies
 
 hamiltonian(oh::OpenHamiltonian) = oh.h
 
+default_hamiltonian(oh::OpenHamiltonian) = default_hamiltonian(oh.h)
+
 lattice(oh::OpenHamiltonian) = lattice(oh.h)
 
 zerocell(h::OpenHamiltonian) = zerocell(parent(h))
@@ -1908,6 +1893,9 @@ Base.size(oh::OpenHamiltonian, i...) = size(oh.h, i...)
 
 Base.parent(oh::OpenHamiltonian) = oh.h
 
+boundingbox(oh::OpenHamiltonian) =
+    boundingbox(tupleflatten(boundingbox.(orbslice.(selfenergies(oh)))...))
+
 #endregion
 #endregion
 
@@ -2089,8 +2077,13 @@ Base.parent(i::DiagIndices) = i.inds
 
 kernel(i::DiagIndices) = i.kernel
 
+# returns the Hamiltonian field
 hamiltonian(g::Union{GreenFunction,GreenSolution,GreenSlice}) = hamiltonian(g.parent)
 
+# Like the above, but it may not be === the field (it can be a copy with parameters applied)
+# needed for qplot(g(; params...))
+default_hamiltonian(g::GreenFunction) = default_hamiltonian(parent(g))  # default params
+
 lattice(g::Union{GreenFunction,GreenSolution,GreenSlice}) = lattice(g.parent)
 
 latslice(g::GreenFunction, i) = orbslice(g.contacts, i)
@@ -2114,6 +2107,7 @@ ncontacts(g::Union{GreenSolution,GreenSlice}) = ncontacts(parent(g))
 
 slicer(g::GreenSolution) = g.slicer
 
+selfenergies(g::GreenFunction) = selfenergies(contacts(g))
 selfenergies(g::GreenSolution) = g.contactΣs
 
 has_selfenergy(g::Union{GreenFunction,GreenSlice,GreenSolution}) =
@@ -2168,8 +2162,8 @@ function similar_Matrix(gs::GreenSlice{T}) where {T}
 end
 
 boundaries(g::GreenFunction) = boundaries(solver(g))
-# fallback
-boundaries(g::AppliedGreenSolver) = ()
+# fallback (for solvers without boundaries, or for OpenHamiltonian)
+boundaries(_) = ()
 
 boundingbox(g::GreenFunction) = isempty(contacts(g)) ?
     (zerocell(lattice(g)), zerocell(lattice(g))) : boundingbox(contacts(g))

test/test_greenfunction.jl

@@ -1,4 +1,5 @@
-using Quantica: GreenFunction, GreenSlice, GreenSolution, zerocell, CellOrbitals, ncontacts
+using Quantica: GreenFunction, GreenSlice, GreenSolution, zerocell, CellOrbitals, ncontacts,
+    solver
 
 function testgreen(h, s; kw...)
     ω = 0.2
@@ -110,6 +111,17 @@ end
     @test (@allocations view(gs, 1, 2)) <= 2  # should be 1, but in some platforms/versions it could be 2
 end
 
+@testset "GreenFunction partial evaluation" begin
+    g = LP.linear() |> hamiltonian(@hopping((; q = 1) -> q*I), orbitals = 2) |> attach(@onsite((ω; p = 0) ->p*SA[0 1; 1 0]), cells = 1) |> greenfunction
+    g´ = g(p = 2)
+    h´, h = hamiltonian(g´), hamiltonian(g)
+    @test h´ !== h
+    @test h´ == h(p = 2)
+    @test g(0.2, p = 2)[] == g´(0.2)[]
+    @test g(0.2, p = 1)[] != g´(0.2, p = 1)[]       # p = 1 is ignored by g´, fixed to p = 2
+    @test solver(g) !== solver(parent(solver(g´)))  # no aliasing
+end
+
 @testset "GreenSolvers applicability" begin
     h = HP.graphene()
     @test_throws ArgumentError greenfunction(h, GS.Schur())

test/test_plots.jl

@@ -46,6 +46,13 @@ end
         attach(nothing, region = RP.circle(1, SA[2,3])) |> attach(nothing, region = RP.circle(1, SA[3,-3])) |>
         greenfunction(GS.Bands(subdiv(-π, π, 13), subdiv(-π, π, 13), boundary = 2=>-3))
     @test qplot(g) isa Figure
+    # Issue 200
+    g = LP.linear() |> hamiltonian(@hopping((; q = 1) -> q*I), orbitals = 2) |> attach(@onsite((ω; p = 0) ->p*SA[0 1; 1 0]), cells = 1) |> greenfunction
+    @test qplot(g) isa Figure
+    @test qplot(g(p = 3)) isa Figure
+    # Issue 243
+    oh = LP.linear() |> hopping(1) |> attach(@onsite((ω; p = 1) -> p), cells = 1) |> attach(@onsite((ω; p = 1) -> p), cells = 3)
+    @test qplot(oh) isa Figure
 end
 
 @testset "plot bands" begin