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/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..f4cfcbfd 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{T}) where {T<:Number} + check_unflatten_dst_dims(vflat, dimh(o)) + return T.(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´}