Alternative representation of GreenFunction boundaries in plots (#255)

* distinguish between null and unspecified selectors

more fastpaths and tests

* try boundary as cell overlay

* boundaries plotted as cells

* fix 1D and green_boundingbox, add boxes over boundary sites,

docs/src/tutorial/greenfunctions.md

@@ -191,6 +191,8 @@ julia> qplot(g, children = (; selector = siteselector(; cells = 1:5), sitecolor
 
 Note that since we did not specify the `solver` in `greenfunction`, the `L=0` default `GS.SparseLU()` was taken.
 
+We can also visualize `glead`, which is defined on a 1D lattice with a boundary. Boundary cells are shown by default in red
+
 !!! tip "The GreenFunction <-> AbstractHamiltonian relation"
     Its important un appreciate that a `g::GreenFunction` represents the retarded Green function between sites in a given `AbstractHamiltonian`, but not on sites of the coupled `AbstractHamiltonians` of its attached self-energies. Therefore, `gcentral` above cannot yield observables in the leads (blue sites above), only on the red sites. To obtain observables in a given lead, its `GreenFunction` must be constructed, with an attached self-energy coming from the central region plus the other three leads.
 

docs/src/tutorial/observables.md

@@ -247,7 +247,7 @@ Integrator: Complex-plane integrator
     Contact                 : 1
     Number of phase shifts  : 0
 
-julia> qplot(g)
+julia> qplot(g, children = (; sitecolor = :blue))
 ```
 ```@raw html
 <img src="../../assets/josephson_lat.png" alt="Josephson junction" width="400" class="center"/>

ext/QuanticaMakieExt/docstrings.jl

@@ -56,13 +56,20 @@ Render a finite subset of sites in a lattice and its Bravais vectors.
 Render the lattice on which `h` is defined, together with all hoppings in the form of
 inter-site links.
 
+    plotlattice(h::GreenFunction; kw...)
+
+Render the lattice on which `h` is defined, together with all hoppings in the form of
+inter-site links, and a representation of contacts.
+
 ## Keywords
 
 - `flat = true`: whether to render sites and hops as flat circles and lines, or as 3D meshes
-- `force_transparency = false`: whether to disable occlusion of all non-transparent elements
+- `force_transparency = false`: whether to disable occlusion of all non-transparent elements. Useful in case `inspector` tooltips are occluded.
 - `shellopacity = 0.07`: opacity of elements surrounding the unitcell or the set of selected sites (dubbed "shell")
 - `cellopacity = 0.03`: opacity of the unitcell's boundingbox
 - `cellcolor = RGBAf(0,0,1)`: color of the unitcell's boundingbox
+- `boundarycolor = RGBAf(1,0,0)`: color of boundary cells for GreenFunction plots
+- `boundaryopacity = 0.07`: opacity of boundary cells for GreenFunction plots
 - `sitecolor = missing`: color of sites, as a index in `sitecolormap`, a named color, a `Makie.Colorant`, or as a site shader (see below). If `missing`, cycle through `sitecolormap`.
 - `sitecolormap = :Spectral_9`: colormap to use for `sitecolor` (see options in https://tinyurl.com/cschemes)
 - `siteopacity = missing`: opacity of sites, as a real between 0 and 1, or as a site shader (see below). If `missing`, obey `shellopacity`.
@@ -79,7 +86,8 @@ inter-site links.
 - `minmaxhopradius = (0, 0.1)`: if `hopdradius` is a shader, minimum and maximum hop radius.
 - `hopdarken = 0.85`: darkening factor for hops
 - `selector = missing`: an optional `siteselector(; sites...)` to filter which sites are shown (see `siteselector`)
-- `hide = (:cell,)`: collection of elements to hide, to choose from `(:hops, :sites, :hops, :bravais, :cell, :axes, :shell, :all)`. It can also be an empty collection or `nothing` to show all elements.
+- `hide = (:cell,)`: collection of elements to hide, to choose from `(:hops, :sites, :hops, :bravais, :cell, :axes, :shell, :boundary, :contacts, :all)`. It can also be an empty collection or `nothing` to show all elements.
+- `children = missing`: collection of `NamedTuple`s of any of the above keywords to be applied (cyclically) to contacts in GreenFunction plots
 
 ## Shaders
 

ext/QuanticaMakieExt/plotlattice.jl

@@ -9,10 +9,13 @@
         shellopacity = 0.07,
         cellopacity = 0.03,
         cellcolor = RGBAf(0,0,1),
+        boundarycolor = RGBAf(1,0,0),
+        boundaryopacity = 0.07,
         sitecolor = missing,         # accepts (i, r) -> float, IndexableObservable
         siteopacity = missing,       # accepts (i, r) -> float, IndexableObservable
         minmaxsiteradius = (0.0, 0.5),
         siteradius = 0.2,            # accepts (i, r) -> float, IndexableObservable
+        siteradiusfactor = 1.0,      # independent multiplier to apply to siteradius
         siteoutline = 2,
         siteoutlinedarken = 0.6,
         sitedarken = 0.0,
@@ -25,7 +28,7 @@
         hopcolormap = :Spectral_9,
         hoppixels = 2,
         selector = missing,
-        hide = :cell,               # :hops, :sites, :bravais, :cell, :axes, :shell, :all
+        hide = :cell,               # :hops, :sites, :bravais, :cell, :axes, :shell, :boundary, :contacts, :all
         isAxis3 = false,            # for internal use, to fix marker scaling
         marker = :auto,
         children = missing,
@@ -50,32 +53,6 @@ end
 
 Quantica.qplot!(x::PlotLatticeArgumentType; kw...) = plotlattice!(x; kw...)
 
-function green_selector(g)
-    mincell, maxcell = green_bounding_box(g)
-    s = siteselector(cells = n -> all(mincell .<= n .<= maxcell))
-    return s
-end
-
-not_boundary_selector(g) = siteselector(cells = n -> !isboundarycell(n, g))
-
-green_bounding_box(g) =
-    broaden_bounding_box(Quantica.boundingbox(g), Quantica.boundaries(g)...)
-
-function broaden_bounding_box((mincell, maxcell)::Tuple{SVector{N},SVector{N}}, (dir, cell), bs...) where {N}
-    isfinite(cell) || return (mincell, maxcell)
-    mincell´ = SVector(ntuple(i -> i == dir ? min(mincell[i], cell + 1) : mincell[i], Val(N)))
-    maxcell´ = SVector(ntuple(i -> i == dir ? max(maxcell[i], cell - 1) : maxcell[i], Val(N)))
-    return mincell´, maxcell´
-end
-
-broaden_bounding_box(bb) = bb
-
-isboundarycell(n, g) = any(((dir, cell),) -> n[dir] == cell, Quantica.boundaries(g))
-
-get_plottables_and_kws(Σp::Tuple) = (Σp, NamedTuple())
-get_plottables_and_kws(Σp::Quantica.FrankenTuple) = (Tuple(Σp), NamedTuple(Σp))
-get_plottables_and_kws(Σp) = ((Σp,), NamedTuple())
-
 #endregion
 
 ############################################################################################
@@ -510,28 +487,46 @@ end
 
 function Makie.plot!(plot::PlotLattice{Tuple{G}}) where {G<:GreenFunction}
     g = to_value(plot[1])
-    Σkws = Iterators.cycle(parse_children(plot[:children]))
-    Σs = Quantica.selfenergies(Quantica.contacts(g))
-    bsel = not_boundary_selector(g)
-    # plot lattice
     gsel = haskey(plot, :selector) && plot[:selector][] !== missing ?
         plot[:selector][] : green_selector(g)
     h = hamiltonian(g)
     latslice = lattice(h)[gsel]
     latslice´ = Quantica.growdiff(latslice, h)
-    latslice´ = latslice´[bsel]
-    latslice = latslice[bsel]
-    hideh = Quantica.tupleflatten(:cell, :bravais, plot[:hide][])
-    plotkw´´ = (; hopopacity = 0.2, siteopacity = 0.2, shellopacity = 0.07, plot.attributes...,
-        hide = hideh)
-    plotlattice!(plot, h, latslice, latslice´; plotkw´´...)
+    L = Quantica.latdim(h)
+    squaremarker = !plot[:flat][] ? Rect3f(Vec3f(-0.5), Vec3f(1)) : Rect2
+
+    # plot boundary as cells
+    if !ishidden((:boundary, :boundaries), plot)
+        bsel = boundary_selector(g)
+        blatslice = latslice[bsel]
+        blatslice´ = latslice´[bsel]
+        if L > 1
+            bkws = (; plot.attributes..., hide = (:axes, :sites, :hops), cellcolor = :boundarycolor, cellopacity = :boundaryopacity)
+            isempty(blatslice) || plotlattice!(plot, blatslice; bkws...)
+            isempty(blatslice´) || plotlattice!(plot, blatslice´; bkws...)
+        end
+        hideB = Quantica.tupleflatten(:bravais, plot[:hide][])
+        bkws´ = (; plot.attributes..., hide = hideB, sitecolor = :boundarycolor, siteopacity = :boundaryopacity,
+            siteradiusfactor = sqrt(2), marker = squaremarker)
+        isempty(blatslice) || plotlattice!(plot, blatslice; bkws´...)
+        isempty(blatslice´) || plotlattice!(plot, blatslice´; bkws´...)
+
+    end
+
+    # plot cells
+    plotlattice!(plot, h, latslice, latslice´; plot.attributes...)
+
     # plot contacts
-    for (Σ, Σkw) in zip(Σs, Σkws)
-        Σplottables = Quantica.selfenergy_plottables(Σ, bsel)
-        marker = !plot[:flat][] ? Rect3f(Vec3f(-0.5), Vec3f(1)) : Rect2
-        for Σp in Σplottables
-            plottables, kws = get_plottables_and_kws(Σp)
-            plotlattice!(plot, plottables...; plot.attributes..., marker, kws..., Σkw...)
+    if !ishidden((:contact, :contacts), plot)
+        Σkws = Iterators.cycle(parse_children(plot[:children]))
+        Σs = Quantica.selfenergies(Quantica.contacts(g))
+        hideΣ = Quantica.tupleflatten(:bravais, plot[:hide][])
+        for (Σ, Σkw) in zip(Σs, Σkws)
+            Σplottables = Quantica.selfenergy_plottables(Σ)
+            for Σp in Σplottables
+                plottables, kws = get_plottables_and_kws(Σp)
+                plotlattice!(plot, plottables...; plot.attributes..., hide = hideΣ, marker = squaremarker, kws..., Σkw...)
+            end
         end
     end
     return plot
@@ -557,9 +552,54 @@ function joint_colors_radii_update!(p, p´, plot)
     return nothing
 end
 
+############################################################################################
+# green_selector and boundary_selector: sites plotted for green functions
+#region
+
+function green_selector(g)
+    mincell, maxcell = green_bounding_box(g)
+    s = siteselector(cells = n -> all(mincell .<= n .<= maxcell))
+    return s
+end
+
+boundary_selector(g) = siteselector(cells = n -> isboundarycell(n, g))
+
+function green_bounding_box(g)
+    L = Quantica.latdim(hamiltonian(g))
+    return isempty(Quantica.contacts(g)) ?
+        boundary_bounding_box(Val(L), Quantica.boundaries(g)...) :
+        broaden_bounding_box(Quantica.boundingbox(g), Quantica.boundaries(g)...)
+end
+
+function broaden_bounding_box((mincell, maxcell)::Tuple{SVector{N},SVector{N}}, (dir, cell), bs...) where {N}
+    isfinite(cell) || return (mincell, maxcell)
+    mincell´ = SVector(ntuple(i -> i == dir ? min(mincell[i], cell + 1) : mincell[i], Val(N)))
+    maxcell´ = SVector(ntuple(i -> i == dir ? max(maxcell[i], cell - 1) : maxcell[i], Val(N)))
+    return broaden_bounding_box((mincell´, maxcell´), bs...)
+end
+
+broaden_bounding_box(mm::Tuple) = mm
+
+function boundary_bounding_box(::Val{L}, (dir, cell), bs...) where {L}
+    m = SVector(ntuple(i -> i == dir ? cell + 1 : 0, Val(L)))
+    return broaden_bounding_box((m, m), bs...)
+end
+
+isboundarycell(n, g) = any(((dir, cell),) -> n[dir] == cell, Quantica.boundaries(g))
+
+get_plottables_and_kws(Σp::Tuple) = (Σp, NamedTuple())
+get_plottables_and_kws(Σp::Quantica.FrankenTuple) = (Tuple(Σp), NamedTuple(Σp))
+get_plottables_and_kws(Σp) = ((Σp,), NamedTuple())
+
+#endregion
+
+############################################################################################
+# core plotting methods
+#region
+
 function plotsites_shaded!(plot::PlotLattice, sp::SitePrimitives, transparency)
     inspector_label = (self, i, r) -> sp.tooltips[i]
-    sizefactor = 1.0
+    sizefactor = plot[:siteradiusfactor][]
     if plot[:marker][] == :auto  # circle markers
         marker = (;)
     else
@@ -574,10 +614,11 @@ function plotsites_shaded!(plot::PlotLattice, sp::SitePrimitives, transparency)
     return plot
 end
 
-function plotsites_flat!(plot::PlotLattice, sp::SitePrimitives, transparency)
+function plotsites_flat!(plot::PlotLattice, sp::SitePrimitives{E}, transparency) where {E}
     inspector_label = (self, i, r) -> sp.tooltips[i]
-    sizefactor = ifelse(plot[:isAxis3][] && plot[:flat][], 0.5, sqrt(2))
-    if plot[:marker][] == :auto  # circle markers
+    sizefactor = ifelse(plot[:isAxis3][] && plot[:flat][], 0.5, sqrt(2)) * plot[:siteradiusfactor][]
+
+    if plot[:marker][] == :auto  # default circular markers
         marker = (;)
     else
         marker = (; marker = plot[:marker][])
@@ -632,9 +673,8 @@ function plotbravais!(plot::PlotLattice, lat::Lattice{<:Any,E,L}, latslice) wher
             arrows!(plot, [r0], [v]; color, inspectable = false)
         end
     end
-
     if !ishidden((:cell, :cells), plot) && L > 1
-        cellcolor = plot[:cellcolor][]
+        cellcolor = parse(RGBAf, plot[:cellcolor][])
         colface = transparent(cellcolor, plot[:cellopacity][])
         coledge = transparent(cellcolor, 5 * plot[:cellopacity][])
         rect = Rect{L,Int}(Point{L,Int}(0), Point{L,Int}(1))
@@ -646,7 +686,7 @@ function plotbravais!(plot::PlotLattice, lat::Lattice{<:Any,E,L}, latslice) wher
             mrect = GeometryBasics.mesh(rect, pointtype=Point{E,Float32}, facetype=QuadFace{Int})
             vertices = mrect.position
             vertices .= Ref(r0) .+ Ref(mat) .* (vertices0 .+ Ref(cell))
-            mesh!(plot, mrect; color = colface, transparency = true, inspectable = false)
+            mesh!(plot, mrect; color = colface, transparency = true, shading = NoShading, inspectable = false)
             wireframe!(plot, mrect; color = coledge, transparency = true, strokewidth = 1, inspectable = false)
         end
     end

ext/QuanticaMakieExt/tools.jl

@@ -11,7 +11,7 @@ end
 
 darken(colors::Vector, v = 0.66) = darken.(colors, Ref(v))
 
-transparent(rgba::RGBAf, v = 0.5) = RGBAf(rgba.r, rgba.g, rgba.b, rgba.alpha * v)
+transparent(rgba::RGBAf, v = 0.5) = RGBAf(rgba.r, rgba.g, rgba.b, clamp(rgba.alpha * v, 0f0, 1f0))
 
 maybedim(color, dn, dimming) = iszero(dn) ? color : transparent(color, 1 - dimming)
 
@@ -37,7 +37,8 @@ has_transparencies(::Missing) = false
 has_transparencies(x) = true
 
 function resolve_cross_references!(plot::PlotLattice)
-    names = (:shellopacity, :siteopacity, :hopopacity, :cellopacity, :sitecolor, :hopcolor, :siteradius, :hopradius)
+    names = (:shellopacity, :siteopacity, :hopopacity, :cellcolor, :cellopacity,
+        :boundarycolor, :boundaryopacity, :sitecolor, :hopcolor, :siteradius, :hopradius)
     for name in names
         property = plot[name][]
         if property isa Symbol && property in names

src/show.jl

@@ -96,20 +96,23 @@ function print_selector(io::IO, s::SiteSelector)
     i = get(io, :indent, "")
     print(io,
 "$(i)  Region            : $(s.region === missing ? "any" : "Function")
-$(i)  Sublattices       : $(s.sublats === missing ? "any" : s.sublats)
-$(i)  Cells             : $(s.cells === missing ? "any" : s.cells)")
+$(i)  Sublattices       : $(s.sublats === missing ? "any" : display_maybe_function(s.sublats))
+$(i)  Cells             : $(s.cells === missing ? "any" : display_maybe_function(s.cells))")
 end
 
 function print_selector(io::IO, s::HopSelector)
     i = get(io, :indent, "")
     print(io,
 "$(i)  Region            : $(s.region === missing ? "any" : "Function")
-$(i)  Sublattice pairs  : $(s.sublats === missing ? "any" : s.sublats)
-$(i)  Cell distances    : $(s.dcells === missing ? "any" : s.dcells)
+$(i)  Sublattice pairs  : $(s.sublats === missing ? "any" : display_maybe_function(s.sublats))
+$(i)  Cell distances    : $(s.dcells === missing ? "any" : display_maybe_function(s.dcells))
 $(i)  Hopping range     : $(displayrange(s.range))
 $(i)  Reverse hops      : $(s.adjoint)")
 end
 
+display_maybe_function(::Function) = "Function"
+display_maybe_function(x) = x
+
 #endregion
 
 ############################################################################################

src/solvers/selfenergy.jl

@@ -11,7 +11,7 @@
 #       applied (e.g. SelfEnergyModelSolver). Otherwise one must assume that any parent
 #       ParametricHamiltonian to GreenFunction has already been call!-ed before calling s.
 #   Optional: AbstractSelfEnergySolver's can also implement `selfenergy_plottables`
-#     - selfenergy_plottables(s::AbstractSelfEnergySolver, parent_latslice, boundary_sel)
+#     - selfenergy_plottables(s::AbstractSelfEnergySolver, parent_latslice)
 #       -> collection of tuples to be passed to plotlattice!(axis, tup...) for visualization
 ############################################################################################
 
@@ -22,12 +22,10 @@
 #   SelfEnergy wraps the corresponding SelfEnergySolver, be it Regular or Extended
 ############################################################################################
 
-selfenergy_plottables(s::SelfEnergy, boundaryselector) =
-    selfenergy_plottables(s.solver, orbslice(s), boundaryselector)
+selfenergy_plottables(s::SelfEnergy) = selfenergy_plottables(s.solver, orbslice(s))
 
 # fallback
-selfenergy_plottables(s::AbstractSelfEnergySolver, parent_latslice, boundaryselector) =
-    (parent_latslice[boundaryselector],)
+selfenergy_plottables(s::AbstractSelfEnergySolver, parent_latslice) = (parent_latslice,)
 
 include("selfenergy/nothing.jl")
 include("selfenergy/model.jl")

src/solvers/selfenergy/generic.jl

@@ -71,9 +71,9 @@ minimal_callsafe_copy(s::SelfEnergyGenericSolver) =
         minimal_callsafe_copy(s.V),
         copy(s.V´g), copy(s.Σ))
 
-function selfenergy_plottables(s::SelfEnergyGenericSolver, ls::LatticeSlice, bsel::SiteSelector)
-    p1 = ftuple(s.hcoupling)
-    p2 = ls[bsel]
+function selfenergy_plottables(s::SelfEnergyGenericSolver, ls::LatticeSlice)
+    p1 = s.hcoupling
+    p2 = ls
     return (p1, p2)
 end
 

src/solvers/selfenergy/schur.jl

@@ -111,9 +111,9 @@ minimal_callsafe_copy(s::SelfEnergySchurSolver) =
     SelfEnergySchurSolver(minimal_callsafe_copy(s.fsolver), minimal_callsafe_copy(s.hlead),
         s.isleftside, s.boundary, s.leadtoparent)
 
-function selfenergy_plottables(s::SelfEnergySchurSolver, ls::LatticeSlice, bsel::SiteSelector)
-    p1 = s.hlead, lattice(s.hlead)[cells = SA[1]]
-    p2 = (ls[bsel], )
+function selfenergy_plottables(s::SelfEnergySchurSolver, ls::LatticeSlice)
+    p1 = ftuple(s.hlead; selector = siteselector(cells = SA[1]))
+    p2 = (ls, )
     return p1, p2
 end
 
@@ -233,10 +233,10 @@ minimal_callsafe_copy(s::SelfEnergyCouplingSchurSolver) =
         minimal_callsafe_copy(s.g⁻¹),
         minimal_callsafe_copy(s.V))
 
-function selfenergy_plottables(s::SelfEnergyCouplingSchurSolver, ls::LatticeSlice, bsel::SiteSelector)
+function selfenergy_plottables(s::SelfEnergyCouplingSchurSolver, ls::LatticeSlice)
     p1 = ftuple(s.hcoupling; hide = :sites)
     p2 = hamiltonian(s.gunit)
-    p3 = ls[bsel]
+    p3 = ls
     return (p1, p2, p3)
 end
 

test/test_plots.jl

@@ -38,6 +38,13 @@ end
     glead = LP.honeycomb() |> hopping(1, range = 1) |> supercell((1,-1), region = r -> 0<=r[2]<=5) |> attach(nothing, cells = SA[5]) |> greenfunction(GS.Schur(boundary = 0));
     g = LP.honeycomb() |> hopping(1) |> supercell(region = r -> -6<=r[1]<=6 && 0<=r[2]<=5) |> attach(glead, region = r -> r[1] > 4.9) |> greenfunction;
     @test qplot(g, shellopacity = 0.3) isa Figure
+    hlead = LP.square() |> supercell((1,0), region = r -> 0 <= r[2] < 2) |> hopping(1)
+    glead = LP.square() |> onsite(4) - hopping(1) |> supercell((1, 0), region = r -> abs(r[2]) <= 5/2) |> attach(nothing, cells = SA[2]) |> greenfunction(GS.Schur(boundary = -2))
+    @test qplot(glead, siteradius = 0.25, children = (; sitecolor = :blue)) isa Figure
+    g = LP.honeycomb() |> hopping(1, range = 1) |>
+        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
 end
 
 @testset "plot bands" begin