diff --git a/docs/devdocs/singleshot.md b/docs/devdocs/singleshot.md new file mode 100644 index 00000000..d2a70c36 --- /dev/null +++ b/docs/devdocs/singleshot.md @@ -0,0 +1,80 @@ +# SingleShot1DGreensSolver + +Goal: solve the semi-infinite Green's function [g₀⁻¹ -V'GV]G = 1, where GV = ∑ₖφᵢᵏ λᵏ (φ⁻¹)ᵏⱼ. +It involves solving retarded solutions to + [(ωI-H₀)λ - V - V'λ²] φ = 0 +We do a SVD of + V = WSU' + V' = US'W' +where U and W are unitary, and S is diagonal with perhaps some zero +singlular values. We write [φₛ, φᵦ] = U'φ, where the β sector has zero singular values, + SPᵦ = 0 + U' = [Pₛ; Pᵦ] + [φₛ; φᵦ] = U'φ = [Pₛφ; Pᵦφ] + [Cₛₛ Cₛᵦ; Cᵦₛ Cᵦᵦ] = U'CU for any operator C +Then + [λU'(ωI-H₀)U - U'VU - λ²U'VU] Uφ = 0 + U'VU = U'WS = [Vₛₛ 0; Vᵦₛ 0] + U'V'U= S'W'U= [Vₛₛ' Vᵦₛ; 0 0] + U'(ωI-H₀)U = U'(g₀⁻¹)U = [g₀⁻¹ₛₛ g₀⁻¹ₛᵦ; g₀⁻¹ᵦₛ g₀⁻¹ᵦᵦ] = ωI - [H₀ₛₛ H₀ₛᵦ; H₀ᵦₛ H₀ᵦᵦ] + +The [λ(ωI-H₀) - V - λ²V'] φ = 0 problem can be solved by defining φχ = [φ,χ] = [φ, λφ] and +solving eigenpairs A*φχ = λ B*φχ, where + A = [0 I; -V g₀⁻¹] + B = [I 0; 0 V'] +To eliminate zero or infinite λ, we reduce the problem to the s sector + Aₛₛ [φₛ,χₛ] = λ Bₛₛ [φₛ,χₛ] + Aₛₛ = [0 I; -Vₛₛ g₀⁻¹ₛₛ] + [0 0; -Σ₁ -Σ₀] + Bₛₛ = [I 0; 0 Vₛₛ'] + [0 0; 0 Σ₁'] +where + Σ₀ = g₀⁻¹ₛᵦ g₀ᵦᵦ g₀⁻¹ᵦₛ + V'ₛᵦ g₀ᵦᵦ Vᵦₛ = H₀ₛᵦ g₀ᵦᵦ H₀ᵦₛ + Vₛᵦ g₀ᵦᵦ Vᵦₛ + Σ₁ = - g₀⁻¹ₛᵦ g₀ᵦᵦ Vᵦₛ = H₀ₛᵦ g₀ᵦᵦ Vᵦₛ + g₀⁻¹ₛₛ = ωI - H₀ₛₛ + g₀⁻¹ₛᵦ = - H₀ₛᵦ +Here g₀ᵦᵦ is the inverse of the bulk part of g₀⁻¹, g₀⁻¹ᵦᵦ = ωI - H₀ᵦᵦ. To compute this inverse +efficiently, we store the Hessenberg factorization `hessbb = hessenberg(-H₀ᵦᵦ)` and use shifts. +Then, g₀ᵦᵦ H₀ᵦₛ = (hess + ω I) \ H₀ᵦₛ and g₀ᵦᵦ Vᵦₛ = (hess + ω I) \ Vᵦₛ. +Diagonalizing (Aₛₛ, Bₛₛ) we obtain the surface projection φₛ = Pₛφ of eigenvectors φ. The +bulk part is + φᵦ = g₀ᵦᵦ (λ⁻¹Vᵦₛ - g₀⁻¹ᵦₛ) φₛ = g₀ᵦᵦ(λ⁻¹Vᵦₛ + H₀ᵦₛ) φₛ +so that the full φ's with non-zero λ read + φ = U[φₛ; φᵦ] = [Pₛ' Pᵦ'][φₛ; φᵦ] = [Pₛ' + Pᵦ' g₀ᵦᵦ(λ⁻¹Vᵦₛ + H₀ᵦₛ)] φₛ = Z[λ] φₛ + Z[λ] = [Pₛ' + Pᵦ' g₀ᵦᵦ(λ⁻¹Vᵦₛ + H₀ᵦₛ)] +We can complete the set with the λ=0 solutions, Vφᴿ₀ = 0 and the λ=∞ solutions V'φᴬ₀ = 0. +Its important to note that U'φᴿ₀ = [0; φ₀ᵦᴿ] and W'φᴬ₀ = [0; φ₀ᵦ´ᴬ] + +By computing the velocities vᵢ = im * φᵢ'(V'λᵢ-Vλᵢ')φᵢ / φᵢ'φᵢ = 2*imag(χᵢ'Vφᵢ)/φᵢ'φᵢ we +classify φ into retarded (|λ|<1 or v > 0) and advanced (|λ|>1 or v < 0). Then + φᴿ = Z[λᴿ]φₛᴿ = [φₛᴿ; φᵦᴿ] + U'φᴿ = [φₛᴿ; φᵦᴿ] + φᴬ = Z[λᴬ]φₛᴬ = [φₛᴬ; φᵦᴬ] + Wφᴬ = [φₛ´ᴬ; φᵦ´ᴬ] (different from [φₛᴬ; φᵦᴬ]) +The square matrix of all retarded and advanced modes read [φᴿ φᴿ₀] and [φᴬ φᴬ₀]. +We now return to the Green's functions. The right and left semiinfinite GrV and GlV' read + GrV = [φᴿ φ₀ᴿ][λ 0; 0 0][φᴿ φ₀ᴿ]⁻¹ = φᴿ λ (φₛᴿ)⁻¹Pₛ + (used [φᴿ φ₀ᴿ] = U[φₛᴿ 0; φᵦᴿ φᴿ₀ᵦ] => [φᴿ φ₀ᴿ]⁻¹ = [φₛᴿ⁻¹ 0; -φᵦᴿ(φₛᴿ)⁻¹ (φᴿ₀ᵦ)⁻¹]U') + GlV´ = [φᴬ φ₀ᴬ][(λᴬ)⁻¹ 0; 0 0][φᴬ φ₀ᴬ]⁻¹ = φᴬ λ (φₛ´ᴬ)⁻¹Pₛ´ + (used [φᴬ φᴬ₀]⁻¹ = W[φₛ´ᴬ 0; φᵦ´ᴬ φ₀ᵦ´ᴬ] => [φᴬ φ₀ᴬ]⁻¹ = [(φₛ´ᴬ)⁻¹ 0; -φᵦ´ᴬ(φₛ´ᴬ)⁻¹ (φ₀ᵦ´ᴬ)⁻¹]W') + +We can then write the local Green function G∞_{0} of the semiinfinite and infinite lead as + Gr_{1,1} = Gr = [g₀⁻¹ - V'GrV]⁻¹ + Gl_{-1,-1} = Gl = [g₀⁻¹ - VGlV']⁻¹ + G∞_{0} = [g₀⁻¹ - VGlV' - V'GrV]⁻¹ + (GrV)ᴺ = φᴿ (λᴿ)ᴺ (φₛᴿ)⁻¹ Pₛ = χᴿ (λᴿ)ᴺ⁻¹ (φₛᴿ)⁻¹ Pₛ + (GlV´)ᴺ = φᴬ (λᴬ)⁻ᴺ (φₛᴬ)⁻¹ Pₛₚ = φᴬ (λᴬ)¹⁻ᴺ(χₛᴬ)⁻¹ Pₛₚ + φᴿ = [Pₛ' + Pᵦ' g₀ᵦᵦ ((λᴿ)⁻¹Vᵦₛ + H₀ᵦₛ)]φₛᴿ + φᴬ = [Pₛ' + Pᵦ' g₀ᵦᵦ ((λᴬ)Vᵦₛ + H₀ᵦₛ)]φₛᴬ +where φₛᴿ and φₛᴬ are retarded/advanced eigenstates with eigenvalues λᴿ and λᴬ. Here Pᵦ is +the projector onto the uncoupled V sites, VPᵦ = 0, Pₛ = 1 - Pᵦ is the surface projector of V +and Pₛₚ is the surface projector of V'. + +Defining + GVᴺ = (GrV)ᴺ for N>0 and (GlV´)⁻ᴺ for N<0 +we have + Semiinifinite: G_{N,M} = (GVᴺ⁻ᴹ - GVᴺGV⁻ᴹ)G∞_{0} + Infinite: G∞_{N} = GVᴺ G∞_{0} +Spelling it out + Semiinfinite right (N,M > 0): G_{N,M} = G∞_(N-M) - (GrV)ᴺ G∞_(-M) = [GVᴺ⁻ᴹ - (GrV)ᴺ (GlV´)ᴹ]G∞_0. + At surface (N=M=1), G_{1,1} = Gr = (1- GrVGlV´)G∞_0 = [g₀⁻¹ - V'GrV]⁻¹ + Semiinfinite left (N,M < 0): G_{N,M} = G∞_(N-M) - (GlV´)⁻ᴺ G∞_(-M) = [GVᴺ⁻ᴹ - (GlV´)⁻ᴺ(GrV)⁻ᴹ]G∞_0. + At surface (N=M=-1), G_{1,1} = Gl = (1- GrVGlV´)G∞_0 = [g₀⁻¹ - VGlV']⁻¹ diff --git a/src/Quantica.jl b/src/Quantica.jl index d89474d6..a5ff72b2 100644 --- a/src/Quantica.jl +++ b/src/Quantica.jl @@ -29,7 +29,7 @@ export sublat, bravais, lattice, dims, supercell, unitcell, bands, vertices, energies, states, degeneracy, momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM, - greens, greensolver + greens, greensolver, SingleShot1D export RegionPresets, RP, LatticePresets, LP, HamiltonianPresets, HP diff --git a/src/bandstructure.jl b/src/bandstructure.jl index 3f4cb1a0..ed4130e9 100644 --- a/src/bandstructure.jl +++ b/src/bandstructure.jl @@ -164,7 +164,7 @@ subspace(s::Spectrum, rngs) = subspace.(Ref(s), rngs) function subspace(s::Spectrum, rng::AbstractUnitRange) ε = mean(j -> s.energies[j], rng) ψs = view(s.states, :, rng) - Subspace(s.diag.orbstruct, ε, ψs) + return Subspace(s.diag.orbstruct, ε, ψs) end nsubspaces(s::Spectrum) = length(s.subs) diff --git a/src/diagonalizer.jl b/src/diagonalizer.jl index 8d232d19..3f5d010f 100644 --- a/src/diagonalizer.jl +++ b/src/diagonalizer.jl @@ -14,9 +14,10 @@ end struct Diagonalizer{M<:AbstractDiagonalizeMethod,F,O<:Union{OrbitalStructure,Missing}} method::M - perm::Vector{Int} # reusable permutation vector - matrixf::F # functor or function matrixf(φs) that produces matrices to be diagonalized - orbstruct::O # store structure of original Hamiltonian if available (to allow unflattening eigenstates) + matrixf::F # functor or function matrixf(φs) that produces matrices to be diagonalized + orbstruct::O # store structure of original Hamiltonian if available (to allow unflattening eigenstates) + perm::Vector{Int} # reusable permutation vector + perm´::Vector{Int} # reusable permutation vector end struct NoUnflatten end @@ -57,17 +58,17 @@ function diagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}; method = Line matrixf = HamiltonianBlochFunctor(h, matrix, mapping) perm = Vector{Int}(undef, size(matrix, 2)) orbstruct = parent(h).orbstruct - return Diagonalizer(method, perm, matrixf, orbstruct) + return Diagonalizer(method, matrixf, orbstruct, perm, copy(perm)) end function diagonalizer(matrixf::Function, dimh; method = LinearAlgebraPackage()) perm = Vector{Int}(undef, dimh) - return Diagonalizer(method, perm, matrixf, missing) + return Diagonalizer(method, matrixf, missing, perm, copy(perm)) end @inline function (d::Diagonalizer)(φs, ::NoUnflatten) ϵ, ψ = diagonalize(d.matrixf(φs), d.method) - issorted(ϵ, by = real) || sorteigs!(d.perm, ϵ, ψ) + issorted(ϵ, by = real) || sorteigs!(ϵ, ψ, d) fixphase!(ψ) return ϵ, ψ end @@ -82,11 +83,13 @@ end @inline (d::Diagonalizer)() = d(()) -function sorteigs!(perm, ϵ::AbstractVector, ψ::AbstractMatrix) - resize!(perm, length(ϵ)) - p = sortperm!(perm, ϵ, by = real, alg = Base.DEFAULT_UNSTABLE) - permute!(ϵ, p) - Base.permutecols!!(ψ, p) +function sorteigs!(ϵ::AbstractVector, ψ::AbstractMatrix, d) + p, p´ = d.perm, d.perm´ + resize!(p, length(ϵ)) + resize!(p´, length(ϵ)) + sortperm!(p, ϵ, by = real, alg = Base.DEFAULT_UNSTABLE) + Base.permute!!(ϵ, copy!(p´, p)) + Base.permutecols!!(ψ, copy!(p´, p)) return ϵ, ψ end diff --git a/src/greens.jl b/src/greens.jl index cdfa5d95..be9e108d 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -1,35 +1,44 @@ ####################################################################### # Green's function ####################################################################### -abstract type GreenSolver end +abstract type AbstractGreensSolver end -struct GreensFunction{S<:GreenSolver,H} - h::H +struct GreensFunction{S<:AbstractGreensSolver,L,B<:NTuple{L,Union{Int,Missing}},H<:Hamiltonian} solver::S + h::H + boundaries::B end """ - greens(h::Hamiltonian, solveobject) + greens(h::Hamiltonian, solveobject; boundaries::NTuple{L,Integer} = missing) -Construct the Green's function `g::GreensFunction` of `h` using the provided `solveobject`. -Currently valid `solveobject`s are +Construct the Green's function `g::GreensFunction` of `L`-dimensional Hamiltonian `h` using +the provided `solveobject`. Currently valid `solveobject`s are - the `Bandstructure` of `h` (for an unbounded `h` or an `Hamiltonian{<:Superlattice}}`) - the `Spectrum` of `h` (for a bounded `h`) +- `SingleShot1D(; direct = false)` (single-shot generalized [or direct if `direct = true`] eigenvalue approach for 1D Hamiltonians) + +If a `boundaries = (n₁, n₂, ...)` is provided, a reflecting boundary is assumed for each +non-missing `nᵢ` perpendicular to Bravais vector `i` at a cell distance `nᵢ` from the +origin. + + h |> greens(h -> solveobject(h), args...) - h |> greens(h -> solveobject(h)) +Curried form equivalent to the above, giving `greens(h, solveobject(h), args...)`. -Curried form equivalent to the above, giving `greens(h, solveobject(h))` (see -example below). + g(ω, cells::Pair) - g([m,] ω, cells::Pair = missing) +From a constructed `g::GreensFunction`, obtain the retarded Green's function matrix at +frequency `ω` between unit cells `src` and `dst` by calling `g(ω, src => dst)`, where `src, +dst` are `::NTuple{L,Int}` or `SVector{L,Int}`. If not provided, `cells` default to +`(1, 1, ...) => (1, 1, ...)`. -From a constructed `g::GreensFunction`, obtain the retarded Green's function -matrix at frequency `ω` between unit cells `src` and `dst` by calling `g(ω, src -=> dst)`, where `src, dst` are `::NTuple{L,Int}` or `SVector{L,Int}`. If -`cells` is missing, `src` and `dst` are assumed to be zero vectors. For -performance, one can use a preallocated matrix `m` (e.g. `m = -similarmatrix(h)`) by calling `g(m, ω, cells)`. + g(ω, missing) + +If allowed by the used `solveobject`, build an efficient function `cells -> g(ω, cells)` +that can produce the Greens function between different cells at fixed `ω` without repeating +cell-independent parts of the computation. # Examples @@ -49,16 +58,462 @@ julia> m = similarmatrix(g); g(m, 0.2) 6.663377810046025 - 24.472789025006396im ``` +# See also + `greens!`, `SingleShot1D` +""" +greens(h::Hamiltonian{<:Any,L}, solverobject; boundaries = filltuple(missing, Val(L))) where {L} = + GreensFunction(greensolver(solverobject, h), h, boundaries) +greens(solver::Function, args...; kw...) = h -> greens(solver(h), h, args...; kw...) + +# solver fallback +greensolver(s::AbstractGreensSolver) = s + +# call API fallback +(g::GreensFunction)(ω, cells = default_cells(g)) = greens!(similarmatrix(g), g, ω, cells) + +similarmatrix(g::GreensFunction, type = Matrix{blocktype(g.h)}) = similarmatrix(g.h, type) + +greens!(matrix, g, ω, cells) = greens!(matrix, g, ω, sanitize_cells(cells, g)) + +default_cells(::GreensFunction{S,L}) where {S,L} = filltuple(1, Val(L)) => filltuple(1, Val(L)) + +sanitize_cells((cell0, cell1)::Pair{<:Integer,<:Integer}, ::GreensFunction{S,1}) where {S} = + SA[cell0] => SA[cell1] +sanitize_cells((cell0, cell1)::Pair{<:NTuple{L,Integer},<:NTuple{L,Integer}}, ::GreensFunction{S,L}) where {S,L} = + SVector(cell0) => SVector(cell1) +sanitize_cells(cells, g::GreensFunction{S,L}) where {S,L} = + throw(ArgumentError("Cells should be of the form `cᵢ => cⱼ`, with each `c` an `NTuple{$L,Integer}`")) + +const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} + +####################################################################### +# SingleShot1DGreensSolver +####################################################################### +""" + SingleShot1D(; direct = true) + +Return a Greens function solver using the generalized eigenvalue approach, whereby given the +energy `ω`, the eigenmodes of the infinite 1D Hamiltonian, and the corresponding infinite +and semi-infinite Greens function can be computed by solving the generalized eigenvalue +equation + + A⋅φχ = λ B⋅φχ + A = [0 I; V ω-H0] + B = [I 0; 0 V'] + +This is the matrix form of the problem `λ(ω-H0)φ - Vφ - λ²V'φ = 0`, where `φχ = [φ; λφ]`, +and `φ` are `ω`-energy eigenmodes, with (possibly complex) momentum `q`, and eigenvalues are +`λ = exp(-iqa₀)`. The algorithm assumes the Hamiltonian has only `dn = (0,)` and `dn = (±1, +)` Bloch harmonics (`H0`, `V` and `V'`), so its unit cell will be enlarged before applying +the solver if needed. Bound states in the spectrum will yield delta functions in the density +of states that can be resolved by adding a broadening in the form of a small positive +imaginary part to `ω`. + +To avoid singular solutions `λ=0,∞`, the nullspace of `V` is projected out of the problem. +Sometimes, however, some linear algebra implementations of the eigensolver can lead to +incomplete convergence, signaled by a `LAPACKException(n)`. This problem does not arise if +we the generalized eigenproblem into a standard eigenproblem. This is achieved with `direct += true`. The performance is also slightly improved. Such approach, however, is not possible +for some special infinite 1D lattices that have bound states embedded in a continuum even in +the absence of boundaries. A `SingularExcepcion` error is then thrown. In such cases, try +`direct = false`. + +# Examples +```jldoctest +julia> using LinearAlgebra + +julia> h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), (10,10)) |> Quantica.wrap(2); + +julia> g = greens(h, SingleShot1D(), boundaries = (0,)) +GreensFunction{SingleShot1DGreensSolver}: Green's function using the single-shot 1D method + Matrix size : 40 × 40 + Element type : scalar (ComplexF64) + Boundaries : (0,) + Reduced dims : 20 + +julia> tr(g(0.3)) +-32.193416068730784 - 3.4399800712973474im +``` + +# See also + `greens` """ -greens(h, solver) = GreensFunction(h, greensolver(solver)) -greens(solver::Function) = h -> greens(h, solver(h)) +struct SingleShot1D + invert_B::Bool +end + +SingleShot1D(; direct = true) = SingleShot1D(direct) + +struct SingleShot1DTemporaries{T} + b2s::Matrix{T} # size = dim_b, 2dim_s + φ::Matrix{T} # size = dim_H0, 2dim_s + χ::Matrix{T} # size = dim_H0, 2dim_s + ss1::Matrix{T} # size = dim_s, dim_s + ss2::Matrix{T} # size = dim_s, dim_s + Hs1::Matrix{T} # size = dim_H0, dim_s + Hs2::Matrix{T} # size = dim_H0, dim_s + HH0::Matrix{T} # size = dim_H0, dim_H0 + HH1::Matrix{T} # size = dim_H0, dim_H0 + HH2::Matrix{T} # size = dim_H0, dim_H0 + HH3::Matrix{T} # size = dim_H0, dim_H0 + vH::Vector{T} # size = dim_H0 + v2s1::Vector{Int} # size = 2dim_s + v2s2::Vector{Int} # size = 2dim_s +end + +function SingleShot1DTemporaries{T}(dim_H, dim_s, dim_b) where {T} + b2s = Matrix{T}(undef, dim_b, 2dim_s) + φ = Matrix{T}(undef, dim_H, 2dim_s) + χ = Matrix{T}(undef, dim_H, 2dim_s) + ss1 = Matrix{T}(undef, dim_s, dim_s) + ss2 = Matrix{T}(undef, dim_s, dim_s) + Hs1 = Matrix{T}(undef, dim_H, dim_s) + Hs2 = Matrix{T}(undef, dim_H, dim_s) + HH0 = Matrix{T}(undef, dim_H, dim_H) + HH1 = Matrix{T}(undef, dim_H, dim_H) + HH2 = Matrix{T}(undef, dim_H, dim_H) + HH3 = Matrix{T}(undef, dim_H, dim_H) + vH = Vector{T}(undef, dim_H) + v2s1 = Vector{Int}(undef, 2dim_s) + v2s2 = Vector{Int}(undef, 2dim_s) + return SingleShot1DTemporaries{T}(b2s, φ, χ, ss1, ss2, Hs1, Hs2, HH0, HH1, HH2, HH3, vH, v2s1, v2s2) +end + +struct SingleShot1DGreensSolver{T<:Complex,R<:Real,H<:Hessenberg{T}} <: AbstractGreensSolver + invert_B::Bool + A::Matrix{T} + B::Matrix{T} + minusH0::SparseMatrixCSC{T,Int} + V::SparseMatrixCSC{T,Int} + Pb::Matrix{T} # size = dim_b, dim_H0 + Ps::Matrix{T} # size = dim_s, dim_H0 + Ps´::Matrix{T} # size = dim_s´, dim_H0 + H0ss::Matrix{T} + H0bs::Matrix{T} + Vss::Matrix{T} + Vbs::Matrix{T} + hessbb::H + velocities::Vector{R} # size = 2dim_s = 2num_modes + maxdn::Int + temps::SingleShot1DTemporaries{T} +end + +function Base.show(io::IO, g::GreensFunction{<:SingleShot1DGreensSolver}) + print(io, summary(g), "\n", +" Matrix size : $(size(g.solver.V, 1)) × $(size(g.solver.V, 2)) + Reduced size : $(size(g.solver.Ps, 1)) × $(size(g.solver.Ps, 1)) + Element type : $(displayelements(g.h)) + Boundaries : $(g.boundaries)") +end + +Base.summary(g::GreensFunction{<:SingleShot1DGreensSolver}) = + "GreensFunction{SingleShot1DGreensSolver}: Green's function using the single-shot 1D method" + +hasbulk(gs::SingleShot1DGreensSolver) = !iszero(size(gs.Pb, 1)) + +function greensolver(s::SingleShot1D, h) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) + return SingleShot1DGreensSolver(h, s.invert_B) +end + +## Preparation + +function SingleShot1DGreensSolver(h::Hamiltonian, invert_B) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) + maxdn = max(1, maximum(har -> abs(first(har.dn)), h.harmonics)) + H = flatten(maxdn == 1 ? h : unitcell(h, (maxdn,))) + H0, V, V´ = H[(0,)], H[(1,)], H[(-1,)] + Pb, Ps, Ps´ = bulk_surface_projectors(H0, V, V´) + H0ss = Ps * H0 * Ps' + H0bs = Pb * H0 * Ps' + Vss = Ps * V * Ps' + Vbs = Pb * V * Ps' + hessbb = hessenberg!(Pb * (-H0) * Pb') + dim_s = size(Ps, 1) + T = complex(blockeltype(H)) + A = zeros(T, 2dim_s, 2dim_s) + B = zeros(T, 2dim_s, 2dim_s) + dim_s, dim_b, dim_H = size(Ps, 1), size(Pb, 1), size(H0, 2) + velocities = fill(zero(real(T)), 2dim_s) + temps = SingleShot1DTemporaries{T}(dim_H, dim_s, dim_b) + return SingleShot1DGreensSolver(invert_B, A, B, -H0, V, Pb, Ps, Ps´, H0ss, H0bs, Vss, Vbs, hessbb, velocities, maxdn, temps) +end + +function bulk_surface_projectors(H0::AbstractMatrix{T}, V, V´) where {T} + SVD = svd(Matrix(V), full = true) + W, S, U = SVD.U, SVD.S, SVD.V + dim_b = count(s -> iszero(chop(s)), S) + dim_s = length(S) - dim_b + if iszero(dim_b) + Ps = Matrix{T}(I, dim_s, dim_s) + Ps´ = copy(Ps) + Pb = Ps[1:0, :] + else + Ps = U'[1:dim_s, :] + Pb = U'[dim_s+1:end, :] + Ps´ = W'[1:dim_s, :] + Pb´ = W'[dim_s+1:end, :] + end + return Pb, Ps, Ps´ +end + +## Solver execution + +function (gs::SingleShot1DGreensSolver)(ω) + A, B = single_shot_surface_matrices(gs, ω) + λs, φχs = eigen_funcbarrier(A, B, gs.invert_B) + dim_s = size(φχs, 1) ÷ 2 + φs = view(φχs, 1:dim_s, :) + χs = view(φχs, dim_s+1:2dim_s, :) + φ = gs.temps.φ + χ = gs.temps.χ + if hasbulk(gs) + reconstruct_bulk!(φ, ω, λs, φs, gs) + reconstruct_bulk!(χ, ω, λs, χs, gs) + else + copy!(φ, φs) + copy!(χ, χs) + end + return classify_retarded_advanced(λs, φs, χs, φ, χ, gs) +end + +function eigen_funcbarrier(A::AbstractMatrix{T}, B, invert_B)::Tuple{Vector{T},Matrix{T}} where {T} + if invert_B + factB = lu!(B) + B⁻¹A = ldiv!(factB, A) + λs, φχs = eigen!(B⁻¹A; sortby = abs) + else + λs, φχs = eigen!(A, B; sortby = abs) + end + cleanup_λ!(λs) + return λs, φχs +end + +function cleanup_λ!(λs::AbstractVector{T}) where {T} + λs .= (λ -> ifelse(isnan(λ), T(Inf), ifelse(abs2(λ) < eps(real(T)), zero(T), λ))).(λs) + return λs +end + +function single_shot_surface_matrices(gs::SingleShot1DGreensSolver{T}, ω) where {T} + A, B = gs.A, gs.B + dim_s = size(gs.H0bs, 2) + fill!(A, 0) + fill!(B, 0) + for i in 1:dim_s + A[i, dim_s + i] = B[i, i] = one(T) + A[dim_s + i, dim_s + i] = ω + end + tmp = view(gs.temps.b2s, :, 1:dim_s) # gs.temps.b2s has 2dim_s cols + dim = size(A, 1) ÷ 2 + A21 = view(A, dim+1:2dim, 1:dim) + A22 = view(A, dim+1:2dim, dim+1:2dim) + B22 = view(B, dim+1:2dim, dim+1:2dim) + + # A22 = ωI - H₀ₛₛ - H₀ₛᵦ g₀ᵦᵦ H₀ᵦₛ - V'ₛᵦ g₀ᵦᵦ Vᵦₛ + A22 .-= gs.H0ss # ω already added to diagonal above + if hasbulk(gs) + copy!(tmp, gs.H0bs) + ldiv!(gs.hessbb + ω*I, tmp) + mul!(A22, gs.H0bs', tmp, -1, 1) + copy!(tmp, gs.Vbs) + ldiv!(gs.hessbb + ω*I, tmp) + mul!(A22, gs.Vbs', tmp, -1, 1) + end + + # A21 = - Vₛₛ - H₀ₛᵦ g₀ᵦᵦ Vᵦₛ + A21 .= .- gs.Vss + if hasbulk(gs) + copy!(tmp, gs.Vbs) + ldiv!(gs.hessbb + ω*I, tmp) + mul!(A21, gs.H0bs', tmp, -1, 1) + end + + # B22 = -A21' + B22 .= .- A21' + + chkfinite(A) + chkfinite(B) + + return A, B +end + +function chkfinite(A::AbstractMatrix) + for a in A + if !isfinite(a) + throw(ArgumentError("Matrix contains Infs or NaNs. This may happen when the energy ω exactly matches a bound state in the spectrum. Try adding a small positive imaginary part to ω.")) + end + end + return true +end + +# φ = [Pₛ' + Pᵦ' g₀ᵦᵦ (λ⁻¹Vᵦₛ + H₀ᵦₛ)] φₛ +function reconstruct_bulk!(φ, ω, λs, φs, gs) + tmp = gs.temps.b2s + mul!(tmp, gs.Vbs, φs) + tmp ./= transpose(λs) + mul!(tmp, gs.H0bs, φs, 1, 1) + ldiv!(gs.hessbb + ω*I, tmp) + mul!(φ, gs.Pb', tmp) + mul!(φ, gs.Ps', φs, 1, 1) + return φ +end + +# function classify_retarded_advanced(λs, φs, φ, χ, gs) +function classify_retarded_advanced(λs, φs, χs, φ, χ, gs) + vs = compute_velocities_and_normalize!(φ, χ, λs, gs) + # order for ret-evan, ret-prop, adv-prop, adv-evan + p = sortperm!(gs.temps.v2s1, vs; rev = true) + p´ = gs.temps.v2s2 + Base.permute!!(vs, copy!(p´, p)) + Base.permute!!(λs, copy!(p´, p)) + Base.permutecols!!(φ, copy!(p´, p)) + Base.permutecols!!(χ, copy!(p´, p)) + Base.permutecols!!(φs, copy!(p´, p)) + + ret, adv = nonsingular_ret_adv(λs, vs) + # @show vs + # @show abs.(λs) + # @show ret, adv, length(λs), length(vs) + + λR = view(λs, ret) + χR = view(χ, :, ret) # This resides in part of gs.temps.χ + φR = view(φ, :, ret) # This resides in part of gs.temps.φ + # overwrite output of eigen to preserve normalization of full φ + φRs = mul!(view(φs, :, ret), gs.Ps, φR) + iφRs = issquare(φRs) ? rdiv!(copyto!(gs.temps.ss1, I), lu!(φRs)) : pinv(φRs) + + iλA = view(λs, adv) + iλA .= inv.(iλA) + χA = view(χ, :, adv) # This resides in part of gs.temps.χ + φA = view(φ, :, adv) # This resides in part of gs.temps.φ + # overwrite output of eigen to preserve normalization of full χ + χAs´ = mul!(view(χs, :, adv), gs.Ps´, χA) + iχAs´ = issquare(χAs´) ? rdiv!(copyto!(gs.temps.ss2, I), lu!(χAs´)) : pinv(χAs´) + + return λR, χR, iφRs, iλA, φA, iχAs´ +end + +# The velocity is given by vᵢ = im * φᵢ'(V'λᵢ-Vλᵢ')φᵢ / φᵢ'φᵢ = 2*imag(χᵢ'Vφᵢ)/φᵢ'φᵢ +function compute_velocities_and_normalize!(φ, χ, λs, gs) + vs = gs.velocities + tmp = gs.temps.vH + for (i, λ) in enumerate(λs) + abs2λ = abs2(λ) + if abs2λ ≈ 1 + φi = view(φ, :, i) + χi = view(χ, :, i) + norm2φi = dot(φi, φi) + mul!(tmp, gs.V, φi) + v = 2*imag(dot(χi, tmp))/norm2φi + φi .*= inv(sqrt(norm2φi * abs(v))) + χi .*= inv(sqrt(norm2φi * abs(v))) + vs[i] = v + else + vs[i] = abs2λ < 1 ? Inf : -Inf + end + end # sortperm(vs) would now give the order of adv-evan, adv-prop, ret-prop, ret-evan + return view(vs, 1:length(λs)) +end + +function nonsingular_ret_adv(λs::AbstractVector{T}, vs) where {T} + rmin, rmax = 1, 0 + amin, amax = 1, 0 + ε = eps(real(T)) + for (i, v) in enumerate(vs) + aλ2 = abs2(λs[i]) + if aλ2 < ε + rmin = i + 1 + elseif v > 0 + rmax = i + amin = i + 1 + elseif v < 0 && isfinite(aλ2) + amax = i + end + end + return rmin:rmax, amin:amax +end + +## Greens execution + +(g::GreensFunction{<:SingleShot1DGreensSolver})(ω, cells) = g(ω, missing)(sanitize_cells(cells, g)) + +# Infinite: G∞_{N} = GVᴺ G∞_{0} +function (g::GreensFunction{<:SingleShot1DGreensSolver,1,Tuple{Missing}})(ω, ::Missing) + gs = g.solver + factors = gs(ω) + G∞⁻¹ = inverse_G∞!(gs.temps.HH0, ω, gs, factors) |> lu! + return cells -> G_infinite!(gs, factors, G∞⁻¹, cells) +end + +# Semiinifinite: G_{N,M} = (GVᴺ⁻ᴹ - GVᴺGV⁻ᴹ)G∞_{0} +function (g::GreensFunction{<:SingleShot1DGreensSolver,1,Tuple{Int}})(ω, ::Missing) + gs = g.solver + factors = gs(ω) + G∞⁻¹ = inverse_G∞!(gs.temps.HH0, ω, gs, factors) |> lu! + N0 = only(g.boundaries) + return cells -> G_semiinfinite!(gs, factors, G∞⁻¹, shift_cells(cells, N0)) +end + +function G_infinite!(gs, factors, G∞⁻¹, (src, dst)) + src´ = div(only(src), gs.maxdn, RoundToZero) + dst´ = div(only(dst), gs.maxdn, RoundToZero) + N = dst´ - src´ + GVᴺ = GVᴺ!(gs.temps.HH1, N, gs, factors) + G∞ = rdiv!(GVᴺ, G∞⁻¹) + return G∞ +end + +function G_semiinfinite!(gs, factors, G∞⁻¹, (src, dst)) + M = div(only(src), gs.maxdn, RoundToZero) + N = div(only(dst), gs.maxdn, RoundToZero) + if sign(N) != sign(M) + G∞ = fill!(gs.temps.HH1, 0) + else + GVᴺ⁻ᴹ = GVᴺ!(gs.temps.HH1, N-M, gs, factors) + GVᴺ = GVᴺ!(gs.temps.HH2, N, gs, factors) + GV⁻ᴹ = GVᴺ!(gs.temps.HH3, -M, gs, factors) + mul!(GVᴺ⁻ᴹ, GVᴺ, GV⁻ᴹ, -1, 1) # (GVᴺ⁻ᴹ - GVᴺGV⁻ᴹ) + G∞ = rdiv!(GVᴺ⁻ᴹ , G∞⁻¹) + end + return G∞ +end -# Needed to make similarmatrix work with GreensFunction -matrixtype(g::GreensFunction) = Matrix{eltype(g.h)} -Base.parent(g::GreensFunction) = g.h +shift_cells((src, dst), N0) = (only(src) - N0, only(dst) - N0) + +# G∞⁻¹ = G₀⁻¹ - V´GrV - VGlV´ with V´GrV = V'*χR*φRs⁻¹Ps and VGlV´ = iχA*φAs´⁻¹Ps´ +function inverse_G∞!(matrix, ω, gs, (λR, χR, iφRs, iλA, φA, iχAs´)) + G0⁻¹ = inverse_G0!(matrix, ω, gs) + V´GrV = mul!(gs.temps.Hs2, gs.V', mul!(gs.temps.Hs1, χR, iφRs)) + mul!(G0⁻¹, V´GrV, gs.Ps, -1, 1) + VGlV´ = mul!(gs.temps.Hs2, gs.V, mul!(gs.temps.Hs1, φA, iχAs´)) + G∞⁻¹ = mul!(G0⁻¹, VGlV´, gs.Ps´, -1, 1) + return G∞⁻¹ +end + +# G₀⁻¹ = ω - H₀ +function inverse_G0!(G0⁻¹, ω, gs) + copy!(G0⁻¹, gs.minusH0) + for i in axes(G0⁻¹, 1) + G0⁻¹[i, i] += ω + end + return G0⁻¹ +end + +function GVᴺ!(GVᴺ, N, gs, (λR, χR, iφRs, iλA, φA, iχAs´)) + if N == 0 + copyto!(GVᴺ, I) + elseif N > 0 + χᴿλᴿᴺ´ = N == 1 ? χR : (gs.temps.Hs1 .= χR .* transpose(λR) .^ (N-1)) + mul!(GVᴺ, mul!(gs.temps.Hs2, χᴿλᴿᴺ´, iφRs), gs.Ps) + else + φᴬλᴬᴺ´ = N == -1 ? φA : (gs.temps.Hs1 .= φA .* transpose(iλA) .^ (-N-1)) + mul!(GVᴺ, mul!(gs.temps.Hs2, φᴬλᴬᴺ´, iχAs´), gs.Ps´) + end + return GVᴺ +end ####################################################################### -# BandGreenSolver +# BandGreensSolver ####################################################################### struct SimplexData{D,E,T,C<:SMatrix,DD,SA<:SubArray} ε0::T @@ -75,25 +530,26 @@ struct SimplexData{D,E,T,C<:SMatrix,DD,SA<:SubArray} φs::NTuple{D,SA} end -struct BandGreenSolver{P<:SimplexData,E} <: GreenSolver +struct BandGreensSolver{P<:SimplexData,E,H<:Hamiltonian} <: AbstractGreensSolver simplexdata::Vector{P} indsedges::NTuple{E,Tuple{Int,Int}} # all distinct pairs of 1:V, where V=D+1=num verts + h::H end -function Base.show(io::IO, g::GreensFunction{<:BandGreenSolver}) +function Base.show(io::IO, g::GreensFunction{<:BandGreensSolver}) print(io, summary(g), "\n", " Matrix size : $(size(g.h, 1)) × $(size(g.h, 2)) Element type : $(displayelements(g.h)) Band simplices : $(length(g.solver.simplexdata))") end -Base.summary(g::GreensFunction{<:BandGreenSolver}) = +Base.summary(g::GreensFunction{<:BandGreensSolver}) = "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.h))D bandstructure" -function greensolver(b::Bandstructure{D}) where {D} +function greensolver(b::Bandstructure{D}, h) where {D} indsedges = tuplepairs(Val(D)) v = [SimplexData(simplex, band, indsedges) for band in bands(b) for simplex in band.simplices] - return BandGreenSolver(v, indsedges) + return BandGreensSolver(v, indsedges, h) end edges_per_simplex(L) = binomial(L,2) @@ -168,24 +624,16 @@ end ## Call API -(g::GreensFunction{<:BandGreenSolver})(ω::Number, cells = missing) = g(similarmatrix(g), ω, cells) - -function (g::GreensFunction{<:BandGreenSolver})(matrix::AbstractMatrix, ω::Number, cells = missing) +function greens!(matrix, g::GreensFunction{<:BandGreensSolver,L}, ω::Number, (src, dst)::SVectorPair{L}) where {L} fill!(matrix, zero(eltype(matrix))) - cells´ = sanitize_dn(cells, g.h) + dn = dst - src for simplexdata in g.solver.simplexdata - g0, gjs = green_simplex(ω, cells´, simplexdata, g.solver.indsedges) + g0, gjs = green_simplex(ω, dn, simplexdata, g.solver.indsedges) addsimplex!(matrix, g0, gjs, simplexdata) end return matrix end -sanitize_dn((src, dst)::Pair, ::Hamiltonian{LA,L}) where {LA,L} = - SVector{L}(dst) - SVector{L}(src) - -sanitize_dn(::Missing, ::Hamiltonian{LA,L}) where {LA,L} = - zero(SVector{L,Int}) - function green_simplex(ω, dn, data::SimplexData{L}, indsedges) where {L} dη = data.Δks' * dn phase = cis(dot(dn, data.k0)) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 6f6f9f46..cb2ebe3a 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -71,6 +71,8 @@ sublats(o::OrbitalStructure) = 1:length(o.orbitals) orbitals(o::OrbitalStructure) = o.orbitals +siterange(o::OrbitalStructure, sublat) = (1+o.offsets[sublat]):o.offsets[sublat+1] + # sublat offsets after flattening (without padding zeros) flatoffsets(offsets, norbs) = _flatoffsets((0,), offsets, norbs...) _flatoffsets(offsets´::NTuple{N,Any}, offsets, n, ns...) where {N} = @@ -82,7 +84,7 @@ _flatoffsets(offsets´, offsets) = offsets´ function flatoffsetorbs_site(i, orbstruct) s = sublat_site(i, orbstruct) - N = length(orbstruct.orbitals[s]) + N = length(orbitals(orbstruct)[s]) offset = orbstruct.offsets[s] offset´ = orbstruct.flatoffsets[s] Δi = i - offset @@ -199,7 +201,7 @@ _flatten(src::AbstractArray{T}, orbstruct, ::Type{T}) where {T<:Number} = src _flatten(src::AbstractArray{T}, orbstruct, ::Type{T´}) where {T<:Number,T´<:Number} = T´.(src) function _flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} - norbs = length.(orbstruct.orbitals) + norbs = length.(orbitals(orbstruct)) offsets´ = orbstruct.flatoffsets dim´ = last(offsets´) @@ -225,7 +227,7 @@ function _flatten(src::SparseMatrixCSC{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} end function _flatten(src::StridedMatrix{<:SMatrix{N,N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} - norbs = length.(orbstruct.orbitals) + norbs = length.(orbitals(orbstruct)) offsets´ = orbstruct.flatoffsets dim´ = last(offsets´) matrix = similar(src, T´, dim´, dim´) @@ -245,7 +247,7 @@ end # for Subspace bases function _flatten(src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T) where {N,T,T´} - norbs = length.(orbstruct.orbitals) + norbs = length.(orbitals(orbstruct)) offsets´ = orbstruct.flatoffsets dim´ = last(offsets´) matrix = similar(src, T´, dim´, size(src, 2)) @@ -263,15 +265,15 @@ function _flatten(src::StridedMatrix{<:SVector{N,T}}, orbstruct, ::Type{T´} = T end function flatten(lat::Lattice, orbstruct) - length(orbstruct.orbitals) == nsublats(lat) || throw(ArgumentError("Mismatch between sublattices and orbitals")) - unitcell´ = flatten(lat.unitcell, orbstruct) #length.(orbstruct.orbitals)) + length(orbitals(orbstruct)) == nsublats(lat) || throw(ArgumentError("Mismatch between sublattices and orbitals")) + unitcell´ = flatten(lat.unitcell, orbstruct) bravais´ = lat.bravais lat´ = Lattice(bravais´, unitcell´) end function flatten(unitcell::Unitcell, orbstruct) # orbs::NTuple{S,Any}) where {S} - norbs = length.(orbstruct.orbitals) - nsublats = length(orbstruct.orbitals) + norbs = length.(orbitals(orbstruct)) + nsublats = length(sublats(orbstruct)) offsets´ = orbstruct.flatoffsets ns´ = last(offsets´) sites´ = similar(unitcell.sites, ns´) @@ -285,17 +287,17 @@ function flatten(unitcell::Unitcell, orbstruct) # orbs::NTuple{S,Any}) where {S} return unitcell´ end -function flatten_sparse_copy!(dst, src, h) +function flatten_sparse_copy!(dst, src, o::OrbitalStructure) fill!(dst, zero(eltype(dst))) - norbs = length.(orbitals(h)) + norbs = length.(orbitals(o)) coloffset = 0 - for s´ in sublats(h.lattice) + for s´ in sublats(o) N´ = norbs[s´] - for col in siterange(h.lattice, s´) + for col in siterange(o, s´) for p in nzrange(src, col) val = nonzeros(src)[p] row = rowvals(src)[p] - rowoffset, M´ = flatoffsetorbs_site(row, h.orbstruct) + rowoffset, M´ = flatoffsetorbs_site(row, o) for j in 1:N´, i in 1:M´ dst[i + rowoffset, j + coloffset] = val[i, j] end @@ -306,16 +308,16 @@ function flatten_sparse_copy!(dst, src, h) return dst end -function flatten_sparse_muladd!(dst, src, h, α = I) - norbs = length.(orbitals(h)) +function flatten_sparse_muladd!(dst, src, o::OrbitalStructure, α = I) + norbs = length.(orbitals(o)) coloffset = 0 - for s´ in sublats(h.lattice) + for s´ in sublats(o) N´ = norbs[s´] - for col in siterange(h.lattice, s´) + for col in siterange(o, s´) for p in nzrange(src, col) val = α * nonzeros(src)[p] row = rowvals(src)[p] - rowoffset, M´ = flatoffsetorbs_site(row, h.orbstruct) + rowoffset, M´ = flatoffsetorbs_site(row, o) for j in 1:N´, i in 1:M´ dst[i + rowoffset, j + coloffset] += val[i, j] end @@ -326,16 +328,16 @@ function flatten_sparse_muladd!(dst, src, h, α = I) return dst end -function flatten_dense_muladd!(dst, src, h, α = I) - norbs = length.(orbitals(h)) +function flatten_dense_muladd!(dst, src, o::OrbitalStructure, α = I) + norbs = length.(orbitals(o)) coloffset = 0 - for s´ in sublats(h.lattice) + for s´ in sublats(o) N´ = norbs[s´] - for col in siterange(h.lattice, s´) + for col in siterange(o, s´) rowoffset = 0 - for s in sublats(h.lattice) + for s in sublats(o) M´ = norbs[s] - for row in siterange(h.lattice, s) + for row in siterange(o, s) val = α * src[row, col] for j in 1:N´, i in 1:M´ dst[i + rowoffset, j + coloffset] += val[i, j] @@ -349,9 +351,9 @@ function flatten_dense_muladd!(dst, src, h, α = I) return dst end -function flatten_dense_copy!(dst, src, h) +function flatten_dense_copy!(dst, src, o::OrbitalStructure) fill!(dst, zero(eltype(dst))) - return flatten_dense_muladd!(dst, src, h, I) + return flatten_dense_muladd!(dst, src, o, I) end ## unflatten ## @@ -396,7 +398,7 @@ unflatten(vflat::AbstractMatrix, o::OrbitalStructure{T}) where {T} = unflatten!(similar(vflat, T, dimh(o), size(vflat, 2)), vflat, o) function unflatten!(v::AbstractArray{T}, vflat::AbstractArray, o::OrbitalStructure) where {T} - norbs = length.(o.orbitals) + norbs = length.(orbitals(o)) flatoffsets = o.flatoffsets dimflat = last(flatoffsets) check_unflatten_dst_dims(v, dimh(o)) @@ -420,8 +422,14 @@ unflatten_or_reinterpret(vflat, ::Missing) = vflat # source is already of the correct orbitaltype(h) function unflatten_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{T}) where {T} check_unflatten_dst_dims(vflat, dimh(o)) - vflat + return vflat end + +function unflatten_or_reinterpret(vflat::AbstractArray{<:Number}, o::OrbitalStructure{<:Number}) + check_unflatten_dst_dims(vflat, dimh(o)) + return vflat +end + # source can be reinterpreted, because the number of orbitals is the same M for all N sublattices unflatten_or_reinterpret(vflat::AbstractArray{T}, o::OrbitalStructure{S,N,<:NTuple{N,NTuple{M}}}) where {T,S,N,M} = reinterpret(SVector{M,T}, vflat) @@ -1738,14 +1746,14 @@ bloch!(matrix, h::Hamiltonian, ϕs, axis = 0) = _bloch!(matrix, h, toSVector(ϕs bloch!(matrix, h::Hamiltonian, ϕs::Tuple{SVector,NamedTuple}, args...) = bloch!(matrix, h, first(ϕs), args...) function bloch!(matrix, h::Hamiltonian) - _copy!(parent(matrix), first(h.harmonics).h, h) # faster copy!(dense, sparse) specialization + _copy!(parent(matrix), first(h.harmonics).h, h.orbstruct) # faster copy!(dense, sparse) specialization return matrix end function _bloch!(matrix::AbstractMatrix, h::Hamiltonian{<:Lattice,L,M}, ϕs::SVector{L}, axis::Number) where {L,M} rawmatrix = parent(matrix) if iszero(axis) - _copy!(rawmatrix, first(h.harmonics).h, h) # faster copy!(dense, sparse) specialization + _copy!(rawmatrix, first(h.harmonics).h, h.orbstruct) # faster copy!(dense, sparse) specialization add_harmonics!(rawmatrix, h, ϕs, dn -> 1) else fill!(rawmatrix, zero(M)) # There is no guarantee of same structure @@ -1760,7 +1768,7 @@ function _bloch!(matrix::AbstractMatrix, h::Hamiltonian{<:Lattice,L,M}, ϕs::SVe if iszero(prefactor0) fill!(rawmatrix, zero(eltype(rawmatrix))) else - _copy!(rawmatrix, first(h.harmonics).h, h) + _copy!(rawmatrix, first(h.harmonics).h, h.orbstruct) rmul!(rawmatrix, prefactor0) end add_harmonics!(rawmatrix, h, ϕs, dnfunc) @@ -1779,24 +1787,24 @@ function add_harmonics!(zerobloch, h::Hamiltonian{<:Lattice,L}, ϕs::SVector{L}, prefactor = dnfunc(hh.dn) iszero(prefactor) && continue ephi = prefactor * cis(-ϕs´ * hh.dn) - _add!(zerobloch, hhmatrix, h, ephi) + _add!(zerobloch, hhmatrix, h.orbstruct, ephi) end return zerobloch end ############################################################################################ -######## _copy! and _add! call specialized methods in tools.jl ############################# +######## _copy! and _add! call specialized methods ######################################### ############################################################################################ _copy!(dest, src, h) = copy!(dest, src) -_copy!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, h) = _fast_sparse_copy!(dst, src) -_copy!(dst::StridedMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, h) = _fast_sparse_copy!(dst, src) -_copy!(dst::StridedMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, h) where {N} = _fast_sparse_copy!(dst, src) -_copy!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:SMatrix}, h) = flatten_sparse_copy!(dst, src, h) -_copy!(dst::StridedMatrix{<:Number}, src::StridedMatrix{<:SMatrix}, h) = flatten_dense_copy!(dst, src, h) - -_add!(dest, src, h, α) = _plain_muladd!(dest, src, α) -_add!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, h, α = 1) = _fast_sparse_muladd!(dst, src, α) -_add!(dst::AbstractMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, h, α = I) where {N} = _fast_sparse_muladd!(dst, src, α) -_add!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:SMatrix}, h, α = I) = flatten_sparse_muladd!(dst, src, h, α) -_add!(dst::StridedMatrix{<:Number}, src::StridedMatrix{<:SMatrix}, h, α = I) = flatten_dense_muladd!(dst, src, h, α) \ No newline at end of file +_copy!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, o) = _fast_sparse_copy!(dst, src) +_copy!(dst::StridedMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, o) = _fast_sparse_copy!(dst, src) +_copy!(dst::StridedMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, o) where {N} = _fast_sparse_copy!(dst, src) +_copy!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:SMatrix}, o) = flatten_sparse_copy!(dst, src, o) +_copy!(dst::StridedMatrix{<:Number}, src::StridedMatrix{<:SMatrix}, o) = flatten_dense_copy!(dst, src, o) + +_add!(dest, src, o, α) = _plain_muladd!(dest, src, α) +_add!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, o, α = 1) = _fast_sparse_muladd!(dst, src, α) +_add!(dst::AbstractMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, o, α = I) where {N} = _fast_sparse_muladd!(dst, src, α) +_add!(dst::AbstractMatrix{<:Number}, src::SparseMatrixCSC{<:SMatrix}, o, α = I) = flatten_sparse_muladd!(dst, src, o, α) +_add!(dst::StridedMatrix{<:Number}, src::StridedMatrix{<:SMatrix}, o, α = I) = flatten_dense_muladd!(dst, src, o, α) \ No newline at end of file diff --git a/src/iterators.jl b/src/iterators.jl index 0b3ca559..ff7dec24 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -370,8 +370,8 @@ struct Runs{T,F} istogether::F end -approxruns(xs::Vector{T}) where {T<:Number} = Runs(xs, (x, y) -> isapprox(x, y; atol = sqrt(eps(real(T))))) equalruns(xs) = Runs(xs, ==) +approxruns(xs::Vector{T}) where {T<:Number} = Runs(xs, (x, y) -> isapprox(x, y; atol = sqrt(eps(real(T))))) function last_in_run(xs::Vector{T}, i, istogether) where {T} xi = xs[i] diff --git a/src/plot_vegalite.jl b/src/plot_vegalite.jl index b4a74027..edfcf324 100644 --- a/src/plot_vegalite.jl +++ b/src/plot_vegalite.jl @@ -102,7 +102,16 @@ Plot the the Hamiltonian lattice projected along `axes` using VegaLite. Plot an eigenstate `psi` on the lattice, using one of various possible visualization shaders: `sitesize`, `siteopacity`, `sitecolor`, `linksize`, `linkopacity` or `linkcolor` (see keywords below). If `psi` is obtained in a flattened form, it will be automatically -"unflattened" to restore the orbital structure of `h`. +"unflattened" to restore the orbital structure of `h`. The vector `psi` must have length +equal to the number of sites (after unflattening if required). + +If `h` is defined on an unbounded `L`-dimensional lattice, and `psi` is of the form + + psi = [cell₁ => psi₁, cell₂ => psi₂,...] + +where `cellᵢ::NTuple{L,Int}` are indices of unit cells, said unit cells will all be plotted +together with their corresponding `psiᵢ`. Otherwise, a single cell at the origin will be +plotted. vlplot(h::Hamiltonian, psi::AbstractMatrix; kw...) @@ -130,6 +139,7 @@ Equivalent to `vlplot(h, psi.basis; kw...)`. - `digits = 4`: number of significant digits to show in onsite energy and hopping tooltips - `plotsites = true`: whether to plot sites - `plotlinks = true`: whether to plot links + - `onlyonecell = false`: whether to omit plotting links to neighboring unit cells in unbounded lattices - `sitestroke = :white`: the color of site outlines. If `nothing`, no outline will be plotted. - `colorscheme = "redyellowblue"`: Color scheme from `https://vega.github.io/vega/docs/schemes/` - `discretecolorscheme = "category10"`: Color scheme from `https://vega.github.io/vega/docs/schemes/` @@ -201,24 +211,16 @@ end function VegaLite.vlplot(h::Hamiltonian{LA}, psi = missing; labels = ("x","y"), size = 800, axes::Tuple{Int,Int} = (1,2), xlims = missing, ylims = missing, digits = 4, - plotsites = true, plotlinks = true, + onlyonecell = false, plotsites = true, plotlinks = true, maxdiameter = 20, mindiameter = 0, maxthickness = 0.25, sitestroke = :white, sitesize = maxdiameter, siteopacity = 0.9, sitecolor = missing, linksize = maxthickness, linkopacity = 1.0, linkcolor = missing, colorrange = missing, colorscheme = "redyellowblue", discretecolorscheme = "category10") where {E,LA<:Lattice{E}} - psi´ = unflatten_or_reinterpret_or_missing(psi, h) - checkdims_psi(h, psi´) - directives = (; axes = axes, digits = digits, - sitesize_shader = site_shader(sitesize, h, psi´), siteopacity_shader = site_shader(siteopacity, h, psi´), - sitecolor_shader = site_shader(sitecolor, h, psi´), linksize_shader = link_shader(linksize, h, psi´), - linkopacity_shader = link_shader(linkopacity, h, psi´), linkcolor_shader = link_shader(linkcolor, h, psi´)) - - table = linkstable(h, directives, plotsites, plotlinks) - maxsiteopacity = plotsites ? maximum(s.opacity for s in table if !s.islink) : 0.0 + table = linkstable(h, psi; axes, digits, onlyonecell, plotsites, plotlinks, + sitesize, siteopacity, sitecolor, linksize, linkopacity, linkcolor) maxsitesize = plotsites ? maximum(s.scale for s in table if !s.islink) : 0.0 - maxlinkopacity = plotlinks ? maximum(s.opacity for s in table if s.islink) : 0.0 maxlinksize = plotlinks ? maximum(s.scale for s in table if s.islink) : 0.0 mindiameter > 0 && filter!(s -> s.islink || s.scale >= mindiameter*maxsitesize/maxdiameter, table) @@ -248,7 +250,7 @@ function VegaLite.vlplot(h::Hamiltonian{LA}, psi = missing; scale = {domain = colorrange´, scheme = colorscheme´}, legend = needslegend(linkcolor)}, strokeOpacity = {:opacity, - scale = {range = (0, 1), domain = (0, maxlinkopacity)}, + scale = {range = (0, 1), domain = (0, 1)}, legend = needslegend(linkopacity)}, transform = [{filter = "datum.islink"}], selection = {grid2 = {type = :interval, bind = :scales}}, @@ -268,7 +270,7 @@ function VegaLite.vlplot(h::Hamiltonian{LA}, psi = missing; scale = {domain = colorrange´, scheme = colorscheme´}, legend = needslegend(sitecolor)}, opacity = {:opacity, - scale = {range = (0, 1), domain = (0, maxsiteopacity)}, + scale = {range = (0, 1), domain = (0, 1)}, legend = needslegend(siteopacity)}, selection = {grid1 = {type = :interval, bind = :scales}}, transform = [{filter = "!datum.islink"}], @@ -299,32 +301,58 @@ needslegend(x::Number) = nothing needslegend(x) = true unflatten_or_reinterpret_or_missing(psi::Missing, h) = missing -unflatten_or_reinterpret_or_missing(psi::AbstractArray, h) = unflatten_or_reinterpret(psi, h) - -unflatten_or_reinterpret_or_missing(s::Subspace, h) = unflatten_or_reinterpret(s.basis, h) +unflatten_or_reinterpret_or_missing(psi::AbstractArray, h) = unflatten_or_reinterpret(psi, h.orbstruct) +unflatten_or_reinterpret_or_missing(s::Subspace, h) = unflatten_or_reinterpret(s.basis, h.orbstruct) checkdims_psi(h, psi) = size(h, 2) == size(psi, 1) || throw(ArgumentError("The eigenstate length $(size(psi,1)) must match the Hamiltonian dimension $(size(h, 2))")) checkdims_psi(h, ::Missing) = nothing -function linkstable(h::Hamiltonian, d, plotsites, plotlinks) +function linkstable(h, psi; kw...) + T = numbertype(h.lattice) + table = NamedTuple{(:x, :y, :x2, :y2, :sublat, :tooltip, :scale, :color, :opacity, :islink), + Tuple{T,T,T,T,NameType,String,Float64,Float64,Float64,Bool}}[] + linkstable!(table, h, psi; kw...) + return table +end + +function linkstable!(table, h, cellpsi::AbstractVector{<:Pair}; kw...) + foreach(cellpsi) do (cell, psi) + linkstable!(table, h, psi, cell; kw..., onlyonecell = true) + end + return table +end + +function linkstable!(table, h, psi, cell = missing; + axes, digits, onlyonecell, plotsites, plotlinks, + sitesize, siteopacity, sitecolor, linksize, linkopacity, linkcolor) where {LA,L} + psi´ = unflatten_or_reinterpret_or_missing(psi, h) + checkdims_psi(h, psi´) + d = (; axes = axes, digits = digits, + sitesize_shader = site_shader(sitesize, h, psi´), + siteopacity_shader = site_shader(siteopacity, h, psi´), + sitecolor_shader = site_shader(sitecolor, h, psi´), + linksize_shader = link_shader(linksize, h, psi´), + linkopacity_shader = link_shader(linkopacity, h, psi´), + linkcolor_shader = link_shader(linkcolor, h, psi´)) (a1, a2) = d.axes lat = h.lattice + br = bravais(h) T = numbertype(lat) + r0 = cell === missing ? br * SVector(filltuple(0, latdim(h))) : br * SVector(cell) slats = sublats(lat) rs = allsitepositions(lat) - table = NamedTuple{(:x, :y, :x2, :y2, :sublat, :tooltip, :scale, :color, :opacity, :islink), - Tuple{T,T,T,T,NameType,String,Float64,Float64,Float64,Bool}}[] h0 = h.harmonics[1].h rows = Int[] # list of plotted destination sites for dn != 0 for har in h.harmonics + !iszero(har.dn) && onlyonecell && continue resize!(rows, 0) ridx = 0 for ssrc in slats if iszero(har.dn) for (i, r) in enumeratesites(lat, ssrc) - x = get(r, a1, zero(T)); y = get(r, a2, zero(T)) + x = get(r + r0, a1, zero(T)); y = get(r + r0, a2, zero(T)) ridx += 1 plotsites && push!(table, (x = x, y = y, x2 = x, y2 = y, sublat = sublatname(lat, ssrc), tooltip = matrixstring_inline(i, h0[i, i], d.digits), @@ -336,8 +364,8 @@ function linkstable(h::Hamiltonian, d, plotsites, plotlinks) itr = nonzero_indices(har, siterange(lat, sdst), siterange(lat, ssrc)) for (row, col) in itr iszero(har.dn) && col == row && continue - rsrc = rs[col] - rdst = (rs[row] + bravais(lat) * har.dn) + rsrc = rs[col] + r0 + rdst = (rs[row] + bravais(lat) * har.dn + r0) # sites in neighboring cells connected to dn=0 cell. But add only once. if !iszero(har.dn) && !in(row, rows) x = get(rdst, a1, zero(T)); y = get(rdst, a2, zero(T)) @@ -350,7 +378,7 @@ function linkstable(h::Hamiltonian, d, plotsites, plotlinks) end plotlinks || continue # draw half-links but only intracell - rdst = ifelse(iszero(har.dn), (rdst + rsrc) / 2, rdst) + rdst = ifelse(cell !== missing || iszero(har.dn), (rdst + rsrc) / 2, rdst) x = get(rsrc, a1, zero(T)); y = get(rsrc, a2, zero(T)) x´ = get(rdst, a1, zero(T)); y´ = get(rdst, a2, zero(T)) # Exclude links perpendicular to the screen diff --git a/src/tools.jl b/src/tools.jl index 9d0824ea..3d5c49c1 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -27,7 +27,7 @@ ensuretuple(s) = (s,) indstopair(s::Tuple) = Pair(last(s), first(s)) -filltuple(x, ::Val{L}) where {L} = ntuple(_ -> x, Val(L)) +filltuple(x, L) = ntuple(_ -> x, L) filltuple(x, ::NTuple{N,Any}) where {N} = ntuple(_ -> x, Val(N)) # # toSVector can deal with the L=0 edge case, unlike SVector @@ -135,6 +135,8 @@ displayvectors(mat::SMatrix{E,L,<:AbstractFloat}; kw...) where {E,L} = displayvectors(mat::SMatrix{E,L,<:Integer}; kw...) where {E,L} = ntuple(l -> Tuple(mat[:,l]), Val(L)) +chop(x::T) where {T} = ifelse(abs2(x) < eps(real(T)), zero(T), x) + # pseudoinverse of supercell s times an integer n, so that it is an integer matrix (for accuracy) pinvmultiple(s::SMatrix{L,0}) where {L} = (SMatrix{0,0,Int}(), 0) function pinvmultiple(s::SMatrix{L,L´}) where {L,L´} @@ -155,6 +157,8 @@ function pinverse(m::SMatrix) return inv(qrm.R) * qrm.Q' end +issquare(a::AbstractMatrix) = size(a, 1) == size(a, 2) + # normalize_axis_directions(q::SMatrix{M,N}) where {M,N} = hcat(ntuple(i->q[:,i]*sign(q[i,i]), Val(N))...) padprojector(::Type{S}, ::Val{N}) where {M,N,S<:SMatrix{M,M}} = S(Diagonal(SVector(padright(filltuple(1, Val(N)), Val(M))))) diff --git a/test/test_greens.jl b/test/test_greens.jl index 469e5993..f35d4e75 100644 --- a/test/test_greens.jl +++ b/test/test_greens.jl @@ -1,9 +1,44 @@ using LinearAlgebra: tr -@testset "basic green's function" begin +@testset "basic greens function" begin g = LatticePresets.square() |> hamiltonian(hopping(-1)) |> unitcell((2,0), (0,1)) |> greens(bandstructure(resolution = 17)) @test g(0.2) ≈ transpose(g(0.2)) @test imag(tr(g(1))) < 0 @test Quantica.chop(imag(tr(g(5)))) == 0 -end \ No newline at end of file +end + +@testset "greens functions spectra" begin + h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r->abs(r[2])<3) + g = greens(h, SingleShot1D()) + @test_throws ArgumentError g(0) + dos = [-imag(tr(g(w + 1.0e-6im))) for w in range(-3.2, 3.2, length = 1001)] + @test all(x -> Quantica.chop(x) >= 0, dos) + + h = LP.honeycomb() |> hamiltonian(hopping(1, range = 2)) |> unitcell((1,-1), region = r->abs(r[2])<3) + g = greens(h, SingleShot1D()) + dos = [-imag(tr(g(w))) for w in range(-6, 6, length = 1001)] + @test all(x -> Quantica.chop(x) >= 0, dos) + + h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(3,3)) |> wrap(2) + g = greens(h, SingleShot1D()) + dos = [-imag(tr(g(w))) for w in range(-3.2, 3.2, length = 1001)] + @test all(x -> Quantica.chop(x) >= 0, dos) + + h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((2,-2),(3,3)) |> unitcell(1, 1, indices = not(1)) |> wrap(2) + g = greens(h, SingleShot1D()) + @test_broken imag(tr(g(0.2))) < 0 + g = greens(h, SingleShot1D(direct = false)) + @test imag(tr(g(0.2))) < 0 + dos = [-imag(tr(g(w + 1.0e-6im))) for w in range(-3.2, 3.2, length = 1001)] + @test all(x -> Quantica.chop(x) >= 0, dos) +end + +@testset "greens functions spatial" begin + h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1),(3,3)) |> unitcell((1,0)) + g = greens(h, SingleShot1D(), boundaries = (0,)) + g0 = g(0.01, missing) + ldos = [-imag(tr(g0(n=>n))) for n in 1:200] + @test all(x -> Quantica.chop(x) >= 0, ldos) +end +