From 395502c190643deb06a19d734a543e55ca86e884 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Tue, 5 Jan 2021 15:39:37 +0100 Subject: [PATCH 01/11] wip wip --- src/greens.jl | 194 +++++++++++++++++++++++++++++++++++---------- src/hamiltonian.jl | 86 ++++++++++---------- src/iterators.jl | 2 +- 3 files changed, 198 insertions(+), 84 deletions(-) diff --git a/src/greens.jl b/src/greens.jl index cdfa5d95..8d1a6e30 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -1,35 +1,38 @@ ####################################################################### # 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()` (for 1D Hamiltonians, to use the single-shot generalized eigenvalue approach) - h |> greens(h -> solveobject(h)) +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. -Curried form equivalent to the above, giving `greens(h, solveobject(h))` (see -example below). + h |> greens(h -> solveobject(h), args...) - g([m,] ω, cells::Pair = missing) +Curried form equivalent to the above, giving `greens(h, solveobject(h), args...)`. -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(ω, 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 `cells` is missing, `src` and `dst` are +assumed to be zero vectors. See also `greens!` to use a preallocated matrix. # Examples @@ -49,16 +52,132 @@ julia> m = similarmatrix(g); g(m, 0.2) 6.663377810046025 - 24.472789025006396im ``` +# See also + `greens!` """ -greens(h, solver) = GreensFunction(h, greensolver(solver)) -greens(solver::Function) = h -> greens(h, solver(h)) +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...) + +# fallback +greensolver(s::AbstractGreensSolver) = s + +# call API +(g::GreensFunction)(ω::Number, cells = missing) = 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)) + +sanitize_cells(::Missing, ::GreensFunction{S,L}) where {S,L} = + zero(SVector{L,Int}) => zero(SVector{L,Int}) +sanitize_cells((cell0, cell1)::Pair{Integer,Integer}, ::GreensFunction{S,1}) where {S} = + (cell0,) => (cell1,) +sanitize_cells(cells::Pair{NTuple{L,Integer},NTuple{L,Integer}}, ::GreensFunction{S,L}) where {S,L} = + cells +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}`")) -# Needed to make similarmatrix work with GreensFunction -matrixtype(g::GreensFunction) = Matrix{eltype(g.h)} -Base.parent(g::GreensFunction) = g.h +const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} ####################################################################### -# BandGreenSolver +# SingleShot1DGreensSolver +####################################################################### +struct SingleShot1D end + +struct SingleShot1DGreensSolver{T,O<:OrbitalStructure} <: AbstractGreensSolver + A::Matrix{T} + B::Matrix{T} + Acopy::Matrix{T} + Bcopy::Matrix{T} + maxdn::Int + orbstruct::O + H0block::CartesianIndices{2,Tuple{UnitRange{Int64}, UnitRange{Int64}}} + H0diag::Vector{T} +end +#= +Precomputes A = [0 I 0...; 0 0 I; ...; -V₂ -V₁... ω-H₀ -V₋₁] and B = [I 0 0...; 0 I 0...;...; ... 0 0 V₋₂] +(form for two Harmonics {V₁,V₂}) +These define the eigenvalue problem A*φ = λ B*φ, with λ = exp(i kx a0) and φ = [φ₀, φ₁, φ₂...φₘ], +where φₙ = λ φₙ₋₁ = λⁿ φ₀, and m = max +Since we need at least half the eigenpairs, we use LAPACK, and hence dense matrices. +=# +function SingleShot1DGreensSolver(h::Hamiltonian) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) + maxdn = maximum(har -> abs(first(har.dn)), h.harmonics) + H = unitcell(h, (maxdn,)) + dimh = flatsize(H, 1) + T = complex(blockeltype(H)) + A = zeros(T, 2dimh, 2dimh) + B = zeros(T, 2dimh, 2dimh) + @show H + H0, V´, V = H[(0,)], H[(-1,)], H[(1,)] + orbs = H.orbstruct + block1, block2 = 1:dimh, dimh+1:2dimh + copy!(view(A, block1, block2), I(dimh)) + _add!(view(A, block2, block1), v, orbs, -1) + _add!(view(A, block2, block2), h0, orbs, -1) + copy!(view(B, block1, block1), I(dimh)) + _add!(view(B, block2, block2), v´, orbs, 1) + H0block = CartesianIndices((block2, block2)) + H0diag = [-A[i, j] for (i, j) in zip(block2, block2)] + return SingleShot1DGreensSolver(A, B, copy(A), copy(B), maxdn, orbs, H0block, H0diag) +end + + +function (gs::SingleShot1DGreensSolver{T})(ω) where {T} + A = copy!(gs.Acopy, gs.A) + B = copy!(gs.Bcopy, gs.B) + iη = im * sqrt(eps(real(T))) + for (row, col, h0) in zip(gs.rngrow, gs.rngcol, gs.H0diag) + A[row, col] = ω + iη - h0 + end + λs, χs = eigen(A, B; sortby = abs) # not eigen! because we need the ω-H0 block later + return select_retarded(λs, χs, gs) +end + +function select_retarded(λs, χs, gs) + dimH = length(gs.H0diag) + ret = 1:length(λs)÷2 + adv = length(λs)÷2 + 1:length(λs) + λR = view(λs, ret) + λA = view(λs, adv) + φR = view(χs, 1:dimh, ret) + φA = view(χs, 1:dimh, adv) + φλR = view(χs, dimh+1:2dimh, ret) + @show φR * Diagonal(λR) * inv(φR) + @show φA * Diagonal(inv.(λA)) * inv(φA) + iG0 = view(gs.Acopy, gs.H0block) + tmp = view(gs.Bcopy, 1:dimh, ret) + return λR, φR, φλR, iG0, tmp +end + +function greensolver(::SingleShot1D, h) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $L-dimensional Hamiltonian")) + return SingleShot1DGreensSolver(h) +end + +# SingleShot1DGreensSolver provides the pieces to compute `GV = φR λR φR⁻¹` and from there `G = G_00 = (ω0 - H0 - V´GV)⁻¹` within the first `unit supercell` of a semi-infinite chain +# In an infinite chain we have instead `G_NN = (ω0 - H0 - V'GV - VGV')⁻¹`, for any N, where `VGV'= (V'G'V)'`. Here `GV` and `G'V` involve the retarded and advanced G sectors, respectively +# The retarded/advanced sectors are classified by the eigenmode velocity (positive, negative) if they are propagating, or abs(λ) (< 0, > 0) if they are evanescent +# The velocity is given by vᵢ = im * φᵢ'(V'λᵢ-Vλᵢ')φᵢ / φᵢ'φᵢ +# Unit supercells different from the first. Semi-infinite G_N0 = G_{N-1,0}VG = (GV)^N G. +# Each unit is an maxdn × maxdn block matrix of the actual `g` we want +function greens!(matrix, g::GreensFunction{<:SingleShot1DGreensSolver,1}, ω, (src, dst)::SVectorPair{1}) + λR, φR, φλR, iG0, tmp = g.solver(ω) + N = mod(abs(first(dst - src)), g.maxdn) + if !iszero(N) + λR .^= N + φλR = rmul!(φλR, Diagonal(λR)) + end + iG0φR = mul!(tmp, iG0, φR) + G = (iG0φR' \ φλR')' + copy!(matrix, G) + return matrix +end + +####################################################################### +# BandGreensSolver ####################################################################### struct SimplexData{D,E,T,C<:SMatrix,DD,SA<:SubArray} ε0::T @@ -75,25 +194,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)) +" Matrix size : $(size(g.solver.h, 1)) × $(size(g.h, 2)) + Element type : $(displayelements(g.solver.h)) Band simplices : $(length(g.solver.simplexdata))") end -Base.summary(g::GreensFunction{<:BandGreenSolver}) = - "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.h))D bandstructure" +Base.summary(g::GreensFunction{<:BandGreensSolver}) = + "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.solver.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 +288,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 76e887b2..9df88bc2 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -66,6 +66,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} = @@ -77,7 +79,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 @@ -194,7 +196,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´) @@ -220,7 +222,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´) @@ -240,7 +242,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)) @@ -258,15 +260,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´) @@ -280,17 +282,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 @@ -301,16 +303,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 @@ -321,16 +323,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] @@ -344,9 +346,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 ## @@ -391,7 +393,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)) @@ -1740,7 +1742,7 @@ 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 @@ -1755,7 +1757,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) @@ -1774,24 +1776,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] From 6ae20a27eb3bce5238c842468ca4706ddeb17a5d Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 14 Jan 2021 19:26:19 +0100 Subject: [PATCH 02/11] draft, no tests nonsquare phiR problem seems to work refactor more refactor fix almost, but formally clunky works away from true bound states split off devdocs cleanup removed omega0 heuristics change cell default cleanup show method more show fix test --- docs/devdocs/singleshot.md | 80 +++++++ src/diagonalizer.jl | 25 ++- src/greens.jl | 427 +++++++++++++++++++++++++++++-------- src/hamiltonian.jl | 2 +- src/tools.jl | 2 + 5 files changed, 436 insertions(+), 100 deletions(-) create mode 100644 docs/devdocs/singleshot.md 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/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 8d1a6e30..b26ae643 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -27,12 +27,18 @@ origin. Curried form equivalent to the above, giving `greens(h, solveobject(h), args...)`. - g(ω, cells::Pair = missing) + g(ω, cells::Pair) 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. See also `greens!` to use a preallocated matrix. +dst` are `::NTuple{L,Int}` or `SVector{L,Int}`. If not provided, `cells` default to +`(1, 1, ...) => (1, 1, ...)`. + + 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 @@ -59,22 +65,22 @@ greens(h::Hamiltonian{<:Any,L}, solverobject; boundaries = filltuple(missing, Va GreensFunction(greensolver(solverobject, h), h, boundaries) greens(solver::Function, args...; kw...) = h -> greens(solver(h), h, args...; kw...) -# fallback +# solver fallback greensolver(s::AbstractGreensSolver) = s -# call API -(g::GreensFunction)(ω::Number, cells = missing) = greens!(similarmatrix(g), g, ω, cells) +# 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)) -sanitize_cells(::Missing, ::GreensFunction{S,L}) where {S,L} = - zero(SVector{L,Int}) => zero(SVector{L,Int}) -sanitize_cells((cell0, cell1)::Pair{Integer,Integer}, ::GreensFunction{S,1}) where {S} = - (cell0,) => (cell1,) -sanitize_cells(cells::Pair{NTuple{L,Integer},NTuple{L,Integer}}, ::GreensFunction{S,L}) where {S,L} = - cells +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}`")) @@ -85,95 +91,340 @@ const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} ####################################################################### struct SingleShot1D end -struct SingleShot1DGreensSolver{T,O<:OrbitalStructure} <: AbstractGreensSolver +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 A::Matrix{T} B::Matrix{T} - Acopy::Matrix{T} - Bcopy::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 - orbstruct::O - H0block::CartesianIndices{2,Tuple{UnitRange{Int64}, UnitRange{Int64}}} - H0diag::Vector{T} -end -#= -Precomputes A = [0 I 0...; 0 0 I; ...; -V₂ -V₁... ω-H₀ -V₋₁] and B = [I 0 0...; 0 I 0...;...; ... 0 0 V₋₂] -(form for two Harmonics {V₁,V₂}) -These define the eigenvalue problem A*φ = λ B*φ, with λ = exp(i kx a0) and φ = [φ₀, φ₁, φ₂...φₘ], -where φₙ = λ φₙ₋₁ = λⁿ φ₀, and m = max -Since we need at least half the eigenpairs, we use LAPACK, and hence dense matrices. -=# + temps::SingleShot1DTemporaries{T} +end + +function Base.show(io::IO, g::GreensFunction{<:SingleShot1DGreensSolver}) + print(io, summary(g), "\n", +" Matrix size : $(size(g.h, 1)) × $(size(g.h, 2)) + Element type : $(displayelements(g.h)) + Boundaries : $(g.boundaries) + Reduced dims : $(size(g.solver.Ps, 1))") +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(::SingleShot1D, h) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $L-dimensional Hamiltonian")) + return SingleShot1DGreensSolver(h) +end + +## Preparation + function SingleShot1DGreensSolver(h::Hamiltonian) latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) maxdn = maximum(har -> abs(first(har.dn)), h.harmonics) - H = unitcell(h, (maxdn,)) - dimh = flatsize(H, 1) + H = flatten(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, 2dimh, 2dimh) - B = zeros(T, 2dimh, 2dimh) - @show H - H0, V´, V = H[(0,)], H[(-1,)], H[(1,)] - orbs = H.orbstruct - block1, block2 = 1:dimh, dimh+1:2dimh - copy!(view(A, block1, block2), I(dimh)) - _add!(view(A, block2, block1), v, orbs, -1) - _add!(view(A, block2, block2), h0, orbs, -1) - copy!(view(B, block1, block1), I(dimh)) - _add!(view(B, block2, block2), v´, orbs, 1) - H0block = CartesianIndices((block2, block2)) - H0diag = [-A[i, j] for (i, j) in zip(block2, block2)] - return SingleShot1DGreensSolver(A, B, copy(A), copy(B), maxdn, orbs, H0block, H0diag) -end - - -function (gs::SingleShot1DGreensSolver{T})(ω) where {T} - A = copy!(gs.Acopy, gs.A) - B = copy!(gs.Bcopy, gs.B) - iη = im * sqrt(eps(real(T))) - for (row, col, h0) in zip(gs.rngrow, gs.rngcol, gs.H0diag) - A[row, col] = ω + iη - h0 + 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(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 + Ps = U'[1:dim_s, :] + Pb = U'[dim_s+1:end, :] + Ps´ = W'[1:dim_s, :] + Pb´ = W'[dim_s+1:end, :] + + Pb, Ps = filter_projectors(V, H0, Pb, Ps) + Pb´, Ps´ = filter_projectors(V´, H0, Pb´, Ps´) + + return Pb, Ps, Ps´ +end + +function filter_projectors(V::AbstractMatrix{T}, H0, Pb, Ps) where {T} + Σ1´ = Ps * V' * Pb' * pinv(Pb * H0 * Pb') * Pb * H0 * Ps' + SVD = svd(Matrix(Σ1´), 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´ + Qs = U'[1:dim_s´, :] + Qb = U'[dim_s´+1:end, :] + Pb_filtered = vcat(Pb, Qb * Ps) + Ps_filtered = Qs * Ps + return Pb_filtered, Ps_filtered +end + +## Solver execution + +function eigen_funcbarrier(A::AbstractMatrix{T}, B)::Tuple{Vector{T},Matrix{T}} where {T} + factB = lu!(B) + B⁻¹A = ldiv!(factB, A) + λs, φχs = eigen!(B⁻¹A; sortby = abs) + return λs, φχs +end + +function (gs::SingleShot1DGreensSolver)(ω) + A, B = single_shot_surface_matrices(gs, ω) + λs, φχs = eigen_funcbarrier(A, 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 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 - λs, χs = eigen(A, B; sortby = abs) # not eigen! because we need the ω-H0 block later - return select_retarded(λs, χs, gs) + + # B22 = -A21' + B22 .= .- A21' + + return A, B +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 select_retarded(λs, χs, gs) - dimH = length(gs.H0diag) - ret = 1:length(λs)÷2 - adv = length(λs)÷2 + 1:length(λs) +# 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!!(λs, copy!(p´, p)) + Base.permutecols!!(φ, copy!(p´, p)) + Base.permutecols!!(χ, copy!(p´, p)) + Base.permutecols!!(φs, copy!(p´, p)) + + dim_s = size(φ, 2) ÷ 2 + ret = 1:dim_s + adv = dim_s+1:2dim_s + λR = view(λs, ret) - λA = view(λs, adv) - φR = view(χs, 1:dimh, ret) - φA = view(χs, 1:dimh, adv) - φλR = view(χs, dimh+1:2dimh, ret) - @show φR * Diagonal(λR) * inv(φR) - @show φA * Diagonal(inv.(λA)) * inv(φA) - iG0 = view(gs.Acopy, gs.H0block) - tmp = view(gs.Bcopy, 1:dimh, ret) - return λR, φR, φλR, iG0, tmp + χ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 = rdiv!(copyto!(gs.temps.ss1, I), lu!(φRs)) + + iλA = view(λs, adv) + iλA .= (λ -> ifelse(isnan(λ), 0, 1/λ)).(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´ = rdiv!(copyto!(gs.temps.ss2, I), lu!(χAs´)) + + return λR, χR, iφRs, iλA, φA, iχAs´ end -function greensolver(::SingleShot1D, h) - latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $L-dimensional Hamiltonian")) - return SingleShot1DGreensSolver(h) +# 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 eachindex(vs) + abs2λ = abs2(λs[i]) + 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 vs +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 -# SingleShot1DGreensSolver provides the pieces to compute `GV = φR λR φR⁻¹` and from there `G = G_00 = (ω0 - H0 - V´GV)⁻¹` within the first `unit supercell` of a semi-infinite chain -# In an infinite chain we have instead `G_NN = (ω0 - H0 - V'GV - VGV')⁻¹`, for any N, where `VGV'= (V'G'V)'`. Here `GV` and `G'V` involve the retarded and advanced G sectors, respectively -# The retarded/advanced sectors are classified by the eigenmode velocity (positive, negative) if they are propagating, or abs(λ) (< 0, > 0) if they are evanescent -# The velocity is given by vᵢ = im * φᵢ'(V'λᵢ-Vλᵢ')φᵢ / φᵢ'φᵢ -# Unit supercells different from the first. Semi-infinite G_N0 = G_{N-1,0}VG = (GV)^N G. -# Each unit is an maxdn × maxdn block matrix of the actual `g` we want -function greens!(matrix, g::GreensFunction{<:SingleShot1DGreensSolver,1}, ω, (src, dst)::SVectorPair{1}) - λR, φR, φλR, iG0, tmp = g.solver(ω) - N = mod(abs(first(dst - src)), g.maxdn) - if !iszero(N) - λR .^= N - φλR = rmul!(φλR, Diagonal(λR)) +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 - iG0φR = mul!(tmp, iG0, φR) - G = (iG0φR' \ φλR')' - copy!(matrix, G) - return matrix + return G∞ +end + +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 ####################################################################### @@ -202,13 +453,13 @@ end function Base.show(io::IO, g::GreensFunction{<:BandGreensSolver}) print(io, summary(g), "\n", -" Matrix size : $(size(g.solver.h, 1)) × $(size(g.h, 2)) - Element type : $(displayelements(g.solver.h)) +" 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{<:BandGreensSolver}) = - "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.solver.h))D bandstructure" + "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.h))D bandstructure" function greensolver(b::Bandstructure{D}, h) where {D} indsedges = tuplepairs(Val(D)) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 0f7a766d..28a0051c 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -1740,7 +1740,7 @@ 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 diff --git a/src/tools.jl b/src/tools.jl index 9d0824ea..d9b1d1b7 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -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´} From 1c07d73fe3566023a59476d5bfb0da1eadeef1fb Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 14 Jan 2021 19:38:53 +0100 Subject: [PATCH 03/11] restore heuristic omega0 --- src/greens.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/greens.jl b/src/greens.jl index b26ae643..67667c3d 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -202,7 +202,10 @@ function bulk_surface_projectors(H0::AbstractMatrix{T}, V, V´) where {T} end function filter_projectors(V::AbstractMatrix{T}, H0, Pb, Ps) where {T} - Σ1´ = Ps * V' * Pb' * pinv(Pb * H0 * Pb') * Pb * H0 * Ps' + # Σ1´ = Ps * V' * Pb' * pinv(Pb * H0 * Pb') * Pb * H0 * Ps' + PbH0Pb´ = Pb * H0 * Pb' + ω0 = one(T) + eigmax(PbH0Pb´) + Σ1´ = Ps * V' * Pb' * ldiv!(lu!(ω0*I - PbH0Pb´), Pb * H0 * Ps') SVD = svd(Matrix(Σ1´), full = true) W, S, U = SVD.U, SVD.S, SVD.V dim_b´ = count(s -> iszero(chop(s)), S) From ba9aca4b0aa96d94e025c8485f55ee3ac3d89a7a Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Thu, 14 Jan 2021 19:42:50 +0100 Subject: [PATCH 04/11] fix --- src/greens.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/greens.jl b/src/greens.jl index 67667c3d..1919169b 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -158,7 +158,7 @@ Base.summary(g::GreensFunction{<:SingleShot1DGreensSolver}) = hasbulk(gs::SingleShot1DGreensSolver) = !iszero(size(gs.Pb, 1)) function greensolver(::SingleShot1D, h) - latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $L-dimensional Hamiltonian")) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) return SingleShot1DGreensSolver(h) end From 07356a10158e6087fa005774156ec33558198eaf Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Sat, 16 Jan 2021 17:58:25 +0100 Subject: [PATCH 05/11] linkstable! plus fixes fix fix --- docs/devdocs/singleshot.md | 80 +++++++ src/Quantica.jl | 2 +- src/bandstructure.jl | 2 +- src/diagonalizer.jl | 25 ++- src/greens.jl | 430 +++++++++++++++++++++++++++++-------- src/hamiltonian.jl | 10 +- src/plot_vegalite.jl | 76 ++++--- src/tools.jl | 4 +- 8 files changed, 501 insertions(+), 128 deletions(-) create mode 100644 docs/devdocs/singleshot.md 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 8d1a6e30..1919169b 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -27,12 +27,18 @@ origin. Curried form equivalent to the above, giving `greens(h, solveobject(h), args...)`. - g(ω, cells::Pair = missing) + g(ω, cells::Pair) 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. See also `greens!` to use a preallocated matrix. +dst` are `::NTuple{L,Int}` or `SVector{L,Int}`. If not provided, `cells` default to +`(1, 1, ...) => (1, 1, ...)`. + + 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 @@ -59,22 +65,22 @@ greens(h::Hamiltonian{<:Any,L}, solverobject; boundaries = filltuple(missing, Va GreensFunction(greensolver(solverobject, h), h, boundaries) greens(solver::Function, args...; kw...) = h -> greens(solver(h), h, args...; kw...) -# fallback +# solver fallback greensolver(s::AbstractGreensSolver) = s -# call API -(g::GreensFunction)(ω::Number, cells = missing) = greens!(similarmatrix(g), g, ω, cells) +# 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)) -sanitize_cells(::Missing, ::GreensFunction{S,L}) where {S,L} = - zero(SVector{L,Int}) => zero(SVector{L,Int}) -sanitize_cells((cell0, cell1)::Pair{Integer,Integer}, ::GreensFunction{S,1}) where {S} = - (cell0,) => (cell1,) -sanitize_cells(cells::Pair{NTuple{L,Integer},NTuple{L,Integer}}, ::GreensFunction{S,L}) where {S,L} = - cells +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}`")) @@ -85,95 +91,343 @@ const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} ####################################################################### struct SingleShot1D end -struct SingleShot1DGreensSolver{T,O<:OrbitalStructure} <: AbstractGreensSolver +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 A::Matrix{T} B::Matrix{T} - Acopy::Matrix{T} - Bcopy::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 - orbstruct::O - H0block::CartesianIndices{2,Tuple{UnitRange{Int64}, UnitRange{Int64}}} - H0diag::Vector{T} -end -#= -Precomputes A = [0 I 0...; 0 0 I; ...; -V₂ -V₁... ω-H₀ -V₋₁] and B = [I 0 0...; 0 I 0...;...; ... 0 0 V₋₂] -(form for two Harmonics {V₁,V₂}) -These define the eigenvalue problem A*φ = λ B*φ, with λ = exp(i kx a0) and φ = [φ₀, φ₁, φ₂...φₘ], -where φₙ = λ φₙ₋₁ = λⁿ φ₀, and m = max -Since we need at least half the eigenpairs, we use LAPACK, and hence dense matrices. -=# + temps::SingleShot1DTemporaries{T} +end + +function Base.show(io::IO, g::GreensFunction{<:SingleShot1DGreensSolver}) + print(io, summary(g), "\n", +" Matrix size : $(size(g.h, 1)) × $(size(g.h, 2)) + Element type : $(displayelements(g.h)) + Boundaries : $(g.boundaries) + Reduced dims : $(size(g.solver.Ps, 1))") +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(::SingleShot1D, h) + latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) + return SingleShot1DGreensSolver(h) +end + +## Preparation + function SingleShot1DGreensSolver(h::Hamiltonian) latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $(latdim(h))-dimensional Hamiltonian")) maxdn = maximum(har -> abs(first(har.dn)), h.harmonics) - H = unitcell(h, (maxdn,)) - dimh = flatsize(H, 1) + H = flatten(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, 2dimh, 2dimh) - B = zeros(T, 2dimh, 2dimh) - @show H - H0, V´, V = H[(0,)], H[(-1,)], H[(1,)] - orbs = H.orbstruct - block1, block2 = 1:dimh, dimh+1:2dimh - copy!(view(A, block1, block2), I(dimh)) - _add!(view(A, block2, block1), v, orbs, -1) - _add!(view(A, block2, block2), h0, orbs, -1) - copy!(view(B, block1, block1), I(dimh)) - _add!(view(B, block2, block2), v´, orbs, 1) - H0block = CartesianIndices((block2, block2)) - H0diag = [-A[i, j] for (i, j) in zip(block2, block2)] - return SingleShot1DGreensSolver(A, B, copy(A), copy(B), maxdn, orbs, H0block, H0diag) -end - - -function (gs::SingleShot1DGreensSolver{T})(ω) where {T} - A = copy!(gs.Acopy, gs.A) - B = copy!(gs.Bcopy, gs.B) - iη = im * sqrt(eps(real(T))) - for (row, col, h0) in zip(gs.rngrow, gs.rngcol, gs.H0diag) - A[row, col] = ω + iη - h0 + 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(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 + Ps = U'[1:dim_s, :] + Pb = U'[dim_s+1:end, :] + Ps´ = W'[1:dim_s, :] + Pb´ = W'[dim_s+1:end, :] + + Pb, Ps = filter_projectors(V, H0, Pb, Ps) + Pb´, Ps´ = filter_projectors(V´, H0, Pb´, Ps´) + + return Pb, Ps, Ps´ +end + +function filter_projectors(V::AbstractMatrix{T}, H0, Pb, Ps) where {T} + # Σ1´ = Ps * V' * Pb' * pinv(Pb * H0 * Pb') * Pb * H0 * Ps' + PbH0Pb´ = Pb * H0 * Pb' + ω0 = one(T) + eigmax(PbH0Pb´) + Σ1´ = Ps * V' * Pb' * ldiv!(lu!(ω0*I - PbH0Pb´), Pb * H0 * Ps') + SVD = svd(Matrix(Σ1´), 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´ + Qs = U'[1:dim_s´, :] + Qb = U'[dim_s´+1:end, :] + Pb_filtered = vcat(Pb, Qb * Ps) + Ps_filtered = Qs * Ps + return Pb_filtered, Ps_filtered +end + +## Solver execution + +function eigen_funcbarrier(A::AbstractMatrix{T}, B)::Tuple{Vector{T},Matrix{T}} where {T} + factB = lu!(B) + B⁻¹A = ldiv!(factB, A) + λs, φχs = eigen!(B⁻¹A; sortby = abs) + return λs, φχs +end + +function (gs::SingleShot1DGreensSolver)(ω) + A, B = single_shot_surface_matrices(gs, ω) + λs, φχs = eigen_funcbarrier(A, 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 - λs, χs = eigen(A, B; sortby = abs) # not eigen! because we need the ω-H0 block later - return select_retarded(λs, χs, gs) + return classify_retarded_advanced(λs, φs, χs, φ, χ, gs) end -function select_retarded(λs, χs, gs) - dimH = length(gs.H0diag) - ret = 1:length(λs)÷2 - adv = length(λs)÷2 + 1:length(λs) +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' + + return A, B +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!!(λs, copy!(p´, p)) + Base.permutecols!!(φ, copy!(p´, p)) + Base.permutecols!!(χ, copy!(p´, p)) + Base.permutecols!!(φs, copy!(p´, p)) + + dim_s = size(φ, 2) ÷ 2 + ret = 1:dim_s + adv = dim_s+1:2dim_s + λR = view(λs, ret) - λA = view(λs, adv) - φR = view(χs, 1:dimh, ret) - φA = view(χs, 1:dimh, adv) - φλR = view(χs, dimh+1:2dimh, ret) - @show φR * Diagonal(λR) * inv(φR) - @show φA * Diagonal(inv.(λA)) * inv(φA) - iG0 = view(gs.Acopy, gs.H0block) - tmp = view(gs.Bcopy, 1:dimh, ret) - return λR, φR, φλR, iG0, tmp + χ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 = rdiv!(copyto!(gs.temps.ss1, I), lu!(φRs)) + + iλA = view(λs, adv) + iλA .= (λ -> ifelse(isnan(λ), 0, 1/λ)).(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´ = rdiv!(copyto!(gs.temps.ss2, I), lu!(χAs´)) + + return λR, χR, iφRs, iλA, φA, iχAs´ end -function greensolver(::SingleShot1D, h) - latdim(h) == 1 || throw(ArgumentError("Cannot use a SingleShot1D Green function solver with an $L-dimensional Hamiltonian")) - return SingleShot1DGreensSolver(h) +# 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 eachindex(vs) + abs2λ = abs2(λs[i]) + 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 vs end -# SingleShot1DGreensSolver provides the pieces to compute `GV = φR λR φR⁻¹` and from there `G = G_00 = (ω0 - H0 - V´GV)⁻¹` within the first `unit supercell` of a semi-infinite chain -# In an infinite chain we have instead `G_NN = (ω0 - H0 - V'GV - VGV')⁻¹`, for any N, where `VGV'= (V'G'V)'`. Here `GV` and `G'V` involve the retarded and advanced G sectors, respectively -# The retarded/advanced sectors are classified by the eigenmode velocity (positive, negative) if they are propagating, or abs(λ) (< 0, > 0) if they are evanescent -# The velocity is given by vᵢ = im * φᵢ'(V'λᵢ-Vλᵢ')φᵢ / φᵢ'φᵢ -# Unit supercells different from the first. Semi-infinite G_N0 = G_{N-1,0}VG = (GV)^N G. -# Each unit is an maxdn × maxdn block matrix of the actual `g` we want -function greens!(matrix, g::GreensFunction{<:SingleShot1DGreensSolver,1}, ω, (src, dst)::SVectorPair{1}) - λR, φR, φλR, iG0, tmp = g.solver(ω) - N = mod(abs(first(dst - src)), g.maxdn) - if !iszero(N) - λR .^= N - φλR = rmul!(φλR, Diagonal(λR)) +## 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 - iG0φR = mul!(tmp, iG0, φR) - G = (iG0φR' \ φλR')' - copy!(matrix, G) - return matrix + return G∞ +end + +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 ####################################################################### @@ -202,13 +456,13 @@ end function Base.show(io::IO, g::GreensFunction{<:BandGreensSolver}) print(io, summary(g), "\n", -" Matrix size : $(size(g.solver.h, 1)) × $(size(g.h, 2)) - Element type : $(displayelements(g.solver.h)) +" 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{<:BandGreensSolver}) = - "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.solver.h))D bandstructure" + "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.h))D bandstructure" function greensolver(b::Bandstructure{D}, h) where {D} indsedges = tuplepairs(Val(D)) diff --git a/src/hamiltonian.jl b/src/hamiltonian.jl index 0f7a766d..cb2ebe3a 100644 --- a/src/hamiltonian.jl +++ b/src/hamiltonian.jl @@ -422,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) @@ -1740,7 +1746,7 @@ 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 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..c5b2654a 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´} From 3565e2888b8a0849f574f7b8c1658c0348ed4a91 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Wed, 20 Jan 2021 13:10:51 +0100 Subject: [PATCH 06/11] orbstruct bug --- src/plot_vegalite.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/plot_vegalite.jl b/src/plot_vegalite.jl index b4a74027..b30c14dc 100644 --- a/src/plot_vegalite.jl +++ b/src/plot_vegalite.jl @@ -299,9 +299,8 @@ 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))")) From e2bc3504b815525e1af06be07a86879d13e7c09a Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Wed, 20 Jan 2021 18:46:09 +0100 Subject: [PATCH 07/11] pseudoinverse for singular selfenergies --- src/greens.jl | 49 ++++++++++++++++++++----------------------------- src/tools.jl | 2 ++ 2 files changed, 22 insertions(+), 29 deletions(-) diff --git a/src/greens.jl b/src/greens.jl index 1919169b..6c6a549f 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -194,38 +194,22 @@ function bulk_surface_projectors(H0::AbstractMatrix{T}, V, V´) where {T} Pb = U'[dim_s+1:end, :] Ps´ = W'[1:dim_s, :] Pb´ = W'[dim_s+1:end, :] - - Pb, Ps = filter_projectors(V, H0, Pb, Ps) - Pb´, Ps´ = filter_projectors(V´, H0, Pb´, Ps´) - return Pb, Ps, Ps´ end -function filter_projectors(V::AbstractMatrix{T}, H0, Pb, Ps) where {T} - # Σ1´ = Ps * V' * Pb' * pinv(Pb * H0 * Pb') * Pb * H0 * Ps' - PbH0Pb´ = Pb * H0 * Pb' - ω0 = one(T) + eigmax(PbH0Pb´) - Σ1´ = Ps * V' * Pb' * ldiv!(lu!(ω0*I - PbH0Pb´), Pb * H0 * Ps') - SVD = svd(Matrix(Σ1´), 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´ - Qs = U'[1:dim_s´, :] - Qb = U'[dim_s´+1:end, :] - Pb_filtered = vcat(Pb, Qb * Ps) - Ps_filtered = Qs * Ps - return Pb_filtered, Ps_filtered -end - ## Solver execution function eigen_funcbarrier(A::AbstractMatrix{T}, B)::Tuple{Vector{T},Matrix{T}} where {T} - factB = lu!(B) - B⁻¹A = ldiv!(factB, A) - λs, φχs = eigen!(B⁻¹A; sortby = abs) + λs, φχs = eigen!(A, B; sortby = abs) + 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 (gs::SingleShot1DGreensSolver)(ω) A, B = single_shot_surface_matrices(gs, ω) λs, φχs = eigen_funcbarrier(A, B) @@ -307,24 +291,22 @@ function classify_retarded_advanced(λs, φs, χs, φ, χ, gs) Base.permutecols!!(χ, copy!(p´, p)) Base.permutecols!!(φs, copy!(p´, p)) - dim_s = size(φ, 2) ÷ 2 - ret = 1:dim_s - adv = dim_s+1:2dim_s + ret, adv = nonsingular_ret_adv(λs) λ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 = rdiv!(copyto!(gs.temps.ss1, I), lu!(φRs)) + iφRs = issquare(φRs) ? rdiv!(copyto!(gs.temps.ss1, I), lu!(φRs)) : pinv(φRs) iλA = view(λs, adv) - iλA .= (λ -> ifelse(isnan(λ), 0, 1/λ)).(iλA) + 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´ = rdiv!(copyto!(gs.temps.ss2, I), lu!(χAs´)) + 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 @@ -351,6 +333,15 @@ function compute_velocities_and_normalize!(φ, χ, λs, gs) return vs end +function nonsingular_ret_adv(λs::AbstractVector) + dim_s = length(λs) ÷ 2 + i0 = findfirst(!iszero, λs) + i0 === nothing && throw(ArgumentError("Unexpected non-conjugate λ solutions")) + ret = i0:dim_s + adv = dim_s+1:2dim_s+1-i0 + return ret, adv +end + ## Greens execution (g::GreensFunction{<:SingleShot1DGreensSolver})(ω, cells) = g(ω, missing)(sanitize_cells(cells, g)) diff --git a/src/tools.jl b/src/tools.jl index d9b1d1b7..1a51805b 100644 --- a/src/tools.jl +++ b/src/tools.jl @@ -157,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))))) From 05734a1f902cd2005b929e5c2f16821ab298c804 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 22 Jan 2021 12:52:51 +0100 Subject: [PATCH 08/11] direct kwarg --- src/greens.jl | 95 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 77 insertions(+), 18 deletions(-) diff --git a/src/greens.jl b/src/greens.jl index 6c6a549f..f93e8777 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -17,7 +17,7 @@ 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()` (for 1D Hamiltonians, to use the single-shot generalized eigenvalue approach) +- `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 @@ -89,7 +89,56 @@ const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} ####################################################################### # SingleShot1DGreensSolver ####################################################################### -struct SingleShot1D end +""" + SingleShot1D(; direct = false) + +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)`. If that happens one can try the +option `direct = true` that turns the generalized eigenproblem into a standard eigenproblem. + +# 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` +""" +struct SingleShot1D + invert_B::Bool +end + +SingleShot1D(; direct = false) = SingleShot1D(direct) struct SingleShot1DTemporaries{T} b2s::Matrix{T} # size = dim_b, 2dim_s @@ -127,6 +176,7 @@ function SingleShot1DTemporaries{T}(dim_H, dim_s, dim_b) where {T} end struct SingleShot1DGreensSolver{T<:Complex,R<:Real,H<:Hessenberg{T}} <: AbstractGreensSolver + invert_B::Bool A::Matrix{T} B::Matrix{T} minusH0::SparseMatrixCSC{T,Int} @@ -157,14 +207,14 @@ Base.summary(g::GreensFunction{<:SingleShot1DGreensSolver}) = hasbulk(gs::SingleShot1DGreensSolver) = !iszero(size(gs.Pb, 1)) -function greensolver(::SingleShot1D, h) +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) + return SingleShot1DGreensSolver(h, s.invert_B) end ## Preparation -function SingleShot1DGreensSolver(h::Hamiltonian) +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 = maximum(har -> abs(first(har.dn)), h.harmonics) H = flatten(unitcell(h, (maxdn,))) @@ -182,7 +232,7 @@ function SingleShot1DGreensSolver(h::Hamiltonian) 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(A, B, -H0, V, Pb, Ps, Ps´, H0ss, H0bs, Vss, Vbs, hessbb, velocities, maxdn, temps) + 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} @@ -199,20 +249,9 @@ end ## Solver execution -function eigen_funcbarrier(A::AbstractMatrix{T}, B)::Tuple{Vector{T},Matrix{T}} where {T} - λs, φχs = eigen!(A, B; sortby = abs) - 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 (gs::SingleShot1DGreensSolver)(ω) A, B = single_shot_surface_matrices(gs, ω) - λs, φχs = eigen_funcbarrier(A, B) + λ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, :) @@ -228,6 +267,23 @@ function (gs::SingleShot1DGreensSolver)(ω) 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) @@ -265,6 +321,9 @@ function single_shot_surface_matrices(gs::SingleShot1DGreensSolver{T}, ω) where # B22 = -A21' B22 .= .- A21' + A .= chop.(A) + B .= chop.(B) + return A, B end From b0ee16398d5d337e16dd53e77ab6d68151d32beb Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 22 Jan 2021 13:08:16 +0100 Subject: [PATCH 09/11] docstring [ci skip] docstring [ci skip] --- src/greens.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/greens.jl b/src/greens.jl index f93e8777..5e325982 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -90,7 +90,7 @@ const SVectorPair{L} = Pair{SVector{L,Int},SVector{L,Int}} # SingleShot1DGreensSolver ####################################################################### """ - SingleShot1D(; direct = false) + 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 @@ -111,8 +111,12 @@ 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)`. If that happens one can try the -option `direct = true` that turns the generalized eigenproblem into a standard eigenproblem. +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 @@ -138,7 +142,7 @@ struct SingleShot1D invert_B::Bool end -SingleShot1D(; direct = false) = SingleShot1D(direct) +SingleShot1D(; direct = true) = SingleShot1D(direct) struct SingleShot1DTemporaries{T} b2s::Matrix{T} # size = dim_b, 2dim_s From 447e67ce66a014bfc545dc2744bb54467b921c80 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 22 Jan 2021 18:47:50 +0100 Subject: [PATCH 10/11] nonsingular_ret_adv fix + tests --- src/Quantica.jl | 2 +- src/greens.jl | 76 +++++++++++++++++++++++++++++++-------------- test/test_greens.jl | 39 +++++++++++++++++++++-- 3 files changed, 90 insertions(+), 27 deletions(-) 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/greens.jl b/src/greens.jl index 5e325982..4b302dad 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -200,10 +200,10 @@ end function Base.show(io::IO, g::GreensFunction{<:SingleShot1DGreensSolver}) print(io, summary(g), "\n", -" Matrix size : $(size(g.h, 1)) × $(size(g.h, 2)) +" 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) - Reduced dims : $(size(g.solver.Ps, 1))") + Boundaries : $(g.boundaries)") end Base.summary(g::GreensFunction{<:SingleShot1DGreensSolver}) = @@ -220,8 +220,8 @@ end 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 = maximum(har -> abs(first(har.dn)), h.harmonics) - H = flatten(unitcell(h, (maxdn,))) + 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' @@ -244,10 +244,16 @@ function bulk_surface_projectors(H0::AbstractMatrix{T}, V, V´) where {T} W, S, U = SVD.U, SVD.S, SVD.V dim_b = count(s -> iszero(chop(s)), S) dim_s = length(S) - dim_b - Ps = U'[1:dim_s, :] - Pb = U'[dim_s+1:end, :] - Ps´ = W'[1:dim_s, :] - Pb´ = W'[dim_s+1:end, :] + 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 @@ -325,12 +331,21 @@ function single_shot_surface_matrices(gs::SingleShot1DGreensSolver{T}, ω) where # B22 = -A21' B22 .= .- A21' - A .= chop.(A) - B .= chop.(B) + 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 @@ -349,12 +364,16 @@ function classify_retarded_advanced(λs, φs, χ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) + 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.χ @@ -378,8 +397,8 @@ end function compute_velocities_and_normalize!(φ, χ, λs, gs) vs = gs.velocities tmp = gs.temps.vH - for i in eachindex(vs) - abs2λ = abs2(λs[i]) + for (i, λ) in enumerate(λs) + abs2λ = abs2(λ) if abs2λ ≈ 1 φi = view(φ, :, i) χi = view(χ, :, i) @@ -393,16 +412,25 @@ function compute_velocities_and_normalize!(φ, χ, λs, gs) 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 vs -end - -function nonsingular_ret_adv(λs::AbstractVector) - dim_s = length(λs) ÷ 2 - i0 = findfirst(!iszero, λs) - i0 === nothing && throw(ArgumentError("Unexpected non-conjugate λ solutions")) - ret = i0:dim_s - adv = dim_s+1:2dim_s+1-i0 - return ret, adv + 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 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 + From 3dfca42bf0fe5bd669f7f4de8427d5ca6f636980 Mon Sep 17 00:00:00 2001 From: Pablo San-Jose Date: Fri, 22 Jan 2021 20:34:41 +0100 Subject: [PATCH 11/11] docstring [ci skip] --- src/greens.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/greens.jl b/src/greens.jl index 4b302dad..be9e108d 100644 --- a/src/greens.jl +++ b/src/greens.jl @@ -59,7 +59,7 @@ julia> m = similarmatrix(g); g(m, 0.2) ``` # See also - `greens!` + `greens!`, `SingleShot1D` """ greens(h::Hamiltonian{<:Any,L}, solverobject; boundaries = filltuple(missing, Val(L))) where {L} = GreensFunction(greensolver(solverobject, h), h, boundaries)