diff --git a/docs/src/advanced/meanfield.md b/docs/src/advanced/meanfield.md index 5642276f..4b7b9363 100644 --- a/docs/src/advanced/meanfield.md +++ b/docs/src/advanced/meanfield.md @@ -41,7 +41,7 @@ julia> h(SA[]) julia> g = greenfunction(h, GS.Spectrum()); julia> ρ0 = densitymatrix(g[])(0.5, 0) ## density matrix at chemical potential `µ=0.5` and temperature `kBT = 0` on all sites -4×4 OrbitalSliceMatrix{Matrix{ComplexF64}}: +4×4 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.138197+0.0im 0.223607+0.0im 0.223607+0.0im 0.138197+0.0im 0.223607+0.0im 0.361803+0.0im 0.361803+0.0im 0.223607+0.0im 0.223607+0.0im 0.361803+0.0im 0.361803+0.0im 0.223607+0.0im diff --git a/docs/src/tutorial/greenfunctions.md b/docs/src/tutorial/greenfunctions.md index 06870bc9..c623c9b4 100644 --- a/docs/src/tutorial/greenfunctions.md +++ b/docs/src/tutorial/greenfunctions.md @@ -213,7 +213,7 @@ julia> g(0.2) GreenSolution{Float64,2,0}: Green function at arbitrary positions, but at a fixed energy julia> g(0.2)[1, 3] -5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}: +5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: -2.56906+0.000123273im -4.28767+0.00020578im -4.88512+0.000234514im -4.28534+0.00020578im -2.5664+0.000123273im -4.28767+0.00020578im -7.15613+0.00034351im -8.15346+0.000391475im -7.15257+0.00034351im -4.2836+0.000205781im -4.88512+0.000234514im -8.15346+0.000391475im -9.29002+0.000446138im -8.14982+0.000391476im -4.88095+0.000234514im @@ -221,7 +221,7 @@ julia> g(0.2)[1, 3] -2.5664+0.000123273im -4.2836+0.000205781im -4.88095+0.000234514im -4.28211+0.000205781im -2.56469+0.000123273im julia> g(0.2)[siteselector(region = RP.circle(1, (0.5, 0))), 3] -2×5 OrbitalSliceMatrix{Matrix{ComplexF64}}: +2×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.0749214+3.15744e-8im 0.124325+5.27948e-8im 0.141366+6.01987e-8im 0.124325+5.27948e-8im 0.0749214+3.15744e-8im -0.374862+2.15287e-5im -0.625946+3.5938e-5im -0.712983+4.09561e-5im -0.624747+3.59379e-5im -0.37348+2.15285e-5im ``` diff --git a/docs/src/tutorial/observables.md b/docs/src/tutorial/observables.md index adc05882..9440e739 100644 --- a/docs/src/tutorial/observables.md +++ b/docs/src/tutorial/observables.md @@ -42,7 +42,7 @@ DensityMatrix: density matrix on specified sites using solver of type DensityMat julia> @time ρ(4) 6.150548 seconds (57.84 k allocations: 5.670 GiB, 1.12% gc time) -5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}: +5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.5+0.0im -7.34893e-10-3.94035e-15im 0.204478+1.9366e-14im -7.34889e-10-1.44892e-15im -5.70089e-10+5.48867e-15im -7.34893e-10+3.94035e-15im 0.5+0.0im 0.200693-2.6646e-14im -5.70089e-10-1.95251e-15im -7.34891e-10-2.13804e-15im 0.204478-1.9366e-14im 0.200693+2.6646e-14im 0.5+0.0im 0.200693+3.55692e-14im 0.204779-4.27255e-14im @@ -62,7 +62,7 @@ DensityMatrix: density matrix on specified sites with solver of type DensityMatr julia> @time ρ(4) 0.001659 seconds (9 allocations: 430.906 KiB) -5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}: +5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.5+0.0im -2.21437e-15+0.0im 0.204478+0.0im 2.67668e-15+0.0im 3.49438e-16+0.0im -2.21437e-15+0.0im 0.5+0.0im 0.200693+0.0im -1.40057e-15+0.0im -2.92995e-15+0.0im 0.204478+0.0im 0.200693+0.0im 0.5+0.0im 0.200693+0.0im 0.204779+0.0im @@ -81,7 +81,7 @@ DensityMatrix: density matrix on specified sites with solver of type DensityMatr julia> @time ρ(4) 0.006580 seconds (3 allocations: 1.156 KiB) -5×5 OrbitalSliceMatrix{Matrix{ComplexF64}}: +5×5 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.5+0.0im 2.15097e-17+0.0im 0.20456+0.0im 2.15097e-17+0.0im 3.9251e-17+0.0im 2.15097e-17+0.0im 0.5+0.0im 0.200631+0.0im 1.05873e-16+0.0im 1.70531e-18+0.0im 0.20456+0.0im 0.200631+0.0im 0.5+0.0im 0.200631+0.0im 0.20482+0.0im diff --git a/ext/QuanticaMakieExt/plotlattice.jl b/ext/QuanticaMakieExt/plotlattice.jl index 0b56870d..31e61083 100644 --- a/ext/QuanticaMakieExt/plotlattice.jl +++ b/ext/QuanticaMakieExt/plotlattice.jl @@ -150,7 +150,7 @@ function _hoppingprimitives(ls::LatticeSlice{<:Any,E}, ls´, h, opts, siteradii) return hp, hp´ end -maybe_evaluate_observable(o::Quantica.IndexableObservable, ls) = o[ls] +maybe_evaluate_observable(o::Quantica.IndexableObservable, ls) = parent(o[ls]) maybe_evaluate_observable(x, ls) = x maybe_getindex(v::AbstractVector{<:Number}, i) = v[i] @@ -189,9 +189,9 @@ push_siteopacity!(sp, siteopacity::Real, shellopacity, i, r, is_shell) = push!(s push_siteopacity!(sp, siteopacity::Function, shellopacity, i, r, is_shell) = push!(sp.opacities, is_shell ? siteopacity(i, r) : shellopacity) push_siteopacity!(sp, siteopacity, shellopacity, i, r, is_shell) = argerror("Unrecognized siteradius") -push_siteradius!(sp, siteradius::Real, i, r) = push!(sp.radii, siteradius) +push_siteradius!(sp, siteradius::Number, i, r) = push!(sp.radii, siteradius) push_siteradius!(sp, siteradius::Function, i, r) = push!(sp.radii, siteradius(i, r)) -push_siteradius!(sp, siteradius, i, r) = argerror("Unrecognized siteradius") +push_siteradius!(sp, siteradius, i, r) = argerror("Unrecognized siteradius of type $(typeof(siteradius))") push_sitetooltip!(sp, i, r, mat) = push!(sp.tooltips, matrixstring(i, mat)) push_sitetooltip!(sp, i, r) = push!(sp.tooltips, positionstring(i, r)) @@ -230,9 +230,9 @@ push_hopopacity!(hp, hopopacity::Real, shellopacity, ij, rdr, is_shell) = push!( push_hopopacity!(hp, hopopacity::Function, shellopacity, ij, rdr, is_shell) = push!(hp.opacities, hopopacity(ij, rdr)) push_hopopacity!(hp, hopopacity, shellopacity, ij, rdr, is_shell) = argerror("Unrecognized hopradius") -push_hopradius!(hp, hopradius::Real, ij, rdr) = push!(hp.radii, hopradius) +push_hopradius!(hp, hopradius::Number, ij, rdr) = push!(hp.radii, hopradius) push_hopradius!(hp, hopradius::Function, ij, rdr) = push!(hp.radii, hopradius(ij, rdr)) -push_hopradius!(hp, hopradius, ij, rdr) = argerror("Unrecognized hopradius") +push_hopradius!(hp, hopradius, ij, rdr) = argerror("Unrecognized hopradius of type $(typeof(hopradius))") push_hoptooltip!(hp, (i, j), mat) = push!(hp.tooltips, matrixstring(i, j, mat)) diff --git a/src/Quantica.jl b/src/Quantica.jl index 252accdf..ca67908b 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -25,14 +25,15 @@ using DelimitedFiles export sublat, bravais_matrix, lattice, sites, supercell, hamiltonian, hopping, onsite, @onsite, @hopping, @onsite!, @hopping!, pos, ind, cell, position, - plusadjoint, neighbors, siteselector, hopselector, diagonal, + plusadjoint, neighbors, siteselector, hopselector, diagonal, sitepairs, unflat, torus, transform, translate, combine, spectrum, energies, states, bands, subdiv, greenfunction, selfenergy, attach, plotlattice, plotlattice!, plotbands, plotbands!, qplot, qplot!, qplotdefaults, conductance, josephson, ldos, current, transmission, densitymatrix, - OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix, orbaxes, - serializer, serialize, serialize!, deserialize!, deserialize + OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix, orbaxes, siteindexdict, + serializer, serialize, serialize!, deserialize!, deserialize, + meanfield export LatticePresets, LP, RegionPresets, RP, HamiltonianPresets, HP, ExternalPresets, EP export EigenSolvers, ES, GreenSolvers, GS @@ -63,6 +64,7 @@ include("mesh.jl") include("bands.jl") include("greenfunction.jl") include("observables.jl") +include("meanfield.jl") # Plumbing include("apply.jl") include("show.jl") diff --git a/src/apply.jl b/src/apply.jl index cea62a4a..d81bcd3b 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -218,7 +218,7 @@ function push_pointers!(ptrs, h, har0, s::AppliedSiteSelector, shifts = missing, lat = lattice(h) umat = unflat(har0) rows = rowvals(umat) - norbs = norbitals(h) + norbs = norbitals(h, :) for scol in sublats(lat), col in siterange(lat, scol) isinblock(col, block) || continue for p in nzrange(umat, col) @@ -243,7 +243,7 @@ function push_pointers!(ptrs, h, har, s::AppliedHopSelector, shifts = missing, b lat = lattice(h) dn0 = zerocell(lat) dn = dcell(har) - norbs = norbitals(h) + norbs = norbitals(h, :) umat = unflat(har) rows = rowvals(umat) for scol in sublats(lat), col in siterange(lat, scol), p in nzrange(umat, col) diff --git a/src/docstrings.jl b/src/docstrings.jl index 715ba4fb..75f8d2e1 100644 --- a/src/docstrings.jl +++ b/src/docstrings.jl @@ -716,7 +716,7 @@ julia> h[sites(1), sites(2)] 1.0+0.0im 0.0+0.0im julia> ph = h |> @hopping!((t; p = 3) -> p*t); ph[region = RP.square(1)] -4×4 OrbitalSliceMatrix{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}: +4×4 OrbitalSliceMatrix{ComplexF64,SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}: 0.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 0.0+0.0im 0.0+0.0im 3.0+0.0im 0.0+0.0im 0.0+0.0im @@ -1632,6 +1632,21 @@ For any `gω::GreenSolution` and `C<:Union{Colon,Integer}`, obtain a view (of ty For any `g::Union{GreenFunction,GreenSlice}`, produce a new `GreenFunction` or `GreenSlice` with all parameters fixed to `params` (or to their default values if not provided). +## Full evaluation + + gs(ω; params...) + gω[sites...] + +For `gs::GreenSlice` or `gω::GreenSolution`, return a fully evaluated `m::AbstractMatrix`. +If the selected site slice was defined using `sites`, the concrete type of `m` will be will +be a conventional `Matrix`-based type. Otherwise, it will be of type `OrbitalSliceMatrix`, +an `AbstractMatrix` type that supports both conventional indexing and indexing with `sites` +and `siteselectors`. + +Advanced: in addition to the above, an unexported method `Quantica.call!(gs, ω; params...)` +is provided to reuse the output matrix `m` (preallocated inside `gs`). Use with caution, as +it may lead to unexpected aliasing + # Example ```jldoctest julia> g = LP.honeycomb() |> hamiltonian(@hopping((; t = 1) -> t)) |> supercell(region = RP.circle(10)) |> greenfunction(GS.SparseLU()) @@ -1660,10 +1675,13 @@ GreenSlice{Float64,2,0}: Green function at arbitrary energy, but at a fixed latt julia> gω[ss] == gs(0.1; t = 2) true + +julia> summary(gω[ss]) +"14×14 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}" ``` # See also - `GreenSolvers`, `diagonal`, `ldos`, `conductance`, `current`, `josephson` + `GreenSolvers`, `diagonal`, `sitepairs`, `ldos`, `conductance`, `current`, `josephson` """ greenfunction @@ -1696,47 +1714,91 @@ GreenSolvers """ diagonal(i; kernel = missing) -Wrapper over site or orbital indices (used to index into a `g::GreenFunction` or +Wrapper over site or orbital indices `i` (used to index into a `g::GreenFunction` or `g::GreenSolution`) that represent purely diagonal entries. Here `i` can be any index accepted in `g[i,i]`, e.g. `i::Integer` (contact index), `i::Colon` (merged contacts), `i::SiteSelector` (selected sites), etc. +If `kernel = Q` (a matrix) instead of `missing`, each diagonal block for multiorbital site +`i` is replaced with `Tr(gᵢᵢQ)`. + +For a `gω::GreenSolution`, `gω[diagonal(sel)] = diag(gω[sel, sel])`, although where possible +the former computation is done more efficiently internally. + diagonal(kernel = missing, sites...) Equivalent to `diagonal(siteselector(; sites...); kernel)` ## Keywords - - `kernel`: if missing, all orbitals in the diagonal `g[i, i]` are returned when indexing `g[diagonal(i)]`. Otherwise, `tr(g[site, site]*kernel)` for each site included in `i` is returned. See also `ldos` + - `kernel`: if missing, all orbitals in the diagonal `g[i, i]` are returned when indexing `g[diagonal(i)]`. Otherwise, `Tr(g[site, site]*kernel)` for each site included in `i` is returned. # Example -```jldoctest +```julia julia> g = HP.graphene(orbitals = 2) |> attach(nothing, cells = (0,0)) |> greenfunction(); -julia> g(1)[diagonal(:)] # g diagonal on all contact orbitals -4-element Vector{ComplexF64}: - -0.10919028168061964 - 0.08398577667508965im - -0.10919028168734393 - 0.08398577667508968im - -0.10919028169083328 - 0.08398577667508969im - -0.109190281684109 - 0.08398577667508969im +julia> g(1)[diagonal(:)] # g(ω = 1) diagonal on all contact orbitals +4×4 OrbitalSliceMatrix{ComplexF64,LinearAlgebra.Diagonal{ComplexF64, Vector{ComplexF64}}}: + -0.10919-0.0839858im 0.0+0.0im 0.0+0.0im 0.0+0.0im + 0.0+0.0im -0.10919-0.0839858im 0.0+0.0im 0.0+0.0im + 0.0+0.0im 0.0+0.0im -0.10919-0.0839858im 0.0+0.0im + 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.10919-0.0839858im + +julia> g(1)[diagonal(:, kernel = SA[1 0; 0 -1])] # σz spin density of the above +2×2 OrbitalSliceMatrix{ComplexF64,LinearAlgebra.Diagonal{ComplexF64, Vector{ComplexF64}}}: + 5.61885e-12+1.38778e-17im 0.0+0.0im + 0.0+0.0im -5.61882e-12+2.77556e-17im +``` -julia> g(1)[diagonal(:, kernel = SA[1 0; 0 -1])] # σz spin spectral density at ω = 1 -2-element Vector{ComplexF64}: - 6.724287793247186e-12 + 2.7755575615628914e-17im - -6.724273915459378e-12 + 0.0im +# See also + `sitepairs`, `greenfunction`, `ldos`, `densitymatrix` +""" +diagonal + +""" + + sitepairs(s::HopSelector; kernel = missing) + +Create a selection of site pairs `s::SparseIndices` used to sparsely index into a +`g::GreenFunction` or `g::GreenSolution`, as `g[s]`. Of the resulting `OrbitalSliceMatrix` +only the selected pairs of matrix elements will be computed, leaving the rest as zero +(sparse matrix). The sparse matrix spans the minimum number of complete unit cells to +include all site pairs + +If `kernel = Q` (a matrix instead of `missing`), each of these site blocks `gᵢⱼ` will be +replaced by `Tr(kernel * gᵢⱼ)`. + + sitepairs(; kernel = missing, hops...) + +Equivalent to `sitepairs(hopselector(; hops...); kernel)` + +## Keywords + + - `kernel`: if missing, all orbitals blocks `gᵢⱼ = g[i, j]` between selected sites pairs (i,j) are returned when indexing `g[sitepairs(...)]`. Otherwise, `gᵢⱼ` is replaced by `Tr(gᵢⱼ*kernel)`. + +# Example +```julia +julia> g = HP.graphene(orbitals = 2, a0 = 1) |> attach(nothing, cells = (0,0)) |> greenfunction(); + +julia> summary(g(1)[sitepairs(range = 1)]) # g(ω=1) site blocks between all sites in zero cell and all other sites at distance 1 +"28×4 OrbitalSliceMatrix{ComplexF64,SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}" + +julia> summary(g(1)[sitepairs(range = 1, kernel = SA[1 0; 0 -1])]) # σz spin density of the above +"14×2 OrbitalSliceMatrix{ComplexF64,SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}" ``` # See also - `greenfunction`, `ldos` + `diagonal`, `hopselector`, `greenfunction`, `ldos`, `densitymatrix` """ +sitepairs """ - ldos(gs::GreenSlice; kernel = I) + ldos(gs::GreenSlice; kernel = missing) Build `ρs::LocalSpectralDensitySlice`, a partially evaluated object representing the local density of states `ρᵢ(ω)` at specific sites `i` but at arbitrary energy `ω`. - ldos(gω::GreenSolution; kernel = I) + ldos(gω::GreenSolution; kernel = missing) Build `ρω::LocalSpectralDensitySolution`, as above, but for `ρᵢ(ω)` at a fixed `ω` and arbitrary sites `i`. See also `greenfunction` for details on building a `GreenSlice` and @@ -1747,7 +1809,7 @@ the retarded Green function at a given site `i`. ## Keywords -- `kernel` : for multiorbital sites, `kernel` allows to compute a generalized `ldos` `ρᵢ(ω) = -imag(Tr(gᵢᵢ(ω) * kernel))/π`, where `gᵢᵢ(ω)` is the retarded Green function at site `i` and energy `ω`. If `kernel = missing`, the complete, orbital-resolved `ldos` is returned. +- `kernel` : for multiorbital sites, `kernel` allows to compute a generalized `ldos` `ρᵢ(ω) = -imag(Tr(gᵢᵢ(ω) * kernel))/π`, where `gᵢᵢ(ω)` is the retarded Green function at site `i` and energy `ω`. If `kernel = missing`, the complete, orbital-resolved `ldos` is returned. Default: `missing` ## Full evaluation @@ -1784,7 +1846,7 @@ julia> ldos(g(0.2))[1] 0.03493305572265045 0.036802204179317045 -julia> ldos(g(0.2))[1] == -imag.(g[diagonal(1; kernel = I)](0.2)) ./ π +julia> ldos(g(0.2))[1] == -imag.(diag(g[diagonal(1)](0.2))) ./ π true ``` @@ -1996,7 +2058,7 @@ as an `OrbitalSliceMatrix`, see its docstring for further details. If the generic integration algorithm is used with complex `ωpoints`, the following form is also available: - ρ(μ = 0, kBT = 0, override = missing; params...) + ρ(μ = 0, kBT = 0, override_path = missing; params...) Here `override` can be a collection of `ωpoints` that will replace the original ones provided when defining `ρ`. Alternatively, it can be a function that will be applied to the @@ -2029,7 +2091,7 @@ julia> ρ = densitymatrix(g[region = RP.circle(0.5)]) DensityMatrix: density matrix on specified sites with solver of type DensityMatrixSpectrumSolver julia> ρ() # with mu = kBT = 0 by default -2×2 OrbitalSliceMatrix{Matrix{ComplexF64}}: +2×2 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.5+0.0im -0.262865+0.0im -0.262865+0.0im 0.5+0.0im ``` @@ -2057,15 +2119,16 @@ As above, but with `ωpoints = (-ωmax, ωmax)`. ## Full evaluation - J(kBT = 0, override = missing; params...) # where J::Josephson + J(kBT = 0, override_path = missing; params...) # where J::Josephson Evaluate the current `I_J` at chemical potemtial `µ = 0` and temperature `kBT` (in the same units as the Hamiltonian) for the given `g` parameters `params`, if any. -It's possible to `override` a complex integration path at evaluation time. In this case -`override` can be a collection of `ωpoints` that will replace the original ones provided -when defining `J`. Alternatively, it can be a function that will be applied to the original -`ωpoints`. This may be useful when the integration path must depend on `params`. +It's possible to use `override_path` to override a complex integration path at evaluation +time. In this case `override_path` can be a collection of `ωpoints` that will replace the +original ones provided when defining `J`. Alternatively, it can be a function that will be +applied to the original `ωpoints`. This may be useful when the integration path must depend +on `params`. ## Keywords @@ -2194,7 +2257,7 @@ Note: `diagonal` indexing is currently not supported by `OrbitalSliceArray`. julia> g = LP.linear() |> hamiltonian(hopping(SA[0 1; 1 0]) + onsite(I), orbitals = 2) |> supercell(4) |> greenfunction; julia> mat = g(0.2)[region = r -> 2<=r[1]<=4] -6×6 OrbitalSliceMatrix{Matrix{ComplexF64}}: +6×6 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: -1.93554e-9-0.545545im 0.0-0.0im 0.0-0.0im -0.5+0.218218im 0.4+0.37097im 0.0+0.0im 0.0-0.0im -1.93554e-9-0.545545im -0.5+0.218218im 0.0-0.0im 0.0+0.0im 0.4+0.37097im 0.0-0.0im -0.5+0.218218im -1.93554e-9-0.545545im 0.0-0.0im 0.0+0.0im -0.5+0.218218im @@ -2203,7 +2266,7 @@ julia> mat = g(0.2)[region = r -> 2<=r[1]<=4] 0.0+0.0im 0.4+0.37097im -0.5+0.218218im 0.0+0.0im 0.0-0.0im -1.93554e-9-0.545545im julia> mat[(; cells = SA[1]), (; cells = SA[0])] -2×4 OrbitalSliceMatrix{Matrix{ComplexF64}}: +2×4 OrbitalSliceMatrix{ComplexF64,Matrix{ComplexF64}}: 0.4+0.37097im 0.0+0.0im 0.0+0.0im -0.5+0.218218im 0.0+0.0im 0.4+0.37097im -0.5+0.218218im 0.0+0.0im @@ -2215,23 +2278,58 @@ julia> mat[sites(SA[1], 1)] # See also `siteselector`, `cellindices`, `orbaxes` - """ OrbitalSliceArray, OrbitalSliceVector, OrbitalSliceMatrix """ - diagonal(inds; kernel = missing) + orbaxes(A::OrbitalSliceArray) -Wrapper over indices `inds` used to obtain the diagonal of a `gω::GreenSolution`. If `d = -diagonal(sel)`, then `gω[d] = diag(gω[sel, sel])`, although in most cases the computation -is done more efficiently internally. +Return the orbital axes of `A`. This is a tuple of `OrbitalSliceGrouped` objects that can be +used e.g. to index another `OrbitalSliceArray` or to inspect the indices of each site with +`siteindexdict`. -If `kernel = Q` (a matrix) instead of `missing`, a vector with `Tr(gᵢᵢQ)` on -each site `i` is returned (instead of of `vcat(diag(gᵢᵢ)...)`). +# Examples +```jldoctest +julia> g = HP.graphene(orbitals = 2) |> supercell((1,-1)) |> greenfunction; + +julia> d = ldos(g[cells = SA[0]])(2); summary(d) +"8-element OrbitalSliceVector{Float64,Vector{Float64}}" + +julia> a = only(orbaxes(d)) +OrbitalSliceGrouped{Float64,2,1} : collection of subcells of orbitals (grouped by sites) for a 1D lattice in 2D space + Cells : 1 + Cell range : ([0], [0]) + Total sites : 4 + +julia> siteindexdict(a) +4-element Dictionaries.Dictionary{Quantica.CellIndices{1, Int64, Quantica.SiteLike}, UnitRange{Int64}} + CellSites{1,Int64} : 1 site in cell zero + Sites : 1 │ 1:2 + CellSites{1,Int64} : 1 site in cell zero + Sites : 2 │ 3:4 + CellSites{1,Int64} : 1 site in cell zero + Sites : 3 │ 5:6 + CellSites{1,Int64} : 1 site in cell zero + Sites : 4 │ 7:8 +``` + +# See also + `siteindexdict` """ -diagonal +orbaxes + +""" + siteindexdict(axis::OrbitalSliceGrouped) + +Return a dictionary of `CellSite`s representing single sites in an orbital axis of an +`OrbitalSliceArray`, typically obtained with `orbaxes`. See `orbaxes` for an example. + +# See also + `orbaxes`, `OrbitalSliceArray` +""" +siteindexdict """ serializer(T::Type, selectors...; parameter = :stream, encoder = identity, decoder = identity) diff --git a/src/greenfunction.jl b/src/greenfunction.jl index b484d250..5e9025be 100644 --- a/src/greenfunction.jl +++ b/src/greenfunction.jl @@ -60,11 +60,12 @@ function call!(g::GreenFunction{T}, ω::Complex{T}; params...) where {T} return GreenSolution(g, slicer, Σblocks, corbs) end -call!(g::GreenSlice{T}, ω::T; params...) where {T} = - call!(g, retarded_omega(ω, solver(parent(g))); params...) +# real frequency -> maybe complex frequency +call!(g::GreenSlice{T}, ω::T; kw...) where {T} = + call!(g, retarded_omega(ω, solver(parent(g))); kw...) -call!(g::GreenSlice{T}, ω::Complex{T}; params...) where {T} = - call!(greenfunction(g), ω; params...)[orbinds_or_contactinds(g)...] +call!(gs::GreenSlice{T}, ω::Complex{T}; post = identity, params...) where {T} = + getindex!(gs, call!(greenfunction(gs), ω; params...); post) real_or_complex_convert(::Type{T}, ω::Real) where {T<:Real} = convert(T, ω) real_or_complex_convert(::Type{T}, ω::Complex) where {T<:Real} = convert(Complex{T}, ω) @@ -83,134 +84,213 @@ needs_omega_shift(s::AppliedGreenSolver) = true # If we index with CellIndices, we bypass mortaring and return a bare matrix #region -Base.getindex(g::GreenFunction, i, j = i) = GreenSlice(g, i, j) +## GreenFunction -> GreenSlice + Base.getindex(g::GreenFunction; kw...) = g[siteselector(; kw...)] Base.getindex(g::GreenFunction, kw::NamedTuple) = g[siteselector(; kw...)] +Base.getindex(g::GreenFunction, i, j = i) = GreenSlice(g, i, j) # this calls sites_to_orbs -Base.getindex(g::GreenSolution; kw...) = g[getindex(lattice(g); kw...)] +## GreenSolution -> OrbitalSliceMatrix or AbstractMatrix -# g[::Integer, ::Integer] and g[:, :] - intra and inter contacts +# view(g, ::Integer, ::Integer) and view(g, :, :) - intra and inter contacts Base.view(g::GreenSolution, i::CT, j::CT = i) where {CT<:Union{Integer,Colon}} = view(slicer(g), i, j) -# fastpath for intra and inter-contact -Base.getindex(g::GreenSolution, i::CT, j::CT = i) where {CT<:Union{Integer,Colon}} = - OrbitalSliceMatrix(copy(view(g, i, j)), sites_to_orbs.((i,j), Ref(g))) +# general getindex(::GreenSolution,...) entry point. Relies on GreenSlice constructor +# args could be anything: CellSites, CellOrbs, Colon, Integer, DiagIndices... +function Base.getindex(g::GreenSolution, args...; post = identity, kw...) + gs = getindex(parent(g), args...; kw...) # get GreenSlice + output = getindex!(gs, g; post) + return output +end -# conversion down to CellOrbitals. See sites_to_orbs in slices.jl -Base.getindex(g::GreenSolution, i, j) = getindex(g, sites_to_orbs(i, g), sites_to_orbs(j, g)) -Base.getindex(g::GreenSolution, i) = (i´ = sites_to_orbs(i, g); getindex(g, i´, i´)) +Base.getindex(g::GreenSolution, i::AnyCellOrbitals, j::AnyCellOrbitals) = + slicer(g)[sanitize_cellorbs(i), sanitize_cellorbs(j)] -# wrapped matrix for end user consumption -Base.getindex(g::GreenSolution, i::OrbitalSliceGrouped, j::OrbitalSliceGrouped) = - OrbitalSliceMatrix( - mortar((g[si, sj] for si in cellsdict(i), sj in cellsdict(j))), - (i, j)) +# must ensure that orbindices is not a scalar, to consistently obtain a Matrix +sanitize_cellorbs(c::CellOrbitals) = c +sanitize_cellorbs(c::CellOrbital) = CellOrbitals(cell(c), orbindex(c):orbindex(c)) +sanitize_cellorbs(c::CellOrbitalsGrouped) = CellOrbitals(cell(c), orbindices(c)) -Base.getindex(g::GreenSolution, i::AnyOrbitalSlice, j::AnyOrbitalSlice) = - mortar((g[si, sj] for si in cellsdict(i), sj in cellsdict(j))) +## getindex! - preallocated output for getindex -Base.getindex(g::GreenSolution, i::AnyOrbitalSlice, j::AnyCellOrbitals) = - mortar((g[si, sj] for si in cellsdict(i), sj in (j,))) +# fastpath for intra and inter-contact, only for g::GreenSolution +getindex!(output, g::GreenSolution, i::Union{Integer,Colon}, j::Union{Integer,Colon}; post = identity) = + (parent(output) .= post.(view(g, i, j)); output) -Base.getindex(g::GreenSolution, i::AnyCellOrbitals, j::AnyOrbitalSlice) = - mortar((g[si, sj] for si in (i,), sj in cellsdict(j))) +# indexing over single cells. g can be any type that implements g[::AnyCellOrbitals...] +getindex!(output, g, i::AnyCellOrbitals, j::AnyCellOrbitals; post = identity) = + (output .= post.(g[i, j]); output) -Base.getindex(g::GreenSolution, i::AnyCellOrbitals, j::AnyCellOrbitals) = - slicer(g)[sanitize_cellorbs(i), sanitize_cellorbs(j)] +# indexing over several cells +getindex!(output, g, ci::AnyOrbitalSlice, cj::AnyOrbitalSlice; kw...) = + getindex_cells!(output, g, cellsdict(ci), cellsdict(cj); kw...) +getindex!(output, g, ci::AnyOrbitalSlice, cj::AnyCellOrbitals; kw...) = + getindex_cells!(output, g, cellsdict(ci), (cj,); kw...) +getindex!(output, g, ci::AnyCellOrbitals, cj::AnyOrbitalSlice; kw...) = + getindex_cells!(output, g, (ci,), cellsdict(cj); kw...) + +function getindex_cells!(output, g, cis, cjs; kw...) + for ci in cis, cj in cjs + getindex!(output, g, ci, cj; kw...) # will typically call the method below + end + return output +end + +function getindex!(output::OrbitalSliceMatrix, g, i::AnyCellOrbitals, j::AnyCellOrbitals; kw...) + oi, oj = orbaxes(output) + rows, cols = orbrange(oi, cell(i)), orbrange(oj, cell(j)) + v = view(parent(output), rows, cols) + getindex!(v, g, i, j; kw...) + return output +end + +## common getindex! shortcut in terms of GreenSlice + +# index object g over a slice enconded in gs::GreenSlice, using its preallocated output +getindex!(gs::GreenSlice, g; kw...) = getindex!(call!_output(gs), g, orbinds_or_contactinds(gs)...; kw...) + +# ifelse(rows && cols are contacts, (rows, cols), (orbrows, orbcols)) +# I.e: if rows, cols are contact indices retrieve them instead of orbslices. +orbinds_or_contactinds(g) = orbinds_or_contactinds(rows(g), cols(g), orbrows(g), orbcols(g)) +orbinds_or_contactinds( + r::Union{Colon,Integer,DiagIndices{Colon},DiagIndices{<:Integer}}, + c::Union{Colon,Integer,DiagIndices{Colon},DiagIndices{<:Integer}}, _, _) = (r, c) +orbinds_or_contactinds(_, _, or, oc) = (or, oc) + +## GreenSlicer -> Matrix, dispatches to each solver's implementation # fallback conversion to CellOrbitals -Base.getindex(s::GreenSlicer, i::AnyCellOrbitals, j::AnyCellOrbitals = i) = - getindex(s, sanitize_cellorbs(i), sanitize_cellorbs(j)) +Base.getindex(s::GreenSlicer, i::AnyCellOrbitals, j::AnyCellOrbitals = i; kw...) = + getindex(s, sanitize_cellorbs(i), sanitize_cellorbs(j); kw...) # fallback, for GreenSlicers that forgot to implement getindex -Base.getindex(s::GreenSlicer, ::CellOrbitals, ::CellOrbitals) = +Base.getindex(s::GreenSlicer, ::CellOrbitals, ::CellOrbitals; kw...) = argerror("getindex of $(nameof(typeof(s))) not implemented") # fallback, for GreenSlicers that don't implement view Base.view(g::GreenSlicer, args...) = argerror("GreenSlicer of type $(nameof(typeof(g))) doesn't implement view for these arguments") -# must ensure that orbindices is not a scalar, to consistently obtain a Matrix -sanitize_cellorbs(c::CellOrbitals) = c -sanitize_cellorbs(c::CellOrbital) = CellOrbitals(cell(c), orbindex(c):orbindex(c)) -sanitize_cellorbs(c::CellOrbitalsGrouped) = CellOrbitals(cell(c), orbindices(c)) - -#endregion +## Diagonal and sparse indexing -############################################################################################ -# GreenSolution diagonal indexing -# returns a vector of AbstractMatrices or scalars, -# one per site (if kernel is not missing) or one per orbital (if kernel is missing) -#region +# These i,j are already gone through sites_to_orbs, unless they are Colon or Integer (orbinds_or_contactinds) +getindex!(output, g, i::DiagIndices, j::DiagIndices = i; kw...) = + getindex_diagonal!(output, g, parent(i), kernel(i); kw...) -Base.getindex(gω::GreenSolution{T}, i::DiagIndices, ::DiagIndices = i) where {T} = - append_diagonal!(Complex{T}[], gω, parent(i), kernel(i), gω) |> - maybe_OrbitalSliceArray(sites_to_orbs(i, gω)) +getindex!(output, g, i::OrbitalPairIndices, j::OrbitalPairIndices = i; kw...) = + getindex_sparse!(output, g, parent(i), kernel(i); kw...) #endregion ############################################################################################ -# append_diagonal!(d, x, i, kernel, g; kw...) -# append to d the elements Tr(kernel*x[i,i]) for each site encoded in i, or each orbital -# if kernel is missing. g is the GreenFunction, used to convert i to CellOrbitals +# getindex_diagonal! and getindex_sparse! : GreenSolution diagonal and sparse indexing +# They rely on fill_diagonal! and fill_sparse! to fill the output matrix with mat elements +# We can swith `g` with another object `o` that supports getindex_diag or getindex_sparse +# returning `mat::AbstractMatrix` +# `mat` is not an OrbitalSliceGrouped for performance reasons. Instead the object `o` +# should also implement a `blockstructure` method so we know how to map sites to indices. #region -# If i::Union{SiteSlice,CellSites}, convert to orbitals -append_diagonal!(d, x, i, kernel, g; kw...) = - append_diagonal!(d, x, sites_to_orbs(i, g), kernel; kw...) - -# but not if i is a contact. Jumpt to core driver -append_diagonal!(d, gω::GreenSolution, o::Union{Colon,Integer}, kernel, g; kw...) = - append_diagonal!(d, gω[o,o], contactranges(kernel, o, gω), kernel; kw...) +getindex_diagonal!(output, g::GreenSolution, i::Union{Colon,Integer}, ker; post = identity) = + fill_diagonal!(output, view(g, i, i), contact_kernel_ranges(ker, i, g), ker; post) -# decompose any LatticeSlice into CellOrbitals -append_diagonal!(d, x, s::AnyOrbitalSlice, kernel; kw...) = - append_diagonal!(d, x, cellsdict(s), kernel; kw...) +function getindex_diagonal!(output, g, o::CellOrbitalsGrouped, ker; post = identity) + rngs = orbital_kernel_ranges(ker, o) + mat = getindex_diag(g, o) + fill_diagonal!(output, mat, rngs, ker; post) + return output +end -function append_diagonal!(d, x, s::AnyCellOrbitalsDict, kernel; kw...) - sizehint!(d, length(s)) - for sc in s - append_diagonal!(d, x, sc, kernel; kw...) +# g may be a GreenSolution or any other type with a minimal API (e.g. EigenProduct) +function getindex_diagonal!(output, g, i::OrbitalSliceGrouped, ker; post = identity) + offset = 0 + for o in cellsdict(i) + rngs = orbital_kernel_ranges(ker, o) + mat = getindex_diag(g, o) + fill_diagonal!(output, mat, rngs, ker; post, offset) + offset += length(rngs) end - return d + return output end -append_diagonal!(d, gω::GreenSolution, o::AnyCellOrbitals, kernel; kw...) = - append_diagonal!(d, gω[o,o], orbindranges(kernel, o), kernel; kw...) - -# core driver -function append_diagonal!(d, blockmat::AbstractMatrix, blockrngs, kernel; post = identity) - for rng in blockrngs - val = apply_kernel(kernel, blockmat, rng) - push!(d, post(val)) +function getindex_sparse!(output, g, hmask::Hamiltonian, ker; post = identity) + bs = blockstructure(g) # used to find blocks of g for nonzeros in hmask + oj = sites_to_orbs(sites(zerocell(hmask), :), bs) # full cell + for har in harmonics(hmask) + oi = CellIndices(dcell(har), oj) + mat_cell = getindex_sparse(g, oi, oj, har) + fill_sparse!(har, mat_cell, bs, ker; post) end - return d + fill_sparse!(parent(output), hmask) + return output +end + +# non-optimized fallback, can be specialized by solvers +# optimized version could compute only site blocks in diagonal +getindex_diag(g, oi) = g[oi, oi] +# optimized version could compute only nonzero blocks from har::Harmonic +getindex_sparse(g, oi, oj, har) = g[oi, oj] + +# input index ranges for each output index to fill +contact_kernel_ranges(::Missing, o::Colon, gω) = 1:norbitals(contactorbitals(gω)) +contact_kernel_ranges(::Missing, i::Integer, gω) = 1:norbitals(contactorbitals(gω), i) +contact_kernel_ranges(kernel, ::Colon, gω) = orbranges(contactorbitals(gω)) +contact_kernel_ranges(kernel, i::Integer, gω) = orbranges(contactorbitals(gω), i) +orbital_kernel_ranges(::Missing, o::CellOrbitalsGrouped) = eachindex(orbindices(o)) +orbital_kernel_ranges(kernel, o::CellOrbitalsGrouped) = orbranges(o) + +## core drivers +# write to the diagonal of output, with a given offset, the elements Tr(kernel*mat_cell[i,i]) +# for each orbital or site encoded in range i +function fill_diagonal!(output, mat_cell, blockrngs, kernel; post = identity, offset = 0) + for (n, rng) in enumerate(blockrngs) + i = n + offset + output[i, i] = post(apply_kernel(kernel, mat_cell, rng)) + end + return output +end + +# overwrite nonzero h elements with Tr(kernel*mat_cell[irng,jrng]) for each site pair i,j +function fill_sparse!(har::Harmonic, mat_cell, bs, kernel; post = identity) + B = blocktype(har) + hmat = unflat(har) + rows = rowvals(hmat) + nzs = nonzeros(hmat) + for col in axes(hmat, 2), ptr in nzrange(hmat, col) + row = rows[ptr] + orows, ocols = flatrange(bs, row), flatrange(bs, col) + val = post(apply_kernel(kernel, mat_cell, orows, ocols)) + nzs[ptr] = mask_block(B, val) + end + needs_flat_sync!(matrix(har)) + return har +end + +# overwrite sparse output with an vcat of all Harmonics in h +function fill_sparse!(output::SparseMatrixCSC, h::Hamiltonian) + out_nzs = nonzeros(output) + for col in axes(output, 2) + out_ptr_rng = nzrange(output, col) + iptr = 1 + for har in harmonics(h) + hmat = flat(har) + hnzs = nonzeros(hmat) + for ptr in nzrange(hmat, col) + out_nzs[out_ptr_rng[iptr]] = hnzs[ptr] + iptr += 1 + end + end + end + return output end -# If no kernel is provided, we return the whole diagonal -contactranges(::Missing, o::Colon, gω) = 1:norbitals(contactorbitals(gω)) -contactranges(::Missing, i::Integer, gω) = 1:norbitals(contactorbitals(gω), i) -contactranges(kernel, ::Colon, gω) = orbranges(contactorbitals(gω)) -contactranges(kernel, i::Integer, gω) = orbranges(contactorbitals(gω), i) - -orbindranges(::Missing, o) = eachindex(orbindices(o)) -orbindranges(kernel, o) = orbranges(o) - -apply_kernel(kernel, gblock, rng) = apply_kernel(kernel, view_or_scalar(gblock, rng)) - -apply_kernel(kernel::Missing, v::Number) = v -apply_kernel(kernel::AbstractMatrix, v) = tr(kernel * v) -apply_kernel(kernel::UniformScaling, v) = kernel.λ * tr(v) -apply_kernel(kernel::Number, v) = kernel * tr(v) -apply_kernel(kernel::Diagonal, v::AbstractMatrix) = sum(i -> kernel[i] * v[i, i], eachindex(kernel)) -apply_kernel(kernel::Diagonal, v::Number) = only(kernel) * v - -view_or_scalar(gblock, rng::UnitRange) = view(gblock, rng, rng) -view_or_scalar(gblock, orb::Integer) = gblock[orb, orb] +apply_kernel(kernel, mat_cell, rows, cols = rows) = apply_kernel(kernel, view_or_scalar(mat_cell, rows, cols)) +apply_kernel(kernel::Missing, v) = v +apply_kernel(kernel, v) = trace_prod(kernel, v) -maybe_scalarize(s::OrbitalSliceGrouped, kernel::Missing) = s -maybe_scalarize(s::OrbitalSliceGrouped, kernel) = scalarize(s) +view_or_scalar(mat_cell, rows::UnitRange, cols::UnitRange) = view(mat_cell, rows, cols) +view_or_scalar(mat_cell, orb::Integer, orb´::Integer) = mat_cell[orb, orb´] #endregion diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 4036b0e9..3e8c0e01 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -148,8 +148,10 @@ hamiltonian(lat::Lattice, m::Interblock{<:ParametricModel}; kw...) = parametric( hamiltonian(lat, basemodel(parent(m)), block(m); kw...), modifier.(terms(parent(m)))...) -function hamiltonian(lat::Lattice{T}, m::TightbindingModel = TightbindingModel(), block = missing; orbitals = Val(1)) where {T} - b = IJVBuilder(lat, orbitals) +hamiltonian(lat::Lattice, m::TightbindingModel = TightbindingModel(), block = missing; orbitals = Val(1)) = + hamiltonian!(IJVBuilder(lat, orbitals), m, block) + +function hamiltonian!(b::IJVBuilder, m::TightbindingModel, block = missing) add!(b, m, block) return hamiltonian(b) end @@ -452,7 +454,7 @@ Base.getindex(h::AbstractHamiltonian{<:Any,<:Any,L}, ::HybridInds{Tuple{}}) wher Base.getindex(h::ParametricHamiltonian{<:Any,<:Any,L}, dn::HybridInds{SVector{L,Int}}) where {L} = getindex(hamiltonian(h), dn) -function Base.getindex(h::Hamiltonian{<:Any,<:Any,L}, dn::HybridInds{SVector{L,Int}}) where {L} +Base.@propagate_inbounds 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 @@ -481,6 +483,9 @@ Base.getindex(h::Hamiltonian, i, j) = getindex(h, sites_to_orbs(i, h), sites_to_ # we need AbstractHamiltonian here to avoid ambiguities with dn above Base.getindex(h::AbstractHamiltonian, i) = (i´ = sites_to_orbs(i, h); getindex(h, i´, i´)) +Base.getindex(h::AbstractHamiltonian, ::SparseIndices...) = + argerror("Sparse indexing not yet supported for AbstractHamiltonian") + # wrapped matrix for end user consumption Base.getindex(h::Hamiltonian, i::OrbitalSliceGrouped, j::OrbitalSliceGrouped) = OrbitalSliceMatrix( diff --git a/src/meanfield.jl b/src/meanfield.jl index 7f9fa7c5..6428c455 100644 --- a/src/meanfield.jl +++ b/src/meanfield.jl @@ -1,8 +1,84 @@ -############################################################################################ -# hartree constructors -#region +# ############################################################################################ +# # meanfield +# # designed to construct hfield and ffield, such that +# # hfield[i] = ν * Σ_k v_H(r_i-r_k) * tr(ρ[k,k]*Q) +# # ffield[i,j] = v_F(r_i-r_j) * Q * ρ[i,j] * Q +# # where ν = ifelse(nambu, 1/2, 1), and Q is the charge matrix or [q 0; 0 -q] if nambu. +# # we precompute v_H^{ik} = \sum_n v_H(r_{i0} - r_{kn}), exploiting ρ translation symmetry +# #region -# function hartree +# struct MeanField{T,O<:OrbitalSliceMatrix{T},H,F,Q,B} +# Vhmat::SparseMatrixCSC{T,Int} +# Vfmat::O +# hrho::H +# frho::F +# charge::Q +# blocktype::B +# end +# struct EvaluatedMeanField{H<:OrbitalSliceVector,F<:OrbitalSliceMatrix} +# hfield::H +# ffield::F +# end -#endregion +# #region ## Constructors ## + +# function meanfield(g::GreenFunction{T}, args...; +# potential = Returns(0), hartreepotential = potential, fockpotential = hartreepotential, +# charge = I, nambu::Bool = false, selector = (; range = 0), kw...) where {T} + +# Vh = sanitize_potential(hartreepotential) +# Vf = sanitize_potential(fockpotential) + +# B = blocktype(hamiltonian(g)) +# lat = lattice(hamiltonian(g)) +# hFock = lat |> hopping((r, dr) -> Vf(dr); selector..., includeonsite = true) +# hHartree = Vh === Vf ? hFock : +# lat |> hopping((r, dr) -> Vh(dr); selector..., includeonsite = true) + +# ldst, lsrc = ham_to_latslices(hFock) +# odst, osrc = sites_to_orbs(ldst, g), sites_to_orbs(lsrc, g) +# ρhartree = densitymatrix(g[diagonal(osrc, kernel = charge)], args...; kw...) +# ρfock = densitymatrix(g[odst, osrc], args...; kw...) + +# Vhmat = sum(unflat, harmonics(hHartree)) +# nambu && (nonzeros(Vhmat) .*= T(1/2)) +# Vfmat = hFock[ldst, lsrc] + +# return MeanField(Vhmat, Vfmat, ρhartree, ρfock, charge, B) +# end + +# sanitize_potential(x::Number) = Returns(x) +# sanitize_potential(x::Function) = x + +# sanitize_charge(x, nambu) = nambu ? SA[x 0*x; 0*x -x] : x + +# #endregion + +# #region ## API ## + + +# function (m::MeanField{T})(args...; params...) where {T} +# Q, B = m.charge, m.blocktype +# trρQ = m.hrho(args...; params...) +# QvtrρQ = (m.Vhmat * parent(trρQ)) .* Ref(Q) +# hfield = OrbitalSliceVector(QvtrρQ, orbaxes(trρQ)) +# ffield = m.frho(args...; params...) +# if !isempty(ffield) +# oi, oj = orbaxes(ffield) +# for sj in sites(oj, B), si in sites(oi, B) # si, sj are CellSitePos +# vij = T(only(view(m.Vfmat, CellSite(si), CellSite(sj)))) # a scalar +# view(ffield, si, sj) .= vij * (Q * ffield[si, sj] * Q) +# end +# end +# return EvaluatedMeanField(hfield, ffield) +# end + +# Base.view(m::EvaluatedMeanField, i) = view(m.hfield, i) +# Base.view(m::EvaluatedMeanField, i, j) = view(m.ffield, i, j) + +# Base.getindex(m::EvaluatedMeanField, i) = m.hfield[i] +# Base.getindex(m::EvaluatedMeanField, i, j) = m.ffield[i, j] +# #endregion + +# #endregion diff --git a/src/observables.jl b/src/observables.jl index a88d5cb4..c5b472a4 100644 --- a/src/observables.jl +++ b/src/observables.jl @@ -145,7 +145,7 @@ options(I::Integrator) = I.quadgk_opts # scalar version function call!(I::Integrator{<:Any,Missing}; params...) fx = x -> begin - y = call!(I.integrand, x; params...) + y = call!(I.integrand, x; params...) # should be a scalar I.callback(x, y) return y end @@ -157,12 +157,13 @@ end # nonscalar version function call!(I::Integrator; params...) fx! = (y, x) -> begin - y .= call!(I.integrand, x; params...) + y .= serialize(call!(I.integrand, x; params...)) I.callback(x, y) return y end - result, err = quadgk!(fx!, I.result, I.points...; I.quadgk_opts...) - result´ = I.post(result) # note: post-processing is not element-wise & can be in-place + result, err = quadgk!(fx!, serialize(I.result), I.points...; I.quadgk_opts...) + # note: post-processing is not element-wise & can be in-place + result´ = deserialize(I.result, I.post(result)) return result´ end @@ -173,8 +174,8 @@ end ############################################################################################ # ldos: local spectral density -# d = ldos(::GreenSolution; kernel = I) -> d[sites...]::Vector -# d = ldos(::GreenSlice; kernel = I) -> d(ω; params...)::Vector +# d = ldos(::GreenSolution; kernel = missing) -> d[sites...]::Vector +# d = ldos(::GreenSlice; kernel = missing) -> d(ω; params...)::Vector # Here ldos is given as Tr(ρᵢᵢ * kernel) where ρᵢᵢ is the spectral function at site i # Here is the generic fallback that uses G. Any more specialized methods need to be added # to each GreenSolver @@ -182,24 +183,32 @@ end struct LocalSpectralDensitySolution{T,E,L,G<:GreenSolution{T,E,L},K} <: IndexableObservable gω::G - kernel::K # should return a float when applied to gω[CellSite(n,i)] + kernel::K end struct LocalSpectralDensitySlice{T,E,L,G<:GreenSlice{T,E,L},K} gs::G - kernel::K # should return a float when applied to gω[CellSite(n,i)] + kernel::K # also inside gs end #region ## Constructors ## -ldos(gω::GreenSolution; kernel = I) = LocalSpectralDensitySolution(gω, kernel) +ldos(gω::GreenSolution; kernel = missing) = LocalSpectralDensitySolution(gω, kernel) -function ldos(gs::GreenSlice{T}; kernel = I) where {T} +function ldos(gs::GreenSlice{T}; kernel = missing) where {T} rows(gs) === cols(gs) || argerror("Cannot take ldos of a GreenSlice with rows !== cols") - return LocalSpectralDensitySlice(gs, kernel) + g = parent(gs) + i = ensure_diag_axes(rows(gs), kernel) + gs´ = GreenSlice(g, i, i, T) # forces output eltype to be T<:Real + return LocalSpectralDensitySlice(gs´, kernel) end +ensure_diag_axes(inds, kernel) = diagonal(inds; kernel) +ensure_diag_axes(inds::DiagIndices, kernel) = diagonal(parent(inds); kernel) +ensure_diag_axes(::SparseIndices, _) = + argerror("Unexpected sparse indices in GreenSlice") + #endregion #region ## API ## @@ -208,27 +217,22 @@ greenfunction(d::LocalSpectralDensitySlice) = d.g kernel(d::Union{LocalSpectralDensitySolution,LocalSpectralDensitySlice}) = d.kernel -ldos_kernel(g, kernel::UniformScaling) = - kernel.λ * imag(tr(g)) / π -ldos_kernel(g, kernel) = -imag(tr(g * kernel)) / π - Base.getindex(d::LocalSpectralDensitySolution; selectors...) = - d[getindex(lattice(d.gω); selectors...)] - -Base.getindex(d::LocalSpectralDensitySolution, s::SiteSelector) = - d[getindex(lattice(d.gω), s)] + d[siteselector(; selectors...)] -Base.getindex(d::LocalSpectralDensitySolution{T}, i) where {T} = - append_diagonal!(T[], d.gω, i, d.kernel, d.gω; post = x -> -imag(x)/π) |> - maybe_OrbitalSliceArray(sites_to_orbs(i, d.gω)) # see greenfunction.jl +function Base.getindex(d::LocalSpectralDensitySolution{T}, i) where {T} + di = diagonal(i, kernel = d.kernel) + gs = GreenSlice(parent(d.gω), di, di, T) # get GreenSlice with real output + output = getindex!(gs, d.gω; post = x -> -imag(x)/π) + return diag(output) +end (d::LocalSpectralDensitySlice)(ω; params...) = copy(call!(d, ω; params...)) # fallback through LocalSpectralDensitySolution - overload to allow a more efficient path function call!(d::LocalSpectralDensitySlice{T}, ω; params...) where {T} - orbslice_or_contact = first(orbinds_or_contactinds(d.gs)) - gω = call!(parent(d.gs), ω; params...) - l = ldos(gω; kernel = d.kernel) - return l[orbslice_or_contact] + output = call!(d.gs, ω; post = x -> -imag(x)/π, params...) + return diag(output) end #endregion @@ -245,6 +249,9 @@ end # We use a default charge = -I, corresponding to normal currents densities in units of e/h #region +## TODO: this could probably be refactored using sparse indexing +## removing GreenSolutionCache and unflat_sparse_slice in the process + struct CurrentDensitySolution{T,E,L,G<:GreenSolution{T,E,L},K,V<:Union{Missing,SVector}} <: IndexableObservable gω::G charge::K # should return a float when traced with gʳᵢⱼHᵢⱼ @@ -483,12 +490,13 @@ end # generic fallback (for other solvers) (ρ::DensityMatrix)(mu = 0, kBT = 0; params...) = - ρ.solver(mu, kBT; params...) |> maybe_OrbitalSliceArray(axes(ρ.gs)) + ρ.solver(mu, kBT; params...) # special case for integrator solver -(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0, override = missing; params...) = - ρ.solver(mu, kBT, override)(; params...) |> maybe_OrbitalSliceArray(axes(ρ.gs)) +(ρ::DensityMatrix{<:DensityMatrixIntegratorSolver})(mu = 0, kBT = 0, override_path = missing; params...) = + ρ.solver(mu, kBT, override_path)(; params...) -(s::DensityMatrixIntegratorSolver)(mu, kBT, override = missing) = s.ifunc(mu, kBT, override); +(s::DensityMatrixIntegratorSolver)(mu, kBT, override_path = missing) = + s.ifunc(mu, kBT, override_path); # redirects to specialized method densitymatrix(gs::GreenSlice; kw...) = @@ -503,11 +511,11 @@ densitymatrix(gs::GreenSlice, ωmax::Number; opts...) = densitymatrix(gs, (-ωma function densitymatrix(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), imshift = missing, atol = 1e-7, opts...) where {T} # check_nodiag_axes(gs) - result = similar_Array(gs) + result = copy(call!_output(gs)) opts´ = (; imshift, slope = 1, post = gf_to_rho!, atol, opts...) ωpoints_vec = collect(promote_type(T, typeof.(ωpoints)...), ωpoints) - function ifunc(mu, kBT, override) - ωpoints´ = override_path!(override, ωpoints_vec, ωpoints) + function ifunc(mu, kBT, override_path) + ωpoints´ = override_path!(override_path, ωpoints_vec, ωpoints) ρd = DensityMatrixIntegrand(gs, T(mu), T(kBT), omegamap) pts = maybe_insert_mu!(ωpoints_vec, ωpoints´, mu, kBT) return Integrator(result, ρd, pts; opts´...) @@ -548,7 +556,7 @@ function override_path!(pts´, ptsvec::Vector{<:Complex}, pts) return pts´ end -override_path!(override, ptsvec, pts) = +override_path!(override_path, ptsvec, pts) = argerror("Override of real ωpoints not supported, use complex ωpoints upon construction") #endregion @@ -558,10 +566,9 @@ override_path!(override, ptsvec, pts) = (gf::DensityMatrixIntegrand)(ω; params...) = copy(call!(gf, ω; params...)) function call!(gf::DensityMatrixIntegrand, ω; params...) - gω = call!(gf.gs, ω; gf.omegamap(ω)..., params...) - f = fermi(ω - gf.mu, inv(gf.kBT)) - gω .*= f - return gω + output = call!(gf.gs, ω; gf.omegamap(ω)..., params...) + serialize(output) .*= fermi(ω - gf.mu, inv(gf.kBT)) + return output end function gf_to_rho!(x::AbstractMatrix) @@ -585,6 +592,10 @@ temperature(D::DensityMatrixIntegrand) = D.kBT chemicalpotential(D::DensityMatrixIntegrand) = D.mu +Base.parent(ρ::DensityMatrix) = ρ.gs + +call!_output(ρ::DensityMatrix) = call!_output(ρ.gs) + #endregion ############################################################################################ @@ -631,10 +642,10 @@ end # generic fallback (for other solvers) (j::Josephson)(kBT = 0; params...) = j.solver(kBT; params...) # special case for integrator solver -(j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0, override = missing; params...) = - j.solver(kBT, override)(; params...) +(j::Josephson{<:JosephsonIntegratorSolver})(kBT = 0, override_path = missing; params...) = + j.solver(kBT, override_path)(; params...) -(s::JosephsonIntegratorSolver)(kBT, override = missing) = s.ifunc(kBT, override) +(s::JosephsonIntegratorSolver)(kBT, override_path = missing) = s.ifunc(kBT, override_path) josephson(gs::GreenSlice{T}, ωmax::Number; kw...) where {T} = josephson(gs, (-ωmax, ωmax); kw...) @@ -651,8 +662,8 @@ function josephson(gs::GreenSlice{T}, ωpoints; omegamap = Returns((;)), phases phases´, traces = sanitize_phases_traces(phases, T) opts´ = (; imshift, slope = 1, post = real, atol, opts...) ωpoints_vec = collect(promote_type(T, typeof.(ωpoints)...), ωpoints) - function ifunc(kBT, override) - ωpoints´ = override_path!(override, ωpoints_vec, ωpoints) + function ifunc(kBT, override_path) + ωpoints´ = override_path!(override_path, ωpoints_vec, ωpoints) jd = JosephsonIntegrand(g, T(kBT), contact, tauz, phases´, omegamap, traces, Σfull, Σ, similar(Σ), similar(Σ), similar(Σ), similar(tauz, Complex{T})) pts = maybe_insert_mu!(ωpoints_vec, ωpoints´, zero(T), kBT) @@ -700,7 +711,7 @@ function call!(J::JosephsonIntegrand, ω; params...) end function josephson_traces(J, gω, f) - gr = gω[J.contactind, J.contactind] + gr = view(gω, J.contactind, J.contactind) Σi = selfenergy!(J.Σ, gω, J.contactind) return josephson_traces!(J, gr, Σi, f) end diff --git a/src/serializer.jl b/src/serializer.jl index 4df2dad7..b25d6c45 100644 --- a/src/serializer.jl +++ b/src/serializer.jl @@ -1,4 +1,3 @@ - ############################################################################################ # AppliedSerializer # support for serialization of Hamiltonians and ParametricHamiltonians @@ -6,7 +5,7 @@ ## serialize and serialize! -function serialize!(v, s::AppliedSerializer; kw...) +Base.@propagate_inbounds function serialize!(v, s::AppliedSerializer; kw...) @boundscheck check_serializer_length(s, v) h0 = parent_hamiltonian(s) h = call!(h0; kw...) @@ -60,7 +59,7 @@ deserialize(s::AppliedSerializer, v; kw...) = deserialize!(s(; kw...), v) deserialize!(s::AppliedSerializer, v; kw...) = deserialize!(call!(s; kw...), v) -function deserialize!(s::AppliedSerializer{<:Any,<:Hamiltonian}, v) +Base.@propagate_inbounds function deserialize!(s::AppliedSerializer{<:Any,<:Hamiltonian}, v) @boundscheck check_serializer_length(s, v) h = parent_hamiltonian(s) dec = decoder(s) @@ -126,3 +125,27 @@ function check(s::AppliedSerializer{T}; params...) where {T} end #endregion + + +############################################################################################ +# serialize and deserialize for OrbitalSliceArray +# extract underlying data and reconstruct an object using such data +# no check of correct array dimensionality is performed for performance +#region + +serialize(a::OrbitalSliceArray) = serialize_array(parent(a)) +serialize(a::AbstractArray) = serialize_array(a) + +serialize_array(a::Array) = a +serialize_array(a::Diagonal{<:Any,<:Vector}) = a.diag +serialize_array(a::SparseMatrixCSC) = nonzeros(a) + +deserialize(a::OrbitalSliceArray, v) = OrbitalSliceArray(deserialize_array(parent(a), v), orbaxes(a)) +deserialize(a::AbstractArray, v) = deserialize_array(a, v) + +deserialize_array(::A, v::A´) where {N,T,T´,A<:Array{T,N},A´<:Array{T´,N}} = v +deserialize_array(::Diagonal, v::Vector) = Diagonal(v) +deserialize_array(a::SparseMatrixCSC, v::Vector) = SparseMatrixCSC(a.m, a.n, a.colptr, a.rowval, v) + + +#endregion diff --git a/src/show.jl b/src/show.jl index a0544538..0d861e65 100644 --- a/src/show.jl +++ b/src/show.jl @@ -57,7 +57,8 @@ displayname(s::Sublat) = sublatname(s) == Symbol(:_) ? "pending" : string(":", s function Base.show(io::IO, c::CellSites) i = get(io, :indent, "") - print(io, i, summary(c), "\n") + print(io, i, summary(c), "\n", +"$i Sites : $(siteindices(c))") end Base.summary(c::CellSites{L,V}) where {L,V} = @@ -195,7 +196,7 @@ 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 Orbitals : $(norbitals(h, :)) $i Element type : $(displaytype(blocktype(h))) $i Onsites : $(nonsites(h)) $i Hoppings : $(nhoppings(h)) @@ -378,6 +379,23 @@ Base.summary(g::GreenSlice{T,E,L}) where {T,E,L} = #endregion +############################################################################################ +# SparseIndices +#region + +function Base.show(io::IO, s::SparseIndices) + i = get(io, :indent, "") + print(io, i, summary(s)) +end + +Base.summary(::PairIndices) = + "PairIndices: a form of SparseIndices used to sparsely index an object such as a GreenFunction on a selection of site pairs" + +Base.summary(::DiagIndices) = + "DiagIndices: a form of SparseIndices used to index the diagonal of an object such as a GreenFunction" + +#endregion + ############################################################################################ # Conductance #region @@ -527,14 +545,14 @@ Base.summary(::CurrentDensitySlice{T}) where {T} = # For simplified printing of the array typename -function Base.showarg(io::IO, ::OrbitalSliceMatrix{<:Any,M}, toplevel) where {M} +function Base.showarg(io::IO, ::OrbitalSliceMatrix{T,M}, toplevel) where {T,M} toplevel || print(io, "::") - print(io, "OrbitalSliceMatrix{$M}") + print(io, "OrbitalSliceMatrix{$T,$M}") end -function Base.showarg(io::IO, ::OrbitalSliceVector{<:Any,M}, toplevel) where {M} +function Base.showarg(io::IO, ::OrbitalSliceVector{T,M}, toplevel) where {T,M} toplevel || print(io, "::") - print(io, "OrbitalSliceVector{$M}") + print(io, "OrbitalSliceVector{$T,$M}") end #endregion diff --git a/src/slices.jl b/src/slices.jl index aa3254ac..9d5f84fe 100644 --- a/src/slices.jl +++ b/src/slices.jl @@ -47,7 +47,7 @@ Base.getindex(ls::LatticeSlice, ss::SiteSelector, args...) = getindex(ls, apply(ss, ls), args...) # return cell, siteindex of the i-th site of LatticeSlice -function Base.getindex(l::LatticeSlice, i::Integer) +Base.@propagate_inbounds function Base.getindex(l::LatticeSlice, i::Integer) offset = 0 for scell in cellsdict(l) ninds = length(siteindices(scell)) @@ -148,7 +148,7 @@ function unflat_sparse_slice(h::Hamiltonian, lsrows::LS, lscols::LS = lsrows, po # @assert lattice(h) === lattice(ls) # TODO: fails upon plotting a current density (see tutorial) cszero = zerocellsites(h, 1) # like `sites(1)`, but with explicit cell B = typeof(post(zero(blocktype(h)), (cszero, cszero))) - nrows, ncols = length(lsrows), length(lscols) + nrows, ncols = nsites(lsrows), nsites(lscols) builder = CSC{B}(ncols) hars = harmonics(h) for colcs in cellsites(lscols) @@ -185,6 +185,24 @@ end #endregion + +############################################################################################ +# ham_to_orbslices: build OrbitalSliceGrouped for source and destination cells of h +#region + +function ham_to_orbslices(h::Hamiltonian) + lat = lattice(h) + cssrc = CellSites(zerocell(lat), :) + osrc = sites_to_orbs(LatticeSlice(lat, [cssrc]), h) + csdst = [CellSites(dcell(har), :) for har in harmonics(h)] + odst = sites_to_orbs(LatticeSlice(lat, csdst), h) + return odst, osrc +end + +ham_to_orbslices(h::ParametricHamiltonian) = ham_to_orbslices(parent(h)) + +#endregion + ############################################################################################ # combine, combine! and intersect! #region @@ -336,19 +354,39 @@ end sites_to_orbs(s::AnyOrbitalSlice, _) = s sites_to_orbs(c::AnyCellOrbitalsDict, _) = c sites_to_orbs(c::AnyCellOrbitals, _) = c +sites_to_orbs(c::SparseIndices{<:Any,<:Hamiltonian}, _) = c # unused # sites_to_orbs_nogroups(cs::CellOrbitals, _) = cs -## DiagIndices +## DiagIndices -> DiagIndices -sites_to_orbs(d::DiagIndices, g) = DiagIndices(sites_to_orbs(parent(d), g), kernel(d)) +function sites_to_orbs(d::DiagIndices, g) + ker = kernel(d) + inds = sites_to_orbs(parent(d), g) + return SparseIndices(inds, ker) +end + +## PairIndices -> OrbitalPairIndices + +# creates a Hamiltonian with same blocktype as g, or complex if kernel::Missing +# this may be used as an intermediate to build sparse versions of g[i,j] +function sites_to_orbs(d::PairIndices{<:HopSelector}, g) + hg = hamiltonian(g) + ker = kernel(d) + bs = maybe_scalarize(blockstructure(hg), ker) + b = IJVBuilder(lattice(hg), bs) + hopsel = parent(d) + h = hamiltonian!(b, hopping(I, hopsel)) + return SparseIndices(h, ker) # OrbitalPairIndices +end ## convert SiteSlice -> OrbitalSliceGrouped sites_to_orbs(kw::NamedTuple, g) = sites_to_orbs(siteselector(; kw...), g) sites_to_orbs(s::SiteSelector, g) = sites_to_orbs(lattice(g)[s], g) -sites_to_orbs(i::Union{Colon,Integer}, g) = orbslice(contacts(g), i) +sites_to_orbs(i::Union{Colon,Integer}, g) = + orbslice(contacts(g), i) sites_to_orbs(l::SiteSlice, g) = OrbitalSliceGrouped(lattice(l), sites_to_orbs(cellsdict(l), blockstructure(g))) sites_to_orbs(l::SiteSlice, s::OrbitalSliceGrouped) = s[l] @@ -365,12 +403,14 @@ end ## convert CellSites -> CellOrbitalsGrouped -sites_to_orbs(c::CellSites, g) = sites_to_orbs(sanitize_cellindices(c, g), blockstructure(g)) +sites_to_orbs(c::CellSites, g, ker...) = + sites_to_orbs(sanitize_cellindices(c, g), blockstructure(g), ker...) -function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure) +function sites_to_orbs(cs::CellSites, os::OrbitalBlockStructure, ker...) + os´ = maybe_scalarize(os, ker...) sites = siteindices(cs) - groups = _groups(sites, os) # sites, orbranges - orbinds = _orbinds(sites, groups, os) + groups = _groups(sites, os´) # sites, orbranges + orbinds = _orbinds(sites, groups, os´) return CellOrbitalsGrouped(cell(cs), orbinds, Dictionary(groups...)) end @@ -393,7 +433,7 @@ function sites_to_orbs_nogroups(cellsdict::CellSitesDict{L}, os::OrbitalBlockStr return cellinds_to_dict(co) end -## convert CellSites -> CellOrbitals +## convert CellSites -> CellOrbitals (no groups) sites_to_orbs_nogroups(cs::CellSites, g) = sites_to_orbs_nogroups(sanitize_cellindices(cs, g), blockstructure(g)) @@ -404,7 +444,7 @@ function sites_to_orbs_nogroups(cs::CellSites, os::OrbitalBlockStructure) return CellOrbitals(cell(cs), orbinds) end -## convert CellOrbitalsGrouped -> CellOrbitals +## convert CellOrbitalsGrouped -> CellOrbitals (no groups) # unused # sites_to_orbs_nogroups(cs::CellOrbitalsGrouped, _) = CellOrbitals(cell(cs), orbindices(cs)) diff --git a/src/solvers/green/bands.jl b/src/solvers/green/bands.jl index 528dea55..4ee555ce 100644 --- a/src/solvers/green/bands.jl +++ b/src/solvers/green/bands.jl @@ -734,21 +734,30 @@ Base.getindex(s::BandsGreenSlicer{<:Any,Missing}, i::CellOrbitals, j::CellOrbita Base.getindex(s::BandsGreenSlicer, i::CellOrbitals, j::CellOrbitals) = semi_band_slice(s, i, j) -function inf_band_slice(s::BandsGreenSlicer{C}, i::CellOrbitals, j::CellOrbitals) where {C} +function inf_band_slice(s::BandsGreenSlicer{C}, i::Union{SparseIndices,CellOrbitals}, j::Union{SparseIndices,CellOrbitals}) where {C} solver = s.solver gmat = zeros(C, norbitals(i), norbitals(j)) inf_band_slice!(gmat, s.ω, (i, j), solver.subband, solver.subbandsimps) return gmat end -function inf_band_slice!(gmat, ω, (i, j)::Tuple{CellOrbitals,CellOrbitals}, +function inf_band_slice!(gmat, ω, (i, j)::Tuple{Union{SparseIndices,CellOrbitals},Union{SparseIndices,CellOrbitals}}, subband::Subband, subbandsimps::SubbandSimplices, simpinds = eachindex(simplices(subband))) - dist = cell(i) - cell(j) - orbs = orbindices(i), orbindices(j) + dist = dist_or_zero(i, j) + orbs = orbstuple_or_orbranges(i, j) return inf_band_slice!(gmat, ω, dist, orbs, subband, subbandsimps, simpinds) end +## Fast codepath for diagonal slices +dist_or_zero(i, j) = cell(i) - cell(j) +dist_or_zero(i::S, j::S) where {S<:SparseIndices{<:AnyCellOrbitals}} = zero(cell(parent(i))) + +# can be all orbitals in i⊗j, or a range of orbitals for each site along diagonal +orbstuple_or_orbranges(i, j) = orbindices(i), orbindices(j) +orbstuple_or_orbranges(i::S, j::S) where {S<:SparseIndices{<:AnyCellOrbitals}} = + orbranges(parent(i)) + # main driver function inf_band_slice!(gmat, ω, dist, orbs, subband, subbandsimps, simpinds) ψPdict = projected_states(subband) @@ -849,7 +858,7 @@ function muladd_ψPψ⁺!(gmat, α, ψ, (rows, cols)::Tuple) return gmat end -# fallback for a generic rngs collection (used e.g. by rngs = orbs in append_diagonal! below) +# fill in only rngs blocks in diagonal of gmat (used by fast codepath for diagonal indexing) function muladd_ψPψ⁺!(gmat, α, ψ, rngs) for rng in rngs if length(rng) == 1 @@ -873,28 +882,12 @@ minimal_callsafe_copy(s::BandsGreenSlicer, parentham, parentcontacts) = s # it #endregion ############################################################################################ -# append_diagonal! -# optimized block diagonal of g[::CellOrbitals] for BandsGreenSlicer without boundary -# otherwise we need to do it the normal way +# getindex_diag: optimized cell indexing for BandsGreenSlicer, used by diagonal indexing #region -function append_diagonal!(d, gω::GreenSolution{T,<:Any,<:Any,G}, o::AnyCellOrbitals, kernel; kw...) where {T,G<:BandsGreenSlicer{<:Any,Missing}} - s = slicer(gω) # BandsGreenSlicer - norbs = flatsize(gω) - rngs = orbranges(o) - dist = zero(cell(o)) - C = complex(T) - gmat = Matrix{C}(undef, norbs, norbs) - for rng in rngs - gmat[rng, rng] .= zero(C) - end - subband, subbandsimps = s.solver.subband, s.solver.subbandsimps - simpinds = eachindex(simplices(subband)) - # to main driver with orbs = rngs - inf_band_slice!(gmat, s.ω, dist, rngs, subband, subbandsimps, simpinds) - orbs = orbindranges(kernel, o) - return append_diagonal!(d, gmat, orbs, kernel; kw...) -end +# triggers fast codepath above +getindex_diag(gω::GreenSolution{T,<:Any,<:Any,G}, o::CellOrbitalsGrouped) where {T,G<:BandsGreenSlicer{<:Any,Missing}} = + inf_band_slice(slicer(gω), SparseIndices(o, Missing), SparseIndices(o, Missing)) #endregion diff --git a/src/solvers/green/kpm.jl b/src/solvers/green/kpm.jl index 1996d7ad..4b825c49 100644 --- a/src/solvers/green/kpm.jl +++ b/src/solvers/green/kpm.jl @@ -236,7 +236,8 @@ end # specialized DensityMatrix method for GS.KPM #region -struct DensityMatrixKPMSolver{T,M} +struct DensityMatrixKPMSolver{T,M,G<:GreenSlice} + gs::G momenta::Vector{M} bandCH::Tuple{T,T} end @@ -246,7 +247,7 @@ end function densitymatrix(s::AppliedKPMGreenSolver, gs::GreenSlice{T}) where {T} has_selfenergy(gs) && argerror("The KPM densitymatrix solver currently support only `nothing` contacts") momenta = slice_momenta(s.momenta, gs) - solver = DensityMatrixKPMSolver(momenta, s.bandCH) + solver = DensityMatrixKPMSolver(gs, momenta, s.bandCH) return DensityMatrix(solver, gs) end @@ -266,8 +267,11 @@ _maybe_diagonal(blockmat, gs) = _maybe_diagonal(blockmat, gs, rows(gs)) _maybe_diagonal(blockmat, gs, is) = blockmat function _maybe_diagonal(blockmat::AbstractArray{C}, gs, is::DiagIndices) where {C} - blockrngs = contactranges(kernel(is), parent(is), gs) # see greenfunction.jl - return append_diagonal!(C[], blockmat, blockrngs, kernel(is)) # parent(is) should be Colon or Integer + blockrngs = contact_kernel_ranges(kernel(is), parent(is), gs) # see greenfunction.jl + # note: blockrngs can be a lazy Iterator of unknown length + d = Diagonal(Returns(zero(C)).(blockrngs)) + fill_diagonal!(d, blockmat, blockrngs, kernel(is)) # parent(is) should be Colon or Integer + return d.diag end ## call @@ -276,16 +280,18 @@ function (d::DensityMatrixKPMSolver)(mu, kBT; params...) ω0, Δ = d.bandCH kBT´ = kBT / Δ mu´ = (mu - ω0) / Δ + result = copy(call!_output(d.gs)) + ρ = serialize(result) # obtain underlying array of result if kBT´ ≈ 0 ϕ = acos(mu´) - ρ = copy(first(d.momenta)) * (1-ϕ/π) + ρ .= first(d.momenta) .* (1-ϕ/π) for n in 1:length(d.momenta)-1 @. ρ += d.momenta[n+1] * 2 * (sin(π*n)-sin(n*ϕ))/(π*n) end else throw(argerror("KPM densitymatrix currently doesn't support finite temperatures. Open an issue if you need this feature.")) end - return ρ + return result end #endregion diff --git a/src/solvers/green/schur.jl b/src/solvers/green/schur.jl index 3c9be408..dfe4a007 100644 --- a/src/solvers/green/schur.jl +++ b/src/solvers/green/schur.jl @@ -660,14 +660,11 @@ decay_lengths(g::AbstractHamiltonian1D, µ = 0, args...; params...) = # Does not support non-Nothing contacts #region -struct DensityMatrixSchurSolver{T,A,G<:GreenFunctionSchurLead{T},R,O<:NamedTuple} - g::G # parent of GreenSlice +struct DensityMatrixSchurSolver{T,A,G<:GreenSlice{T},O<:NamedTuple} + gs::G # parent of GreenSlice orbaxes::A # axes of GreenSlice (orbrows, orbcols) - ρmat::R # it spans the full orbitalslices (axes) of GreenSlice. - # can be a Matrix or a Vector (if A::DiagIndices) hmat::Matrix{Complex{T}} # it spans just the unit cell, dense Bloch matrix psis::Matrix{Complex{T}} # it spans just the unit cell, Eigenstates - fmat::Matrix{Complex{T}} # it spans just the unit cell, f(hmat) opts::O # kwargs for quadgk end @@ -678,115 +675,49 @@ function densitymatrix(s::AppliedSchurGreenSolver, gs::GreenSlice; atol = 1e-7, argerror("Boundaries not implemented for DensityMatrixSchurSolver. Consider using the generic integration solver.") has_selfenergy(gs) && argerror("The Schur densitymatrix solver currently support only `nothing` contacts") g = parent(gs) - ρmat = similar_Array(gs) hmat = similar_Array(hamiltonian(g)) psis = similar(hmat) - fmat = similar(hmat) orbaxes = orbrows(gs), orbcols(gs) opts = (; atol, quadgk_opts...) - solver = DensityMatrixSchurSolver(g, orbaxes, ρmat, hmat, psis, fmat, opts) + solver = DensityMatrixSchurSolver(gs, orbaxes, hmat, psis, opts) return DensityMatrix(solver, gs) end -struct FermiEigenstates{T} - psis::Matrix{Complex{T}} - fpsis::Matrix{Complex{T}} -end - # API ## call # use computed Fermi points to integrate spectral function by segments +# returns an AbstractMatrix function (s::DensityMatrixSchurSolver)(µ, kBT; params...) - λs = propagating_eigvals(s.g, µ + 0im, 1e-2; params...) + g = parent(s.gs) + λs = propagating_eigvals(g, µ + 0im, 1e-2; params...) ϕs = @. real(-im*log(λs)) xs = sort!(ϕs ./ (2π)) pushfirst!(xs, -0.5) push!(xs, 0.5) xs = unique!(xs) - integrand(x) = fermi_h!(s, 2π * x, µ, inv(kBT); params...) - result = similar(s.ρmat) - fx! = (y, x) -> (y .= integrand(x)) + result = call!_output(s.gs) + integrand!(x) = fermi_h!(result, s, 2π * x, µ, inv(kBT); params...) + fx! = (y, x) -> (y .= integrand!(x)) quadgk!(fx!, result, xs...; s.opts...) return result end -function fermi_h!(s, ϕ, µ, β = 0; params...) - h = parent(s.g) +function fermi_h!(result, s, ϕ, µ, β = 0; params...) + h = hamiltonian(s.gs) + bs = blockstructure(h) # Similar to spectrum(h, ϕ; params...), but less work (no sort! or sanitization) copy!(s.hmat, call!(h, ϕ; params...)) # sparse to dense ϵs, psis = eigen!(s.hmat) # special-casing β = Inf with views turns out to be slower fs = (@. ϵs = fermi(ϵs - µ, β)) - fpsis = (s.psis .= fs .* psis') - ρmat = fill_rho_blocks!(s, psis, fpsis, ϕ, s.orbaxes...) - return ρmat -end - -# fillblock! consistent with Base.getindex(g::GreenSolution, args...) see greenfunction.jl -# uses views into fmat = f(ϵₙ) * ψₙ * ψₘ' to fill blocks in the ρmat = parent or OrbitalSliceMatrix -function fill_rho_blocks!(s::DensityMatrixSchurSolver, psis, fpsis, ϕ, i, j) - ρmat, fmat = s.ρmat, s.fmat - mul!(fmat, psis, fpsis) - oj = 0 # offsets - for cj in cellsdict_or_single_cellorbs(j) - oi = 0 - for ci in cellsdict_or_single_cellorbs(i) - ρview = view(ρmat, oi + 1:oi + norbitals(ci), oj + 1:oj + norbitals(cj)) - fview = view(fmat, orbindices(ci), orbindices(cj)) - copy!(ρview, fview) - phase = cis(only(cell(ci) - cell(cj)) * ϕ) - ρview .*= phase - oi += norbitals(ci) - end - oj += norbitals(cj) - end - return ρmat -end - -cellsdict_or_single_cellorbs(i::AnyOrbitalSlice) = cellsdict(i) -cellsdict_or_single_cellorbs(i::AnyCellOrbitals) = (i,) - -# fastpath for direct CellOrbitals indexing. Not sure it's worth it. -function fill_rho_blocks!(s::DensityMatrixSchurSolver, psis, fpsis, ϕ, i::AnyCellOrbitals, j::AnyCellOrbitals) - vpsis = view(psis, orbindices(i), :) - vfpsis = view(fpsis, :, orbindices(j)) - phase = cis(only(cell(i) - cell(j)) * ϕ) - mul!(s.ρmat, vpsis, vfpsis, phase, 0) - return s.ρmat -end - -# diagonal indexing -fill_rho_blocks!(s::DensityMatrixSchurSolver, psis, fpsis, ϕ, di::DiagIndices, ::DiagIndices) = - append_diagonal!(empty!(s.ρmat), FermiEigenstates(psis, fpsis), parent(di), kernel(di), s.g) - -# this driver gets called by append_diagonal(d, fe, o, kernel, g), see greenfunction.jl -function append_diagonal!(d, fe::FermiEigenstates, o::AnyCellOrbitals, kernel; post = identity) - for orbs in orbgroups(o) - v, v´ = view(fe.psis, orbs, :), view(fe.fpsis, :, orbs) - append_trace_or_diag!(d, kernel, v, v´, post) - end - return d -end - -function append_trace_or_diag!(d::AbstractVector{T}, ::Missing, v, v´, post) where {T} - for i in axes(v, 1) - val = zero(T) - for j in axes(v, 2) - val += v[i, j] * v´[j, i] - end - push!(d, post(val)) - end - return d + fpsis = (s.psis .= psis .* transpose(fs)) + ρcell = EigenProduct(bs, psis, fpsis, ϕ) + getindex!(result, ρcell, s.orbaxes...) + return result end -append_trace_or_diag!(d, kernel, v, v´, post) = push!(d, post(kernel_trace_multiply(kernel, v, v´))) - -kernel_trace_multiply(kernel::UniformScaling, v, v´) = kernel.λ * trace_prod(v, v´) -kernel_trace_multiply(kernel::Number, v, v´) = kernel * trace_prod(v, v´) -kernel_trace_multiply(kernel, v, v´) = trace_prod(kernel * v, v´) - #endregion #endregion top diff --git a/src/solvers/green/spectrum.jl b/src/solvers/green/spectrum.jl index ccb68d38..caad8336 100644 --- a/src/solvers/green/spectrum.jl +++ b/src/solvers/green/spectrum.jl @@ -74,13 +74,13 @@ end # specialized DensityMatrix method for GS.Spectrum #region -struct DensityMatrixSpectrumSolver{T,G<:GreenSlice,A,P,R} +struct DensityMatrixSpectrumSolver{T,G<:GreenSlice{T},A,P} gs::G - orbs::A + orbaxes::A es::Vector{Complex{T}} + fs::Vector{Complex{T}} psis::P fpsis::P - ρmat::R end ## Constructor @@ -90,37 +90,26 @@ function densitymatrix(s::AppliedSpectrumGreenSolver, gs::GreenSlice) # If rows/cols are contacts, we need their orbrows/orbcols (unlike for gs(ω; params...)) has_selfenergy(gs) && argerror("The Spectrum densitymatrix solver currently support only `nothing` contacts") es, psis = spectrum(s) - fpsis = copy(psis') - ρmat = similar_Array(gs) - orbs = orbrows(gs), orbcols(gs) - solver = DensityMatrixSpectrumSolver(gs, orbs, es, psis, fpsis, ρmat) + fpsis = copy(psis) + fs = copy(es) + orbaxes = orbrows(gs), orbcols(gs) + solver = DensityMatrixSpectrumSolver(gs, orbaxes, es, fs, psis, fpsis) return DensityMatrix(solver, gs) end ## call -function (s::DensityMatrixSpectrumSolver)(mu, kBT; params...) - psis, fpsis, es = s.psis, s.fpsis, s.es +function (s::DensityMatrixSpectrumSolver)(µ, kBT; params...) + bs = blockstructure(s.gs) + psis, fpsis, es, fs = s.psis, s.fpsis, s.es, s.fs β = inv(kBT) - @. fpsis = fermi(es - mu, β) * psis' - oi, oj = s.orbs - ρmat = fill_rho_blocks!(s, psis, fpsis, oi, oj) - return copy(ρmat) + (@. fs = fermi(es - µ, β)) + fpsis .= psis .* transpose(fs) + ρcell = EigenProduct(bs, psis, fpsis) + result = call!_output(s.gs) + getindex!(result, ρcell, s.orbaxes...) + return result end -function fill_rho_blocks!(s::DensityMatrixSpectrumSolver, psis, fpsis, i, j) - oi, oj = zerocellorbinds(i), zerocellorbinds(j) - vpsis = view(psis, oi, :) - vfpsis = view(fpsis, :, oj) - return mul!(s.ρmat, vpsis, vfpsis) -end - -zerocellorbinds(orb::AnyCellOrbitals) = orbindices(orb) -zerocellorbinds(orb::AnyOrbitalSlice) = orbindices(only(cellsdict(orb))) - -# this relies on the append_diagonal! method from the Schur densitymatrix solver, see schur.jl -fill_rho_blocks!(s::DensityMatrixSpectrumSolver, psis, fpsis, i::DiagIndices, ::DiagIndices) = - append_diagonal!(empty!(s.ρmat), FermiEigenstates(psis, fpsis), parent(i), kernel(i), parent(s.gs)) - #endregion #endregion diff --git a/src/solvers/selfenergy/generic.jl b/src/solvers/selfenergy/generic.jl index 6e01b4ef..12f4c2aa 100644 --- a/src/solvers/selfenergy/generic.jl +++ b/src/solvers/selfenergy/generic.jl @@ -47,7 +47,7 @@ function SelfEnergy(hparent::AbstractHamiltonian, gslice::GreenSlice, model::Abs # apply model to lat0 to get hcoupling interblockmodel = interblock(model, 1:nparent, nparent+1:ntotal) hcoupling = hamiltonian(lat0, interblockmodel; - orbitals = vcat(norbitals(hparent), norbitals(gslice))) + orbitals = vcat(norbitals(hparent, :), norbitals(gslice, :))) solver´ = SelfEnergyGenericSolver(gslice, hcoupling, nparent) return SelfEnergy(solver´, lsparent) end diff --git a/src/solvers/selfenergy/model.jl b/src/solvers/selfenergy/model.jl index ce3f0366..266923cb 100644 --- a/src/solvers/selfenergy/model.jl +++ b/src/solvers/selfenergy/model.jl @@ -18,7 +18,7 @@ function SelfEnergyModelSolver(h::AbstractHamiltonian, model::AbstractModel, orb # and fills siteinds::Vector{Int} with the site index for each lat0 site (i.e. for sites ordered by sublattices) lat0 = lattice0D(orbslice, siteinds) # this is a 0D ParametricHamiltonian to build the Σ(ω) as a view over flat(ph(; ...)) - ph = hamiltonian(lat0, modelω; orbitals = norbitals(h)) # WARNING: type-unstable orbs + ph = hamiltonian(lat0, modelω; orbitals = norbitals(h, :)) # WARNING: type-unstable orbs # translation from orbitals in lat0 to orbslice indices # i.e. orbital index on orbslice for each orbital in lat0 (this is just a reordering!) parentinds = reordered_site_orbitals(siteinds, orbslice) diff --git a/src/solvers/selfenergy/schur.jl b/src/solvers/selfenergy/schur.jl index 6591e277..13916295 100644 --- a/src/solvers/selfenergy/schur.jl +++ b/src/solvers/selfenergy/schur.jl @@ -208,7 +208,7 @@ function SelfEnergy(hparent::AbstractHamiltonian, glead::GreenFunctionSchurLead, # apply model to lat0 to get hcoupling interblockmodel = interblock(model, 1:nparent, nparent+1:ntotal) hcoupling = hamiltonian(lat0, interblockmodel; - orbitals = vcat(norbitals(hparent), norbitals(hlead))) + orbitals = vcat(norbitals(hparent, :), norbitals(hlead, :))) gslice = glead[cells = SA[xunit]] Σs = selfenergies(contacts(glead)) diff --git a/src/specialmatrices.jl b/src/specialmatrices.jl index 6e97208a..42abcb73 100644 --- a/src/specialmatrices.jl +++ b/src/specialmatrices.jl @@ -61,7 +61,7 @@ function flat(b::OrbitalBlockStructure{B}, unflat::SparseMatrixCSC{B´}) where { return sparse(builder, n, n) end -function unflat(b::OrbitalBlockStructure{B}, flat::SparseMatrixCSC{<:Number}) where {N,B<:MatrixElementNonscalarType{<:Any,N}} +Base.@propagate_inbounds function unflat(b::OrbitalBlockStructure{B}, flat::SparseMatrixCSC{<:Number}) where {N,B<:MatrixElementNonscalarType{<:Any,N}} @boundscheck(checkblocks(b, flat)) nnzguess = nnz(flat) ÷ (N * N) ncols = unflatsize(b) @@ -131,7 +131,7 @@ Base.getindex(b::HybridSparseMatrix{<:Any,<:SMatrixView}, i::Integer, j::Integer Base.getindex(b::HybridSparseMatrix, i::Integer, j::Integer) = unflat(b)[i, j] # only allowed for elements that are already stored -function Base.setindex!(b::HybridSparseMatrix{<:Any,B}, val::AbstractVecOrMat, i::Integer, j::Integer) where {B<:SMatrixView} +Base.@propagate_inbounds function Base.setindex!(b::HybridSparseMatrix{<:Any,B}, val::AbstractVecOrMat, i::Integer, j::Integer) where {B<:SMatrixView} @boundscheck(checkstored(unflat(b), i, j)) val´ = mask_block(B, val, blocksize(blockstructure(b), i, j)) unflat(b)[i, j] = val´ @@ -139,7 +139,7 @@ function Base.setindex!(b::HybridSparseMatrix{<:Any,B}, val::AbstractVecOrMat, i return val´ end -function Base.setindex!(b::HybridSparseMatrix, val::AbstractVecOrMat, i::Integer, j::Integer) +Base.@propagate_inbounds function Base.setindex!(b::HybridSparseMatrix, val::AbstractVecOrMat, i::Integer, j::Integer) @boundscheck(checkstored(unflat(b), i, j)) unflat(b)[i, j] = val needs_flat_sync!(b) @@ -512,6 +512,13 @@ Base.@propagate_inbounds Base.getindex(a::OrbitalSliceArray, I::Vararg{Int, N}) Base.@propagate_inbounds Base.setindex!(a::OrbitalSliceArray, v, i::Int) = setindex!(parent(a), v, i) Base.@propagate_inbounds Base.setindex!(a::OrbitalSliceArray, v, I::Vararg{Int, N}) where {N} = setindex!(parent(a), v, I...) +Base.:*(x::Number, a::OrbitalSliceArray) = OrbitalSliceArray(parent(a) * x, orbaxes(a)) +Base.:*(a::OrbitalSliceArray, x::Number) = x * a + +Base.:+(a::OrbitalSliceArray, b::OrbitalSliceArray) = parent(a) + parent(b) + +Base.:-(a::OrbitalSliceArray, b::OrbitalSliceArray) = parent(a) - parent(b) + # Additional indexing over sites Base.getindex(a::OrbitalSliceMatrix; sites...) = getindex(a, siteselector(; sites...)) Base.getindex(a::OrbitalSliceMatrix, i::NamedTuple, j::NamedTuple = i) = @@ -531,7 +538,7 @@ function Base.getindex(a::OrbitalSliceMatrix, i::SiteSelector, j::SiteSelector = return OrbitalSliceMatrix(m, (rowslice´, colslice´)) end -# CellSites: return an unwrapped Matrix or a view thereof (non-allocating) +# CellSites: return an unwrapped Matrix Base.getindex(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) = copy(view(a, i, j)) @@ -546,21 +553,41 @@ function Base.view(a::OrbitalSliceMatrix, i::AnyCellSites, j::AnyCellSites = i) return view(parent(a), rows, cols) end +# Minimal getindex for OrbitalSliceVector +Base.getindex(a::OrbitalSliceVector, i::AnyCellSites) = copy(view(a, i)) +Base.getindex(a::OrbitalSliceVector, i::C) where {B,C<:CellSitePos{<:Any,<:Any,<:Any,B}} = + sanitize_block(B, view(a, i)) + +function Base.view(a::OrbitalSliceVector, i::AnyCellSites) + rowslice = only(orbaxes(a)) + i´ = apply(i, lattice(rowslice)) + rows = indexcollection(rowslice, i´) + return view(parent(a), rows) +end + +LinearAlgebra.diag(a::OrbitalSliceMatrix) = + OrbitalSliceVector(diag(parent(a)), (first(orbaxes(a)),)) + +LinearAlgebra.norm(a::OrbitalSliceMatrix) = norm(parent(a)) + ## broadcasting # following the manual: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array +# disabled for now, since it seems to cause issues with OrbitalSliceMatrix{<:SparseMatrixCSC} -Broadcast.BroadcastStyle(::Type{<:OrbitalSliceArray}) = Broadcast.ArrayStyle{OrbitalSliceArray}() +# Broadcast.BroadcastStyle(::Type{O}) where {O<:OrbitalSliceArray} = Broadcast.ArrayStyle{O}() -Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{OrbitalSliceArray}}, ::Type{ElType}) where {ElType} = - OrbitalSliceArray(similar(Array{ElType}, axes(bc)), orbaxes(find_osa(bc))) +# Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{O}}, ::Type{ElType}) where {ElType,O<:OrbitalSliceMatrix{<:Any,<:Array}} = +# OrbitalSliceArray(similar(Array{ElType}, axes(bc)), orbaxes(find_osa(bc))) +# Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{O}}, ::Type{ElType}) where {ElType,O<:OrbitalSliceMatrix{<:Any,<:SparseMatrixCSC}} = +# OrbitalSliceArray(similar(parent(find_osa(bc)), ElType), orbaxes(find_osa(bc))) -find_osa(bc::Base.Broadcast.Broadcasted) = find_osa(bc.args) -find_osa(args::Tuple) = find_osa(find_osa(args[1]), Base.tail(args)) -find_osa(x) = x -find_osa(::Tuple{}) = nothing -find_osa(a::OrbitalSliceArray, rest) = a -find_osa(::Any, rest) = find_osa(rest) +# find_osa(bc::Base.Broadcast.Broadcasted) = find_osa(bc.args) +# find_osa(args::Tuple) = find_osa(find_osa(args[1]), Base.tail(args)) +# find_osa(x) = x +# find_osa(::Tuple{}) = nothing +# find_osa(a::OrbitalSliceArray, rest) = a +# find_osa(::Any, rest) = find_osa(rest) # taken from https://github.com/JuliaArrays/OffsetArrays.jl/blob/756e839563c88faa4ebe4ff971286747863aaff0/src/OffsetArrays.jl#L469 @@ -568,29 +595,97 @@ Base.dataids(A::OrbitalSliceArray) = Base.dataids(parent(A)) Broadcast.broadcast_unalias(dest::OrbitalSliceArray, src::OrbitalSliceArray) = parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src) -## conversion +# ## conversion: If i contains an OrbitalSliceGrouped, wrap it in an OrbitalSliceArray. Otherwise, no-op -maybe_OrbitalSliceArray(i) = x -> maybe_OrbitalSliceArray(x, i) +# maybe_OrbitalSliceArray(i) = x -> maybe_OrbitalSliceArray(x, i) -maybe_OrbitalSliceArray(x::AbstractVector, i) = maybe_OrbitalSliceVector(x, i) -maybe_OrbitalSliceArray(x::AbstractVector, (i,_)::Tuple) = maybe_OrbitalSliceVector(x, i) -maybe_OrbitalSliceArray(x::AbstractMatrix, i) = maybe_OrbitalSliceMatrix(x, i) +# maybe_OrbitalSliceArray(x::AbstractMatrix, i) = maybe_OrbitalSliceMatrix(x, i) -maybe_OrbitalSliceVector(x, i::DiagIndices{Missing,<:OrbitalSliceGrouped}) = - maybe_OrbitalSliceVector(x, parent(i)) -maybe_OrbitalSliceVector(x, i::DiagIndices{<:Any,<:OrbitalSliceGrouped}) = - maybe_OrbitalSliceVector(x, scalarize(parent(i))) -maybe_OrbitalSliceVector(x, i::OrbitalSliceGrouped) = OrbitalSliceVector(x, (i,)) +# maybe_OrbitalSliceMatrix(x::Diagonal, i::OrbitalSliceGrouped) = +# maybe_OrbitalSliceMatrix(x, (i, i)) +# maybe_OrbitalSliceMatrix(x, (i, j)::Tuple{OrbitalSliceGrouped,OrbitalSliceGrouped}) = +# OrbitalSliceMatrix(x, (i, j)) +# maybe_OrbitalSliceMatrix(x, i) = x -# fallback -maybe_OrbitalSliceVector(x, i) = x +## currently unsused +# maybe_OrbitalSliceArray(x::AbstractVector, i) = maybe_OrbitalSliceVector(x, i) +# maybe_OrbitalSliceVector(x, i::DiagIndices{Missing,<:OrbitalSliceGrouped}) = +# maybe_OrbitalSliceVector(x, parent(i)) +# maybe_OrbitalSliceVector(x, i::DiagIndices{<:Any,<:OrbitalSliceGrouped}) = +# maybe_OrbitalSliceVector(x, scalarize(parent(i))) +# maybe_OrbitalSliceVector(x, i::OrbitalSliceGrouped) = OrbitalSliceVector(x, (i,)) +# maybe_OrbitalSliceVector(x, i) = x + + +#endregion + +############################################################################################ +## EigenProduct +# A matrix P = X * U * V' that is stored in terms of matrices U and V without computing +# the product. x is a scalar. Inspired by LowRankMatrices.jl +#region + +struct EigenProduct{T,M<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} <: AbstractMatrix{T} + blockstruct::O + U::M + V::M + phi::T + x::Complex{T} + function EigenProduct{T,M,O}(blockstruct::O, U::M, V::M, phi::T, x::Complex{T}) where {T,M<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} + axes(U, 2) != axes(V, 2) && argerror("U and V must have identical column axis") + return new(blockstruct, U, V, phi, x) + end +end -maybe_OrbitalSliceMatrix(x, i::OrbitalSliceGrouped) = - maybe_OrbitalSliceMatrix(x, (i, i)) -maybe_OrbitalSliceMatrix(x, (i, j)::Tuple{OrbitalSliceGrouped,OrbitalSliceGrouped}) = - OrbitalSliceMatrix(x, (i, j)) +#region ## Constructors ## -# fallback -maybe_OrbitalSliceMatrix(x, i) = x +EigenProduct(blockstruct::O, U::M, V::M, phi = zero(T), x = one(Complex{T})) where {T,M<:AbstractMatrix{Complex{T}},O<:OrbitalBlockStructure} = + EigenProduct{T,M,O}(blockstruct, U, V, T(phi), Complex{T}(x)) #endregion + +#region ## API ## + +blockstructure(P::EigenProduct) = P.blockstruct + +Base.@propagate_inbounds function Base.getindex(P::EigenProduct, i::Integer, j::Integer) + ret = zero(eltype(P)) + @boundscheck checkbounds(P.U, i, :) + @boundscheck checkbounds(P.V, j, :) + @inbounds for k in axes(P.U, 2) + ret = muladd(P.U[i, k], conj(P.V[j, k]), ret) + end + return P.x * ret +end + +Base.@propagate_inbounds function Base.getindex(P::EigenProduct, i::Integer) + ij = Tuple(CartesianIndices(axes(P))[i]) + return P[ij...] +end + +function Base.getindex(P::EigenProduct{T}, ci::AnyCellOrbitals, cj::AnyCellOrbitals) where {T} + phase = isempty(cell(ci)) ? one(T) : cis(only(cell(ci) - cell(cj)) * P.phi) + return view(phase * P, orbindices(ci), orbindices(cj)) +end + +Base.view(P::EigenProduct, is, js) = + EigenProduct(P.blockstruct, view(P.U, is, :), view(P.V, js, :), P.phi, P.x) + +# AbstractArray interface +Base.axes(P::EigenProduct) = axes(P.U, 1), axes(P.V, 1) +Base.size(P::EigenProduct) = size(P.U, 1), size(P.V, 1) +# Base.iterate(a::EigenProduct, i...) = iterate(parent(a), i...) +Base.length(P::EigenProduct) = prod(size(P)) +Base.IndexStyle(::Type{T}) where {M,T<:EigenProduct{<:Any,M}} = IndexStyle(M) +Base.similar(P::EigenProduct) = EigenProduct(P.blockstruct, similar(P.U), similar(P.V), P.phi, P.x) +Base.similar(P::EigenProduct, t::Type) = EigenProduct(P.blockstruct, similar(P.U, t), similar(P.V, t), P.phi, P.x) +Base.copy(P::EigenProduct) = EigenProduct(P.blockstruct, copy(P.U), copy(P.V), P.phi, P.x) +Base.eltype(::EigenProduct{T}) where {T} = Complex{T} + +Base.:*(x::Number, P::EigenProduct) = EigenProduct(P.blockstruct, P.U, P.V, P.phi, P.x * x) +Base.:*(P::EigenProduct, x::Number) = x*P + +LinearAlgebra.tr(P::EigenProduct) = sum(xy -> first(xy) * conj(last(xy)), zip(P.U, P.V)) * P.x + +#endregion +#endregion diff --git a/src/supercell.jl b/src/supercell.jl index a3cc2800..3b9fd12b 100644 --- a/src/supercell.jl +++ b/src/supercell.jl @@ -175,7 +175,7 @@ function supercell(h::Hamiltonian, data::SupercellData; mincoordination = 0) indexlist, offset = supercell_indexlist!(data, h, mincoordination) lat´ = lattice(data) B = blocktype(h) - bs´ = OrbitalBlockStructure{B}(norbitals(h), sublatlengths(lat´)) + bs´ = OrbitalBlockStructure{B}(norbitals(h, :), sublatlengths(lat´)) builder = CSCBuilder(lat´, bs´) har´ = supercell_harmonics(h, data, builder, indexlist, offset) return Hamiltonian(lat´, bs´, har´) diff --git a/src/tools.jl b/src/tools.jl index 4e73407b..2554be67 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -71,7 +71,7 @@ chopsmall(xs, atol) = chopsmall.(xs, Ref(atol)) chopsmall(xs) = chopsmall.(xs) # Flattens matrix of Matrix{<:Number} into a matrix of Number's -function mortar(ms::AbstractMatrix{M}) where {C<:Number,M<:AbstractMatrix{C}} +function mortar(ms::AbstractMatrix{M}) where {C<:Number,M<:Matrix{C}} isempty(ms) && return convert(Matrix{C}, ms) mrows = size.(ms, 1) mcols = size.(ms, 2) @@ -89,7 +89,35 @@ function mortar(ms::AbstractMatrix{M}) where {C<:Number,M<:AbstractMatrix{C}} end return mat end -# faspath for generators or other collections + +#equivalent with SparseMatrixCSC +function mortar(ms::AbstractMatrix{M}) where {C<:Number,M<:SparseMatrixCSC{C}} + isempty(ms) && return convert(SparseMatrixCSC{C,Int}, ms) + mrows = size.(ms, 1) + mcols = size.(ms, 2) + allequal(eachrow(mcols)) && allequal(eachcol(mrows)) || + internalerror("mortar: inconsistent rows or columns") + totalrows = sum(i -> size(ms[i, 1], 1), axes(ms, 1)) + totalcols = sum(i -> size(ms[1, i], 2), axes(ms, 2)) + totalnnz = sum(nnz, ms) + b = CSC{C}(totalcols, totalnnz) + for mj in axes(ms, 2), col in axes(ms[1, mj], 2) + rowoffset = 0 + for mi in axes(ms, 1) + m = ms[mi, mj] + for ptr in nzrange(m, col) + row = rowoffset + rowvals(m)[ptr] + val = nonzeros(m)[ptr] + pushtocolumn!(b, row, val) + end + rowoffset += size(m, 1) + end + finalizecolumn!(b) + end + return sparse(b, totalrows, totalcols) +end + +# # faspath for generators or other collections mortar(ms) = length(ms) == 1 ? only(ms) : mortar(collect(ms)) # equivalent to mat = I[:, cols]. Useful for Green function source @@ -109,8 +137,25 @@ lengths_to_offsets(v) = prepend!(cumsum(v), 0) lengths_to_offsets(f::Function, v) = prepend!(accumulate((i,j) -> i + f(j), v; init = 0), 0) -# fast tr(A*B) without doing the A*B product -trace_prod(A, B) = sum(splat(*), zip(A, transpose(B))) +# fast tr(A*B) +trace_prod(A, B) = (check_sizes(A,B); unsafe_trace_prod(A, B)) + +unsafe_trace_prod(A::AbstractMatrix, B::AbstractMatrix) = sum(splat(*), zip(transpose(A), B)) +unsafe_trace_prod(A::Number, B::AbstractMatrix) = A*tr(B) +unsafe_trace_prod(A::AbstractMatrix, B::Number) = unsafe_trace_prod(B, A) +unsafe_trace_prod(A::UniformScaling, B::AbstractMatrix) = A.λ*tr(B) +unsafe_trace_prod(A::AbstractMatrix, B::UniformScaling) = unsafe_trace_prod(B, A) +unsafe_trace_prod(A::Diagonal, B::Diagonal) = sum(i -> A[i] * B[i], axes(A,1)) +unsafe_trace_prod(A::Diagonal, B::AbstractMatrix) = sum(i -> A[i] * B[i, i], axes(A,1)) +unsafe_trace_prod(A::AbstractMatrix, B::Diagonal) = unsafe_trace_prod(B, A) +unsafe_trace_prod(A::Diagonal, B::Number) = only(A) * B +unsafe_trace_prod(A::Number, B::Diagonal) = unsafe_trace_prod(B, A) +unsafe_trace_prod(A::Union{SMatrix,UniformScaling,Number}, B::Union{SMatrix,UniformScaling,Number}) = + tr(A*B) + +check_sizes(A::AbstractMatrix,B::AbstractMatrix) = size(A,2) == size(B,1) || + throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))")) +check_sizes(_, _) = nothing # Taken from Base julia, now deprecated there function permute!!(a, p::AbstractVector{<:Integer}) @@ -156,7 +201,7 @@ function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer}) a end -function swapcols!(a::AbstractMatrix, i, j) +Base.@propagate_inbounds function swapcols!(a::AbstractMatrix, i, j) i == j && return cols = axes(a,2) @boundscheck i in cols || throw(BoundsError(a, (:,i))) diff --git a/src/types.jl b/src/types.jl index 1d341c77..ee0e91ef 100644 --- a/src/types.jl +++ b/src/types.jl @@ -428,9 +428,12 @@ CellOrbitalsGrouped(cell, inds, groups::Dictionary) = # B <: MatrixElementType{T} CellSitePos(cell, ind, r, B) = CellIndices(sanitize_SVector(Int, cell), ind, SiteLikePos(r, B)) +# changing only cell +CellIndices(cell, c::CellIndices) = CellIndices(cell, c.inds, c.type) + # LatticeSlice from an AbstractVector of CellIndices LatticeSlice(lat::Lattice, cs::AbstractVector{<:CellIndices}) = - LatticeSlice(lat, cellinds_to_dict(cs)) + LatticeSlice(lat, cellinds_to_dict(apply.(cs, Ref(lat)))) # CellIndices to Dictionary(cell=>cellind) cellinds_to_dict(cs::AbstractVector{C}) where {L,C<:CellIndices{L}} = @@ -442,6 +445,11 @@ cellinds_to_dict(cs::CellIndices{L}) where {L} = cellinds_to_dict(SVector(cs)) cellinds_to_dict(cs::AbstractVector{C}) where {L,C<:CellSite{L}} = cellinds_to_dict(CellSites{L,Vector{Int}}.(cs)) +# scalarize only if kernel is not missing +maybe_scalarize(s) = s +maybe_scalarize(s, kernel::Missing) = s +maybe_scalarize(s, kernel) = scalarize(s) + # make OrbitalSliceGrouped scalar (i.e. groups are siteindex => siteindex:siteindex) scalarize(s::OrbitalSliceGrouped) = OrbitalSliceGrouped(s.lat, scalarize(s.cellsdict)) scalarize(d::CellOrbitalsGroupedDict) = scalarize.(d) @@ -535,6 +543,9 @@ ind(s::CellSitePos) = s.inds sites(l::LatticeSlice) = (site(l.lat, i, cell(subcell)) for subcell in cellsdict(l) for i in siteindices(subcell)) +sites(l::LatticeSlice, blocktype::Type) = + (CellIndices(cell(subcell), i, SiteLikePos(site(l.lat, i), blocktype)) + for subcell in cellsdict(l) for i in siteindices(subcell)) # single-site (CellSite) iterator cellsites(cs::Union{SiteSlice,OrbitalSliceGrouped}) = cellsites(cs.cellsdict) @@ -557,10 +568,18 @@ orbranges(cs::OrbitalSliceGrouped) = orbranges(cellsdict(cs)) orbranges(cs::Union{CellOrbitalsGrouped, CellOrbitalsGroupedDict}) = Iterators.accumulate( (rng, rng´) -> maximum(rng)+1:maximum(rng)+length(rng´), orbgroups(cs), init = 0:0) +# range of orbital indices of dn cellorbs +function orbrange(cs::AnyOrbitalSlice, dn::SVector) + offset = get(offsets(cs), dn, 0) + rng = offset+1:offset+norbitals(cs, dn) + return rng +end + Base.isempty(s::LatticeSlice) = isempty(s.cellsdict) Base.isempty(s::CellIndices) = isempty(s.inds) -Base.length(l::LatticeSlice) = nsites(l) +Base.length(l::SiteSlice) = nsites(l) +Base.length(l::AnyOrbitalSlice) = norbitals(l) Base.length(c::CellIndices) = length(c.inds) Base.parent(ls::LatticeSlice) = ls.lat @@ -939,7 +958,7 @@ function sublatorbrange(b::OrbitalBlockStructure, sind::Integer) end # Basic relation: iflat - 1 == (iunflat - soffset - 1) * b + soffset´ -function flatrange(b::OrbitalBlockStructure{<:SMatrixView}, iunflat::Integer) +Base.@propagate_inbounds function flatrange(b::OrbitalBlockStructure{<:SMatrixView}, iunflat::Integer) checkinrange(iunflat, b) soffset = 0 soffset´ = 0 @@ -971,13 +990,13 @@ function flatrange(o::OrbitalBlockStructure, runflat::AbstractUnitRange) return orng end -checkinrange(siteind::Integer, b::OrbitalBlockStructure) = +Base.@propagate_inbounds checkinrange(siteind::Integer, b::OrbitalBlockStructure) = @boundscheck(1 <= siteind <= flatsize(b) || argerror("Requested site $siteind out of range [1, $(flatsize(b))]")) flatindex(b::OrbitalBlockStructure, i) = first(flatrange(b, i)) -function unflatindex_and_blocksize(o::OrbitalBlockStructure{<:SMatrixView}, iflat::Integer) +Base.@propagate_inbounds function unflatindex_and_blocksize(o::OrbitalBlockStructure{<:SMatrixView}, iflat::Integer) soffset = 0 soffset´ = 0 @boundscheck(iflat < 0 && blockbounds_error()) @@ -1004,6 +1023,11 @@ Base.copy(b::OrbitalBlockStructure{B}) where {B} = Base.:(==)(b::OrbitalBlockStructure{B}, b´::OrbitalBlockStructure{B}) where {B} = b.blocksizes == b´.blocksizes && b.subsizes == b´.subsizes + +# make OrbitalBlockStructure scalar +scalarize(b::OrbitalBlockStructure) = + OrbitalBlockStructure{blockeltype(b)}(Returns(1).(b.blocksizes), b.subsizes) + #endregion #endregion @@ -1322,6 +1346,8 @@ flat_unsafe(h::Harmonic) = flat_unsafe(h.h) unflat_unsafe(h::Harmonic) = unflat_unsafe(h.h) +blocktype(::Harmonic{<:Any,<:Any,B}) where {B} = B + Base.size(h::Harmonic, i...) = size(matrix(h), i...) Base.isless(h::Harmonic, h´::Harmonic) = sum(abs2, dcell(h)) < sum(abs2, dcell(h´)) @@ -1491,7 +1517,8 @@ latdim(h::AbstractHamiltonian) = latdim(lattice(h)) bravais_matrix(h::AbstractHamiltonian) = bravais_matrix(lattice(h)) -norbitals(h::AbstractHamiltonian) = blocksizes(blockstructure(h)) +norbitals(h::AbstractHamiltonian, ::Colon) = blocksizes(blockstructure(h)) +norbitals(h::AbstractHamiltonian) = flatsize(h) blockeltype(h::AbstractHamiltonian) = blockeltype(blockstructure(h)) @@ -1532,7 +1559,7 @@ blocktype(::Type{<:AbstractHamiltonian{<:Any,<:Any,<:Any,B}}) where {B} = B # lat must be the result of combining the lattices of h, hs... function blockstructure(lat::Lattice{T}, h::AbstractHamiltonian{T}, hs::AbstractHamiltonian{T}...) where {T} B = blocktype(h, hs...) - orbitals = sanitize_orbitals(vcat(norbitals.((h, hs...))...)) + orbitals = sanitize_orbitals(vcat(norbitals.((h, hs...), Ref(:))...)) subsizes = sublatlengths(lat) return OrbitalBlockStructure{B}(orbitals, subsizes) end @@ -2174,53 +2201,102 @@ struct GreenSolution{T,E,L,S<:GreenSlicer,G<:GreenFunction{T,E,L},Σs} contactorbs::ContactOrbitals{L} end -struct DiagIndices{K,V} # represents Green indices to return only diagonal elements +# represents indices of sparse Green elements +struct SparseIndices{V,K} inds::V kernel::K end +const DiagIndices{V<:Union{SiteSelector,LatticeSlice,CellIndices,Integer,Colon},K} = SparseIndices{V,K} +const PairIndices{V<:HopSelector,K} = SparseIndices{V,K} + +const OrbitalDiagIndices{K} = SparseIndices{<:Union{AnyCellOrbitals,AnyOrbitalSlice},K} +const OrbitalPairIndices{K} = SparseIndices{<:Hamiltonian,K} + struct GreenIndices{I,R} - inds::I - orbinds::R # orbinds = sites_to_orbs(inds) + inds::I # Can be Colon, Integer, SiteSelector, LatticeSlice, CellIndices, SparseIndices... + orbinds::R # orbinds = sites_to_orbs(inds). Can be AnyOrbitalSlice, AnyCellOrbitals or SparseIndices thereof end # Obtained with gs = g[; siteselection...] # Alows call!(gs, ω; params...) or gs(ω; params...) # required to do e.g. h |> attach(g´[sites´], couplingmodel; sites...) -struct GreenSlice{T,E,L,G<:GreenFunction{T,E,L},R<:GreenIndices,C<:GreenIndices} +struct GreenSlice{T,E,L,G<:GreenFunction{T,E,L},R<:GreenIndices,C<:GreenIndices,O<:AbstractArray} parent::G rows::R cols::C + output::O # this will hold the result of call!(gs, ω; params...) end #region ## Constuctors ## -function GreenSlice(parent, rows, cols) +GreenSlice(parent::GreenFunction{T}, rows::GreenIndices, cols::GreenIndices, ::Type{C} = Complex{T}) where {C,T} = + GreenSlice(parent, rows, cols, similar_slice(C, orbindices(rows), orbindices(cols))) + +function GreenSlice(parent, rows, cols, type...) if rows isa DiagIndices || cols isa DiagIndices rows === cols || argerror("Diagonal indices should be identical for rows and columns") end rows´ = greenindices(rows, parent) cols´ = cols === rows ? rows´ : greenindices(cols, parent) - return GreenSlice(parent, rows´, cols´) + return GreenSlice(parent, rows´, cols´, type...) end # see sites_to_orbs in slices.jl greenindices(inds, g) = GreenIndices(inds, sites_to_orbs(inds, g)) greenindices(g::GreenSlice) = g.rows, g.cols +# preallocated AbstractArray for the output of call!(::GreenSlice, ω; params...) +# - OrbitalSliceMatrix for i,j both OrbitalSliceGrouped +# - OrbitalSliceMatrix for i,j both SparseIndices of OrbitalSliceGrouped +# - Matrix otherwise +similar_slice(::Type{C}, i::OrbitalSliceGrouped, j::OrbitalSliceGrouped) where {C} = + OrbitalSliceMatrix{C}(undef, i, j) + +similar_slice(::Type{C}, i::S, j::S) where {C,S<:SparseIndices{<:Union{AnyOrbitalSlice,Hamiltonian}}} = + OrbitalSliceMatrix{C}(undef, maybe_scalarize(i)) + +# If any of rows, cols are not both AnyOrbitalSlice, the return type is a Matrix or Diagonal +# its size is usually norbitals, but may be nsites if orbindices are a SparseIndices with kernel +similar_slice(::Type{C}, i::S, j::S) where {C,S<:SparseIndices{<:AnyCellOrbitals}} = + Diagonal(Vector{C}(undef, length(maybe_scalarize(i)))) + +similar_slice(::Type{C}, i, j) where {C} = + Matrix{C}(undef, length(i), length(j)) + +maybe_scalarize(i::SparseIndices{<:Union{AnyOrbitalSlice,AnyCellOrbitals}}) = + SparseIndices(maybe_scalarize(i.inds, i.kernel), i.kernel) +maybe_scalarize(i::SparseIndices) = i + #endregion #region ## API ## -diagonal(inds; kernel = missing) = DiagIndices(inds, kernel) # exported -diagonal(; kernel = missing, kw...) = DiagIndices(siteselector(; kw...), kernel) # exported +diagonal(inds; kernel = missing) = SparseIndices(inds, kernel) # exported +diagonal(inds::NamedTuple; kernel = missing) = SparseIndices(siteselector(; inds...), kernel) # exported +diagonal(; kernel = missing, kw...) = SparseIndices(siteselector(; kw...), kernel) # exported + +# sitepairs API only accepts HopSelector +sitepairs(h::HopSelector; kernel = missing) = SparseIndices(h, kernel) # exported +sitepairs(; kernel = missing, kw...) = SparseIndices(hopselector(; kw...), kernel) # exported -Base.parent(i::DiagIndices) = i.inds +norbitals(i::SparseIndices) = norbitals(parent(i)) -kernel(i::DiagIndices) = i.kernel +Base.parent(i::SparseIndices) = i.inds -norbitals_or_sites(i::DiagIndices{Missing}) = norbitals(i.inds) -norbitals_or_sites(i::DiagIndices) = nsites(i.inds) +Base.length(i::SparseIndices{<:Any,Missing}) = norbitals(parent(i)) +Base.length(i::SparseIndices) = nsites(parent(i)) + +Base.size(i::OrbitalDiagIndices) = length(i) .* (1, 1) +Base.size(i::OrbitalPairIndices) = length(i) .* (length(harmonics(parent(i))), 1) + +kernel(i::SparseIndices) = i.kernel + +hamiltonian(i::SparseIndices{<:Hamiltonian}) = i.inds + +Base.parent(i::GreenIndices) = i.inds + +orbindices(i::GreenIndices) = i.orbinds # returns the Hamiltonian field hamiltonian(g::Union{GreenFunction,GreenSolution,GreenSlice}) = hamiltonian(g.parent) @@ -2266,8 +2342,8 @@ blockstructure(g::GreenFunction) = blockstructure(hamiltonian(g)) blockstructure(g::GreenSolution) = blockstructure(hamiltonian(g)) blockstructure(g::GreenSlice) = blockstructure(parent(g)) -norbitals(g::GreenFunction) = norbitals(g.parent) -norbitals(g::GreenSlice) = norbitals(g.parent.parent) +norbitals(g::GreenFunction, args...) = norbitals(g.parent, args...) +norbitals(g::GreenSlice, args...) = norbitals(g.parent.parent, args...) contactinds(g::GreenFunction, i...) = contactinds(contacts(g), i...) contactinds(g::Union{GreenSolution,GreenSlice}, i...) = contactinds(contactorbitals(g), i...) @@ -2284,13 +2360,9 @@ orbcols(g::GreenSlice) = g.cols.orbinds Base.axes(g::GreenSlice) = (orbrows(g), orbcols(g)) -# ifelse(rows && cols are contacts, (rows, cols), (orbrows, orbcols)) -# I.e: if rows, cols are contact indices retrieve them instead of orbslices. -orbinds_or_contactinds(g) = orbinds_or_contactinds(rows(g), cols(g), orbrows(g), orbcols(g)) -orbinds_or_contactinds( - r::Union{Colon,Integer,DiagIndices{<:Any,Colon},DiagIndices{<:Any,<:Integer}}, - c::Union{Colon,Integer,DiagIndices{<:Any,Colon},DiagIndices{<:Any,<:Integer}}, _, _) = (r, c) -orbinds_or_contactinds(_, _, or, oc) = (or, oc) +call!_output(g::GreenSlice) = g.output + +replace_output!(g::GreenSlice, output) = GreenSlice(g.parent, g.rows, g.cols, output) Base.parent(g::GreenFunction) = g.parent Base.parent(g::GreenSolution) = g.parent @@ -2302,10 +2374,10 @@ Base.size(g::GreenSolution, i...) = size(g.parent, i...) flatsize(g::GreenFunction, i...) = flatsize(g.parent, i...) flatsize(g::GreenSolution, i...) = flatsize(g.parent, i...) -similar_Array(gs::GreenSlice{T}) where {T} = - matrix_or_vector(Complex{T}, orbrows(gs), orbcols(gs)) -matrix_or_vector(C, r, c) = Matrix{C}(undef, norbitals(r), norbitals(c)) -matrix_or_vector(C, r::DiagIndices, ::DiagIndices) = Vector{C}(undef, norbitals_or_sites(r)) +# similar_Array(gs::GreenSlice{T}) where {T} = +# matrix_or_vector(Complex{T}, orbrows(gs), orbcols(gs)) +# matrix_or_vector(C, r, c) = Matrix{C}(undef, norbitals(r), norbitals(c)) +# matrix_or_vector(C, r::DiagIndices, ::DiagIndices) = Vector{C}(undef, norbitals_or_sites(r)) boundaries(g::GreenFunction) = boundaries(solver(g)) # fallback (for solvers without boundaries, or for OpenHamiltonian) @@ -2318,7 +2390,7 @@ copy_lattice(g::GreenFunction) = GreenFunction(copy_lattice(g.parent), g.solver, copy_lattice(g::GreenSolution) = GreenSolution( copy_lattice(g.parent), g.slicer, g.contactΣs, g.contactbs) copy_lattice(g::GreenSlice) = GreenSlice( - copy_lattice(g.parent), g.rows, g.cols) + copy_lattice(g.parent), g.rows, g.cols, g.output) function minimal_callsafe_copy(g::GreenFunction) parent´ = minimal_callsafe_copy(g.parent) @@ -2337,7 +2409,7 @@ function minimal_callsafe_copy(g::GreenSolution) end minimal_callsafe_copy(g::GreenSlice) = - GreenSlice(minimal_callsafe_copy(g.parent), g.rows, g.cols) + GreenSlice(minimal_callsafe_copy(g.parent), g.rows, g.cols, copy(g.output)) Base.:(==)(g::GreenFunction, g´::GreenFunction) = function_not_defined("==") Base.:(==)(g::GreenSolution, g´::GreenSolution) = function_not_defined("==") @@ -2437,18 +2509,23 @@ const OrbitalSliceMatrix{C,M,A} = OrbitalSliceArray{C,2,M,A} OrbitalSliceVector(v::AbstractVector, axes) = OrbitalSliceArray(v, axes) OrbitalSliceMatrix(m::AbstractMatrix, axes) = OrbitalSliceArray(m, axes) -orbaxes(a::OrbitalSliceArray) = a.orbaxes +OrbitalSliceMatrix{C}(::UndefInitializer, oi, oj) where {C} = + OrbitalSliceArray(Matrix{C}(undef, length(oi), length(oj)), (oi, oj)) -Base.parent(a::OrbitalSliceArray) = a.parent +function OrbitalSliceMatrix{C}(::UndefInitializer, oi::SparseIndices{<:OrbitalSliceGrouped}) where {C} + i = j = parent(oi) + OrbitalSliceArray(Diagonal{C}(undef, length(oi)), (i, j)) +end -#endregion +function OrbitalSliceMatrix{C}(::UndefInitializer, oi::SparseIndices{<:Hamiltonian}) where {C} + h = hamiltonian(oi) + mat = mortar([flat(har) for har in harmonics(h), _ in 1:1]) + i, j = ham_to_orbslices(h) + return OrbitalSliceArray(mat, (i, j)) +end -############################################################################################ -# Hartree and Fock -#region +orbaxes(a::OrbitalSliceArray) = a.orbaxes -# struct Hartree{D<:DensityMatrix} -# rho:: -# end +Base.parent(a::OrbitalSliceArray) = a.parent #endregion diff --git a/test/test_greenfunction.jl b/test/test_greenfunction.jl index e2fcffa4..a9164855 100644 --- a/test/test_greenfunction.jl +++ b/test/test_greenfunction.jl @@ -44,6 +44,11 @@ end g = LP.linear() |> hopping(0) + @onsite((;ω=0) -> ω)|> greenfunction; @test only(g[sites(SA[1],1)](1.0; ω = 0)) ≈ 1.0 atol = 0.000001 @test only(g[sites(SA[1],1)](1.0; ω = -1)) ≈ 0.5 atol = 0.000001 + gs = g[sites(SA[1],1)] + @test gs(0.2) !== gs(0.3) + @test Quantica.call!(gs, 0.2) === Quantica.call!(gs, 0.2) + gs = g[cells = SA[2]] + @test Quantica.call!(gs, 0.2) === Quantica.call!(gs, 0.2) end @testset "greenfunction with contacts" begin @@ -178,7 +183,7 @@ end attach(nothing, region = RP.circle(1)) |> greenfunction(GS.KPM(order = 500)) gωs = g[1](ω) ρflat´ = -imag.(diag(gωs))/pi - ρ´ = ldos(g[1])(ω) + ρ´ = ldos(g[1], kernel = I)(ω) @test all(<(0.01), abs.(ρ .- ρ´)) @test all(<(0.01), abs.(ρflat .- ρflat´)) @test_throws ArgumentError g[sites((), 3:4)](ω) @@ -251,7 +256,7 @@ end # exercise band digagonal fastpath g = h |> greenfunction(GS.Bands(ϕs, ϕs)) - @test g(-0.2)[diagonal(sites(SA[1,2], :))] isa AbstractVector + @test g(-0.2)[diagonal(sites(SA[1,2], :))] isa AbstractMatrix end @testset "greenfunction 32bit" begin @@ -284,52 +289,67 @@ end end end -@testset "greenfunction diagonal slicing" begin +@testset "greenfunction sparse slicing" begin g = HP.graphene(a0 = 1, t0 = 1, orbitals = (2,1)) |> supercell((2,-2), region = r -> 0<=r[2]<=5) |> attach(nothing, cells = 2, region = r -> 0<=r[2]<=2) |> attach(nothing, cells = 3) |> greenfunction(GS.Schur(boundary = -2)) ω = 0.6 - @test g[diagonal(2)](ω) isa OrbitalSliceVector + @test g[diagonal(2)](ω) isa OrbitalSliceMatrix @test g[diagonal(2)](ω) == g(ω)[diagonal(2)] - @test size(g[diagonal(1)](ω)) == (12,) - @test size(g[diagonal(1, kernel = I)](ω)) == (8,) - @test length(g[diagonal(:)](ω)) == length(g[diagonal(1)](ω)) + length(g[diagonal(2)](ω)) == 48 - @test length(g[diagonal(:, kernel = I)](ω)) == length(g[diagonal(1, kernel = I)](ω)) + length(g[diagonal(2, kernel = I)](ω)) == 32 - @test length(g[diagonal(cells = 2:3)](ω)) == 72 - @test length(g[diagonal(cells = 2:3, kernel = I)](ω)) == 48 - @test ldos(g[1])(ω) ≈ -imag.(g[diagonal(1; kernel = I)](ω)) ./ π - + @test g[diagonal(:)](ω) == g(ω)[diagonal(:)] + @test size(g[diagonal(1)](ω)) == (12,12) + @test size(g[diagonal(1, kernel = I)](ω)) == (8,8) + @test size(g[diagonal(:)](ω),1) == size(g[diagonal(1)](ω),1) + size(g[diagonal(2)](ω),1) == 48 + @test size(g[diagonal(:, kernel = I)](ω),1) == size(g[diagonal(1, kernel = I)](ω),1) + size(g[diagonal(2, kernel = I)](ω),1) == 32 + @test size(g[diagonal(cells = 2:3)](ω)) == (72, 72) + @test size(g[diagonal(cells = 2:3, kernel = I)](ω)) == (48, 48) + @test ldos(g[1], kernel = I)(ω) ≈ -imag.(diag(g[diagonal(1; kernel = I)](ω))) ./ π + gs = g[diagonal(cells = 2:3)] + @test Quantica.call!(gs, ω) === Quantica.call!(gs, 2ω) inds = (1, :, (; cells = 1, sublats = :B), sites(2, 1:3), sites(0, 2), sites(3, :)) - for i in inds, j in inds - @test g(0.2)[diagonal(i)] isa Union{OrbitalSliceVector, AbstractVector} + for i in inds + @test g(0.2)[diagonal(i)] isa Union{OrbitalSliceMatrix, AbstractMatrix} end + gg = g(ω)[sitepairs(range = 1/sqrt(3))] + @test gg isa OrbitalSliceMatrix{ComplexF64,Quantica.SparseMatrixCSC{ComplexF64, Int64}} + @test size(gg) == Quantica.flatsize(g) .* (3, 1) + @test_throws DimensionMismatch g(ω)[sitepairs(range = 1, sublats = :B => :B, kernel = SA[1 0; 0 -1])] + gg = g(ω)[sitepairs(range = 1, sublats = :A => :A, kernel = I)] + @test size(gg) == size(g, 1) .* (3, 1) + + g = HP.graphene(a0 = 1, t0 = 1, orbitals = (2,1)) |> supercell((1, -1)) |> + greenfunction(GS.Schur(boundary = -2)) + ω = 0.6 + gg = g(ω)[sitepairs(range = 1)] + gg´ = g(ω)[orbaxes(gg)...] + @test gg .* gg´ ≈ gg .* gg end -@testset "densitymatrix diagonal slicing" begin +@testset "densitymatrix sparse slicing" begin g1 = LP.honeycomb() |> hamiltonian(hopping(0.1I, sublats = (:A,:B) .=> (:A,:B), range = 1) - plusadjoint(@hopping((; t=1) -> t*SA[0.4 1], sublats = :B => :A)), orbitals = (1,2)) |> supercell((1,-1), region = r -> -2<=r[2]<=2) |> attach(nothing, region = RP.circle(1), sublats = :B) |> attach(nothing, sublats = :A, cells = 0, region = r->r[2]<0) |> greenfunction(GS.Schur()) g2 = LP.square() |> hopping(1) - onsite(1) |> supercell(3) |> supercell |> attach(nothing, region = RP.circle(2)) |> greenfunction(GS.Spectrum()) g3 = LP.triangular() |> hamiltonian(hopping(-I) + onsite(1.8I), orbitals = 2) |> supercell(10) |> supercell |> attach(nothing, region = RP.circle(2)) |> attach(nothing, region = RP.circle(2, SA[1,2])) |> greenfunction(GS.KPM(bandrange=(-4,5))) - for g in (g1, g2, g3) - ρ0 = densitymatrix(g[diagonal(1)], 5) - ρ = densitymatrix(g[diagonal(1)]) - @test g[diagonal(1)](0.2) == g(0.2)[diagonal(1)] - @test g[diagonal(:)](0.2) == g(0.2)[diagonal(:)] - @test g[diagonal(1)](0.2) isa OrbitalSliceVector - @test isapprox(ρ0(), ρ()) - @test isapprox(ρ0(0.2), ρ(0.2)) - if g !== g3 # KPM doesn't support finite temperatures yet - @test isapprox(ρ0(0.2, 0.3), ρ(0.2, 0.3)) - end - if g === g1 || g === g3 - ρ0 = densitymatrix(g[diagonal(1, kernel = SA[0 1; 1 0])], 5) - ρ = densitymatrix(g[diagonal(1, kernel = SA[0 1; 1 0])]) - @test isapprox(ρ0(), ρ(); atol = 1e-8) - @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-8) - end + # KPM doesn't support finite temperatures or generic indexing yet, so g3 excluded from this loop + for g in (g1, g2), inds in (diagonal(1), diagonal(:), sitepairs(range = 1)) + ρ0 = densitymatrix(g[inds], 5) + ρ = densitymatrix(g[inds]) + @test isapprox(ρ0(), ρ(); atol = 1e-7) + @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-7) + @test isapprox(ρ0(0.2, 0.3), ρ(0.2, 0.3)) end + for g in (g1, g3) + ρ0 = densitymatrix(g[diagonal(1, kernel = SA[0 1; 1 0])], 5) + ρ = densitymatrix(g[diagonal(1, kernel = SA[0 1; 1 0])]) + @test isapprox(ρ0(), ρ(); atol = 1e-8) + @test isapprox(ρ0(0.2), ρ(0.2); atol = 1e-8) + end + g = HP.graphene(orbitals = 2) |> supercell((1,-1)) |> greenfunction + ρ0 = densitymatrix(g[sitepairs(range = 1)], 6) + mat = ρ0(0.2, 0.3) + @test Quantica.nnz(parent(mat)) < length(parent(mat)) # no broken sparsity end @testset "OrbitalSliceArray slicing" begin @@ -359,6 +379,9 @@ end c = sites(SA[1], 1) view(gmat, c) @test (@allocations view(gmat, c)) <= 2 + i, j = orbaxes(gmat) + @test g(0.2)[i, j] isa Quantica.OrbitalSliceMatrix + @test gmat[c] isa Matrix end function testcond(g0; nambu = false) @@ -466,6 +489,8 @@ end g0 = LP.square() |> hamiltonian(hopping(1)) |> supercell(region = RP.square(10)) |> attach(glead, reverse = true; region = contact2) |> attach(glead; region = contact1) |> greenfunction; testcond(g0) + contact1 = r -> r[1] ≈ 5 && -1 <= r[2] <= 1 + contact2 = r -> r[1] ≈ -5 && -1 <= r[2] <= 1 glead = LP.square() |> hamiltonian(hopping(I) + onsite(SA[0 1; 1 0]), orbitals = 2) |> supercell((1,0), region = r -> -1 <= r[2] <= 1) |> greenfunction(GS.Schur(boundary = 0)); g0 = LP.square() |> hamiltonian(hopping(I), orbitals = 2) |> supercell(region = RP.square(10)) |> attach(glead, reverse = true; region = contact2) |> attach(glead; region = contact1) |> greenfunction; testcond(g0; nambu = true)