diff --git a/docs/src/assets/josephson_lat.png b/docs/src/assets/josephson_lat.png index 54e1047b..75368b06 100644 Binary files a/docs/src/assets/josephson_lat.png and b/docs/src/assets/josephson_lat.png differ diff --git a/docs/src/assets/multiterminal.png b/docs/src/assets/multiterminal.png index 24a03c4a..0da13ed6 100644 Binary files a/docs/src/assets/multiterminal.png and b/docs/src/assets/multiterminal.png differ diff --git a/docs/src/tutorial/greenfunctions.md b/docs/src/tutorial/greenfunctions.md index 7fbc2728..bcdf0e11 100644 --- a/docs/src/tutorial/greenfunctions.md +++ b/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. diff --git a/docs/src/tutorial/observables.md b/docs/src/tutorial/observables.md index f7f9fe69..6f151acd 100644 --- a/docs/src/tutorial/observables.md +++ b/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 Josephson junction diff --git a/ext/QuanticaMakieExt/docstrings.jl b/ext/QuanticaMakieExt/docstrings.jl index f882ded2..b007ff4c 100644 --- a/ext/QuanticaMakieExt/docstrings.jl +++ b/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 diff --git a/ext/QuanticaMakieExt/plotlattice.jl b/ext/QuanticaMakieExt/plotlattice.jl index 202ecc8c..b29dc91f 100644 --- a/ext/QuanticaMakieExt/plotlattice.jl +++ b/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 diff --git a/ext/QuanticaMakieExt/tools.jl b/ext/QuanticaMakieExt/tools.jl index e1e91b1c..85e02b3b 100644 --- a/ext/QuanticaMakieExt/tools.jl +++ b/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 diff --git a/src/apply.jl b/src/apply.jl index 025fcf75..913e9ce3 100644 --- a/src/apply.jl +++ b/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) diff --git a/src/selectors.jl b/src/selectors.jl index 54978141..b98f925f 100644 --- a/src/selectors.jl +++ b/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 diff --git a/src/show.jl b/src/show.jl index 647cb86d..a0c30b55 100644 --- a/src/show.jl +++ b/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 ############################################################################################ diff --git a/src/slices.jl b/src/slices.jl index 86b529b3..92eb7c2f 100644 --- a/src/slices.jl +++ b/src/slices.jl @@ -14,18 +14,20 @@ Base.getindex(lat::Lattice, ::SiteSelectorAll) = lat[siteselector(; cells = zero function Base.getindex(lat::Lattice, as::AppliedSiteSelector) L = latdim(lat) csites = CellSites{L,Vector{Int}}[] - sinds = Int[] - foreach_cell(as) do cell - isempty(sinds) || (sinds = Int[]) - cs = CellSites(cell, sinds) - foreach_site(as, cell) do s, i, r - push!(siteindices(cs), i) - end - if isempty(cs) - return false - else - push!(csites, cs) - return true + if !isnull(as) + sinds = Int[] + foreach_cell(as) do cell + isempty(sinds) || (sinds = Int[]) + cs = CellSites(cell, sinds) + foreach_site(as, cell) do s, i, r + push!(siteindices(cs), i) + end + if isempty(cs) + return false + else + push!(csites, cs) + return true + end end end cellsdict = CellSitesDict{L}(cell.(csites), csites) @@ -60,6 +62,7 @@ function Base.getindex(latslice::SiteSlice, as::AppliedSiteSelector) lat = parent(latslice) L = latdim(lat) cellsdict´ = CellSitesDict{L}() + isnull(as) && return SiteSlice(lat, cellsdict´) sinds = Int[] for subcell in cellsdict(latslice) dn = cell(subcell) @@ -82,6 +85,7 @@ function Base.getindex(latslice::OrbitalSliceGrouped, as::AppliedSiteSelector) lat = parent(latslice) L = latdim(lat) cellsdict´ = CellOrbitalsGroupedDict{L}() + isnull(as) && return OrbitalSliceGrouped(lat, cellsdict´) oinds = Int[] ogroups = Dictionary{Int,UnitRange{Int}}() for subcell in cellsdict(latslice) diff --git a/src/solvers/selfenergy.jl b/src/solvers/selfenergy.jl index 56a1e9c4..8c56cbfe 100644 --- a/src/solvers/selfenergy.jl +++ b/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") diff --git a/src/solvers/selfenergy/generic.jl b/src/solvers/selfenergy/generic.jl index 718e6bfa..343fe56c 100644 --- a/src/solvers/selfenergy/generic.jl +++ b/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 diff --git a/src/solvers/selfenergy/schur.jl b/src/solvers/selfenergy/schur.jl index 0eb202a5..508bd7c8 100644 --- a/src/solvers/selfenergy/schur.jl +++ b/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 diff --git a/src/supercell.jl b/src/supercell.jl index 5551a342..7d9fe2ea 100644 --- a/src/supercell.jl +++ b/src/supercell.jl @@ -88,6 +88,7 @@ function supercell_sitelist!!(sitelist, masklist, smatperp, seedperp, lat, appli masklist´ = copy(masklist) empty!(masklist) empty!(sitelist) + isnull(applied_selector) && return nothing cs = cells(applied_selector) csbool = zeros(Bool, length(cs)) iter = BoxIterator(seedperp) diff --git a/src/types.jl b/src/types.jl index 5348eef6..97a3eddd 100644 --- a/src/types.jl +++ b/src/types.jl @@ -194,6 +194,7 @@ struct AppliedSiteSelector{T,E,L} region::FunctionWrapper{Bool,Tuple{SVector{E,T}}} sublats::Vector{Int} cells::Vector{SVector{L,Int}} + isnull::Bool # if isnull, the selector selects nothing, regardless of other fields end struct HopSelector{F,S,D,R} @@ -210,6 +211,7 @@ struct AppliedHopSelector{T,E,L} sublats::Vector{Pair{Int,Int}} dcells::Vector{SVector{L,Int}} range::Tuple{T,T} + isnull::Bool # if isnull, the selector selects nothing, regardless of other fields end struct Neighbors @@ -252,6 +254,9 @@ iswithinrange(dr, (rmin, rmax)::Tuple{Real,Real}) = ifelse(sign(rmin)*rmin^2 <= isbelowrange(dr, s::AppliedHopSelector) = isbelowrange(dr, s.range) isbelowrange(dr, (rmin, rmax)::Tuple{Real,Real}) = ifelse(dr'dr < rmin^2, true, false) +isnull(s::AppliedSiteSelector) = s.isnull +isnull(s::AppliedHopSelector) = s.isnull + Base.adjoint(s::HopSelector) = HopSelector(s.region, s.sublats, s.dcells, s.range, !s.adjoint) Base.NamedTuple(s::SiteSelector) = @@ -1295,12 +1300,14 @@ Base.size(h::Harmonic, i...) = size(matrix(h), i...) Base.isless(h::Harmonic, h´::Harmonic) = sum(abs2, dcell(h)) < sum(abs2, dcell(h´)) -Base.zero(h::Harmonic{<:Any,<:Any,B}) where B = Harmonic(zero(dcell(h)), zero(matrix(h))) +Base.zero(h::Harmonic{<:Any,<:Any,B}) where {B} = Harmonic(zero(dcell(h)), zero(matrix(h))) Base.copy(h::Harmonic) = Harmonic(dcell(h), copy(matrix(h))) Base.:(==)(h::Harmonic, h´::Harmonic) = h.dn == h´.dn && unflat(h.h) == unflat(h´.h) +Base.iszero(h::Harmonic) = iszero(flat(h)) + #endregion #endregion @@ -1427,6 +1434,7 @@ end Base.size(h::Hamiltonian, i...) = size(bloch(h), i...) Base.axes(h::Hamiltonian, i...) = axes(bloch(h), i...) +Base.iszero(h::Hamiltonian) = all(iszero, harmonics(h)) Base.copy(h::Hamiltonian) = Hamiltonian( copy(lattice(h)), copy(blockstructure(h)), copy.(harmonics(h)), copy(bloch(h))) diff --git a/test/test_hamiltonian.jl b/test/test_hamiltonian.jl index 9d441ec0..6e448f4a 100644 --- a/test/test_hamiltonian.jl +++ b/test/test_hamiltonian.jl @@ -341,7 +341,11 @@ end # non-spatial models h = LP.linear() |> @hopping((i,j) --> ind(i) + ind(j)) + @onsite((i; k = 1) --> pos(i)[k]) @test ishermitian(h()) - + # null selectors + h0 = LP.square() |> onsite(0) + hopping(0) |> supercell(3) |> @onsite!((t, r) -> 1; sublats = Symbol[]) + @test iszero(h0()) + h0 = LP.square() |> hopping(0) |> supercell(3) |> @hopping!((t, r, dr) -> 1; dcells = SVector{2,Int}[]) + @test iszero(h0()) end diff --git a/test/test_lattice.jl b/test/test_lattice.jl index c2e98872..201a5652 100644 --- a/test/test_lattice.jl +++ b/test/test_lattice.jl @@ -183,4 +183,12 @@ end ls = lat[cellsites(SA[1,0], 2)] @test ls isa Quantica.SiteSlice @test nsites(ls) == 1 + # test the difference between a null selector and an unspecified one + ls = lat[cells = SVector{2,Int}[]] + @test isempty(ls) + ls = lat[cells = missing] + @test !isempty(ls) + ls = lat[region = RP.circle(4)] + ls´ = ls[cells = n -> !isreal(n[1])] + @test isempty(ls´) end diff --git a/test/test_plots.jl b/test/test_plots.jl index e386ff57..84d99c06 100644 --- a/test/test_plots.jl +++ b/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