Alternative representation of GreenFunction boundaries in plots

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/apply.jl

@@ -18,7 +18,11 @@ function apply(s::SiteSelector, lat::Lattice{T,E,L}, cells...) where {T,E,L}
     cells = recursive_push!(SVector{L,Int}[], sanitize_cells(s.cells, Val(L)), cells...)
     unique!(sort!(sublats))
     unique!(sort!(cells))
-    return AppliedSiteSelector{T,E,L}(lat, region, sublats, cells)
+    # isnull: to distinguish in a type-stable way between s.cells === missing and no-selected-cells
+    # and the same for sublats
+    isnull = (s.cells !== missing && isempty(cells)) ||
+        (s.sublats !== missing && isempty(sublats))
+    return AppliedSiteSelector{T,E,L}(lat, region, sublats, cells, isnull)
 end
 
 function apply(s::HopSelector, lat::Lattice{T,E,L}, cells...) where {T,E,L}
@@ -36,7 +40,9 @@ function apply(s::HopSelector, lat::Lattice{T,E,L}, cells...) where {T,E,L}
         sublats .= reverse.(sublats)
         dcells .*= -1
     end
-    return AppliedHopSelector{T,E,L}(lat, region, sublats, dcells, (rmin, rmax))
+    isnull = (s.dcells !== missing && isempty(dcells)) ||
+        (s.sublats !== missing && isempty(sublats))
+    return AppliedHopSelector{T,E,L}(lat, region, sublats, dcells, (rmin, rmax), isnull)
 end
 
 sublatindex_or_zero(lat, ::Missing) = missing
@@ -66,7 +72,8 @@ applyrange(r::Real, lat) = r
 padrange(r::Real, m) = isfinite(r) ? float(r) + m * sqrt(eps(float(r))) : float(r)
 
 applied_region(r, ::Missing) = true
-applied_region((r, dr)::Tuple{SVector,SVector}, region::Function) = ifelse(region(r, dr), true, false)
+applied_region((r, dr)::Tuple{SVector,SVector}, region::Function) =
+    ifelse(region(r, dr), true, false)
 applied_region(r::SVector, region::Function) = ifelse(region(r), true, false)
 
 recursive_apply(f, t::Tuple) = recursive_apply.(f, t)
@@ -181,6 +188,7 @@ end
 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.")
     ptrs = Tuple{Int,SVector{E,T},CellSitePos{T,E,L,B},Int}[]
+    isnull(s) && return ptrs
     lat = lattice(h)
     har0 = first(harmonics(h))
     dn0 = zerocell(lat)
@@ -204,6 +212,7 @@ end
 function pointers(h::Hamiltonian{T,E,L,B}, s::AppliedHopSelector{T,E,L}, shifts) where {T,E,L,B}
     hars = harmonics(h)
     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]
+    isnull(s) && return ptrs
     lat = lattice(h)
     dn0 = zerocell(lat)
     norbs = norbitals(h)

src/selectors.jl

@@ -24,12 +24,12 @@ neighbors(n::Int, lat::Lattice) = nrange(n, lat)
 #region
 
 function Base.in((s, r)::Tuple{Int,SVector{E,T}}, sel::AppliedSiteSelector{T,E}) where {T,E}
-    return inregion(r, sel) &&
+    return !isnull(sel) && inregion(r, sel) &&
            insublats(s, sel)
 end
 
 function Base.in((s, r, cell)::Tuple{Int,SVector{E,T},SVector{L,Int}}, sel::AppliedSiteSelector{T,E,L}) where {T,E,L}
-    return incells(cell, sel) &&
+    return !isnull(sel) && incells(cell, sel) &&
            inregion(r, sel) &&
            insublats(s, sel)
 end
@@ -51,7 +51,7 @@ end
 # end
 
 function Base.in(((sj, si), (r, dr), dcell)::Tuple{Pair,Tuple,SVector}, sel::AppliedHopSelector)
-    return !isonsite(dr) &&
+    return !isnull(sel) && !isonsite(dr) &&
             indcells(dcell, sel) &&
             insublats(sj => si, sel) &&
             iswithinrange(dr, sel) &&
@@ -70,6 +70,7 @@ isonsite(dr) = iszero(dr)
 # foreach_cell(f,...) should be called with a boolean function f that returns whether the
 # cell should be mark as accepted when BoxIterated
 function foreach_cell(f, sel::AppliedSiteSelector)
+    isnull(sel) && return nothing
     lat = lattice(sel)
     cells_list = cells(sel)
     if isempty(cells_list)      # no cells specified
@@ -92,6 +93,7 @@ end
 
 
 function foreach_cell(f, sel::AppliedHopSelector)
+    isnull(sel) && return nothing
     lat = lattice(sel)
     dcells_list = dcells(sel)
     if isempty(dcells_list) # no dcells specified
@@ -109,6 +111,7 @@ function foreach_cell(f, sel::AppliedHopSelector)
 end
 
 function foreach_site(f, sel::AppliedSiteSelector, cell::SVector)
+    isnull(sel) && return nothing
     lat = lattice(sel)
     for s in sublats(lat)
         insublats(s, sel) || continue
@@ -122,6 +125,7 @@ function foreach_site(f, sel::AppliedSiteSelector, cell::SVector)
 end
 
 function foreach_hop(f, sel::AppliedHopSelector, kdtrees::Vector{<:KDTree}, ni::SVector = zerocell(lattice(sel)))
+    isnull(sel) && return nothing
     lat = lattice(sel)
     _, rmax = sel.range
     # source cell at origin