Skip to content

Commit

Permalink
bandrangeKPM(::Hamiltonian)
Browse files Browse the repository at this point in the history 
dosKPM(ham) and densityKPM(ham)

typo

rebase
  • Loading branch information
@pablosanjose
pablosanjose committed Jun 17, 2020
1 parent 1f85af0 commit 38c69dd
Show file tree
Hide file tree
Showing 5 changed files with 173 additions and 127 deletions.
270 changes: 143 additions & 127 deletions src/KPM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,17 @@
# Kernel Polynomial Method : momenta
#######################################################################
struct MomentaKPM{T,B<:Tuple}
μlist::Vector{T}
mulist::Vector{T}
bandbracket::B
end

struct KPMBuilder{A,H<:AbstractMatrix,T,K,B}
A::A
h::H
bandbracket::Tuple{B,B}
order::Int
missingket::Bool
μlist::Vector{T}
mulist::Vector{T}
ket::K
ket0::K
ket1::K
Expand All @@ -20,13 +21,12 @@ end
function KPMBuilder(h::AbstractMatrix{Tv}, A = _defaultA(Tv); ket = missing, order = 10, bandrange = missing) where {Tv}
eh = eltype(eltype(h))
aA = eltype(eltype(A))
μlist = zeros(promote_type(eh, aA), order + 1)
mulist = zeros(promote_type(eh, aA), order + 1)
bandbracket = bandbracketKPM(h, bandrange)
missingket = ket === missing
ket´ = missingket ? ketundef(h) : ket
iscompatibleket(h, ket´) || throw(ArgumentError("ket is incompatible with Hamiltonian"))
builder = KPMBuilder(A, h, bandbracket, order, missingket,
μlist, ket´, similar(ket´), similar(ket´))
builder = KPMBuilder(A, h, bandbracket, order, missingket, mulist, ket´, similar(ket´), similar(ket´))
return builder
end

Expand All @@ -45,10 +45,11 @@ _iscompatibleket(::Type{S1}, ::Type{S2}) where {N, S1<:SMatrix{N,N}, S2<:SMatrix
_iscompatibleket(eltype(S1), eltype(S2))
_iscompatibleket(t1, t2) = false

function matrixKPM(h::Hamiltonian{<:Lattice,L}) where {L}
function matrixKPM(h::Hamiltonian{<:Lattice,L}, method = missing) where {L}
iszero(L) ||
throw(ArgumentError("Hamiltonian is defined on an infinite lattice. Convert it to a matrix first using `bloch(h, φs...)`"))
return bloch(h)
m = similarmatrix(h, method)
return bloch!(m, h)
end

matrixKPM(h::AbstractMatrix) = h
Expand All @@ -75,114 +76,133 @@ julia> momentaKPM(bloch(h), bandrange = (-6,6))
Quantica.MomentaKPM{Float64}([0.9594929736144973, -0.005881595972403821, -0.4933354572913581, 0.00359537502632597, 0.09759451291347333, -0.0008081453185250322, -0.00896262538765363, 0.00048205637037715177, -0.0003705198310034668, 9.64901673962623e-20, 9.110915988898614e-18], (0.0, 6.030150753768845))
```
"""
momentaKPM(h::Hamiltonian, A = _defaultA(eltype(h)); kw...) = momentaKPM(matrixKPM(h), matrixKPM(A); kw...)

function momentaKPM(h::AbstractMatrix, A = _defaultA(eltype(h)); randomkets = 1, kw...)
b = KPMBuilder(h, A; kw...)
if b.missingket
pmeter = Progress(b.order * randomkets, "Averaging moments: ")
for n in 1:randomkets
randomize!(b.ket)
addmomentaKPM!(b, pmeter)
end
b.μlist ./= randomkets
else
pmeter = Progress(b.order, "Computing moments: ")
addmomentaKPM!(b, pmeter)
end
return MomentaKPM(jackson!(b.μlist), b.bandbracket)
function momentaKPM(h::Hamiltonian, A = _defaultA(eltype(h)); bandrange = missing, kw...)
bandrange´ = bandrange === missing ? bandrangeKPM(h) : bandrange
momenta = momentaKPM(matrixKPM(h), matrixKPM(A); bandrange = bandrange´, kw...)
return momenta
end

_defaultA(::Type{T}) where {T<:Number} = one(T) * I
_defaultA(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = one(T) * I

# This iterates bras <psi_n| = <psi_0|AT_n(h) instead of kets (faster CSC multiplication)
# In practice we iterate their conjugate |psi_n> = T_n(h') A'|psi_0>, and do the projection
# onto the start ket, |psi_0>
function addmomentaKPM!(b::KPMBuilder{<:AbstractMatrix,<:AbstractSparseMatrix}, pmeter)
μlist, ket, ket0, ket1 = b.μlist, b.ket, b.ket0, b.ket1
h´, A´, bandbracket = b.h', b.A', b.bandbracket
order = length(μlist) - 1
fill!(ket0, zero(eltype(ket0)))
mul!(ket1, A´, ket)
μlist[1] += proj(ket1, ket)
ProgressMeter.next!(pmeter; showvalues = ())
for n in 2:(order+1)
ProgressMeter.next!(pmeter; showvalues = ())
μ = iterateKPM!(ket0, ket1, ket, h´, bandbracket)
μlist[n] += μ
n + 1 > order + 1 && break
ket0, ket1 = ket1, ket0
end
return μlist
end

function iterateKPM!(ket0::A, ket1::A, kini::A, adjh::Adjoint, (center, halfwidth)) where {S,A<:AbstractArray{S}}
h = adjh.parent
nzv = nonzeros(h)
rv = rowvals(h)
T = eltype(S)
μ = zero(T)
α = T(-2 * center / halfwidth)
β = T(2 / halfwidth)
for k in 1:size(ket0, 2)
for col in 1:size(h, 2)
@inbounds begin
tmp = α * ket1[col, k] - ket0[col, k]
for ptr in nzrange(h, col)
tmp += β * adjoint(nzv[ptr]) * ket1[rv[ptr],k]
end
ket0[col, k] = tmp
μ += dot(tmp, kini[col, k])
end
end
end
return μ
end

function addmomentaKPM!(b::KPMBuilder{<:UniformScaling, <:AbstractSparseMatrix}, pmeter)
μlist, ket, ket0, ket1 = b.μlist, b.ket, b.ket0, b.ket1
h´, A, bandbracket = b.h', b.A, b.bandbracket
order = length(μlist) - 1
fill!(ket0, zero(eltype(ket0)))
ket1 .= ket
for n in 1:2:(order+1)
ProgressMeter.next!(pmeter; showvalues = ())
ProgressMeter.next!(pmeter; showvalues = ()) # twice because of 2-step
μ, μ´ = iterateKPM!(ket0, ket1, h´, bandbracket)
μlist[n] += μ
n + 1 > order + 1 && break
μlist[n + 1] += μ´
ket0, ket1 = ket1, ket0
end
A.λ ≈ 1 || (μlist .*= A.λ)
return μlist
end

function iterateKPM!(ket0::A, ket1::A, adjh::Adjoint, (center, halfwidth)) where {S,A<:AbstractArray{S}}
h = adjh.parent
nzv = nonzeros(h)
rv = rowvals(h)
T = eltype(S)
μ = μ´ = zero(T)
α = T(-2 * center / halfwidth)
β = T(2 / halfwidth)
for k in 1:size(ket0, 2)
for col in 1:size(h, 2)
@inbounds begin
k1 = ket1[col, k]
tmp = α * k1 - ket0[col, k]
for ptr in nzrange(h, col)
tmp += β * adjoint(nzv[ptr]) * ket1[rv[ptr],k]
end
ket0[col, k] = tmp
μ += dot(k1, k1)
μ´ += dot(tmp, k1)
end
end
end
return μ, μ´
end
function momentaKPM(h::AbstractMatrix, A = _defaultA(eltype(h)); randomkets = 1, kw...)
b = KPMBuilder(h, A; kw...)
if b.missingket
pmeter = Progress(b.order * randomkets, "Averaging moments: ")
for n in 1:randomkets
randomize!(b.ket)
addmomentaKPM!(b, pmeter)
end
b.mulist ./= randomkets
else
pmeter = Progress(b.order, "Computing moments: ")
addmomentaKPM!(b, pmeter)
end
jackson!(b.mulist)
return MomentaKPM(b.mulist, b.bandbracket)
end

_defaultA(::Type{T}) where {T<:Number} = one(T) * I
_defaultA(::Type{S}) where {N,T,S<:SMatrix{N,N,T}} = one(T) * I

# This iterates bras <psi_n| = <psi_0|AT_n(h) instead of kets (faster CSC multiplication)
# In practice we iterate their conjugate |psi_n> = T_n(h') A'|psi_0>, and do the projection
# onto the start ket, |psi_0>
function addmomentaKPM!(b::KPMBuilder{<:AbstractMatrix,<:AbstractSparseMatrix}, pmeter)
mulist, ket, ket0, ket1 = b.mulist, b.ket, b.ket0, b.ket1
h, A, bandbracket = b.h, b.A, b.bandbracket
order = length(mulist) - 1
mul!(ket0, A', ket)
mulscaled!(ket1, h', ket0, bandbracket)
mulist[1] += proj(ket0, ket)
mulist[2] += proj(ket1, ket)
for n in 3:(order+1)
ProgressMeter.next!(pmeter; showvalues = ())
iterateKPM!(ket0, h', ket1, bandbracket)
mulist[n] += proj(ket0, ket)
ket0, ket1 = ket1, ket0
end
return mulist
end

function addmomentaKPM!(b::KPMBuilder{<:UniformScaling, <:AbstractSparseMatrix,T}, pmeter) where {T}
mulist, ket, ket0, ket1, = b.mulist, b.ket, b.ket0, b.ket1
h, A, bandbracket = b.h, b.A, b.bandbracket
order = length(mulist) - 1
ket0 .= ket
mulscaled!(ket1, h', ket0, bandbracket)
mulist[1] += μ0 = 1.0
mulist[2] += μ1 = proj(ket1, ket0)
# This is not used in the currently activated codepath (BLAS mul!), but is needed in the
# commented out @threads codepath
thread_buffers = (zeros(T, Threads.nthreads()), zeros(T, Threads.nthreads()))
for n in 3:2:(order+1)
ProgressMeter.next!(pmeter; showvalues = ())
ProgressMeter.next!(pmeter; showvalues = ()) # twice because of 2-step
proj11, proj10 = iterateKPM!(ket0, h', ket1, bandbracket, thread_buffers)
mulist[n] += 2 * proj11 - μ0
n + 1 > order + 1 && break
mulist[n + 1] += 2 * proj10 - μ1
ket0, ket1 = ket1, ket0
end
A.λ ≈ 1 || (mulist .*= A.λ)
return mulist
end

function mulscaled!(y, h´, x, (center, halfwidth))
mul!(y, h´, x)
invhalfwidth = 1/halfwidth
@. y = (y - center * x) * invhalfwidth
return y
end

function iterateKPM!(ket0, h´, ket1, (center, halfwidth), buff = ())
α = 2/halfwidth
β = 2center/halfwidth
mul!(ket0, h´, ket1, α, -1)
@. ket0 = ket0 - β * ket1
return proj_or_nothing(buff, ket0, ket1)
end

proj_or_nothing(::Tuple{}, ket0, ket1) = nothing
proj_or_nothing(buff, ket0, ket1) = (proj(ket1, ket1), proj(ket1, ket0))

# function iterateKPM!(ket0, h´, ket1, (center, halfwidth), thread_buffers = ())
# h = parent(h´)
# nz = nonzeros(h)
# rv = rowvals(h)
# α = -2 * center / halfwidth
# β = 2 / halfwidth
# reset_buffers!(thread_buffers)
# @threads for row in 1:size(ket0, 1)
# ptrs = nzrange(h, row)
# @inbounds for col in 1:size(ket0, 2)
# k1 = ket1[row, col]
# tmp = α * k1 - ket0[row, col]
# for ptr in ptrs
# tmp += β * adjoint(nz[ptr]) * ket1[rv[ptr], col]
# end
# # |k0⟩ → (⟨k1|2h - ⟨k0|)' = 2h'|k1⟩ - |k0⟩
# ket0[row, col] = tmp
# update_buffers!(thread_buffers, k1, tmp)
# end
# end
# return sum_buffers(thread_buffers)
# end

# reset_buffers!(::Tuple{}) = nothing
# function reset_buffers!((q, q´))
# fill!(q, zero(eltype(q)))
# fill!(q´, zero(eltype(q´)))
# return nothing
# end

# @inline update_buffers!(::Tuple{}, k1, tmp) = nothing
# @inline function update_buffers!((q, q´), k1, tmp)
# q[threadid()] += dot(k1, k1)
# q´[threadid()] += dot(tmp, k1)
# return nothing
# end

# @inline sum_buffers(::Tuple{}) = nothing
# @inline sum_buffers((q, q´)) = (sum(q), sum(q´))

# This is equivalent to tr(ket1'*ket2) for matrices, and ket1'*ket2 for vectors
proj(ket1, ket2) = dot(vec(ket1), vec(ket2))
Expand All @@ -206,12 +226,14 @@ function jackson!(μ::AbstractVector)
end

function bandbracketKPM(h, ::Missing)
@warn "Computing spectrum bounds... Consider using the `bandrange` kwargs for faster performance."
bandbracketKPM(h, bandrangeKPM(h))
end
bandbracketKPM(h, (ϵmin, ϵmax)::Tuple{T,T}, pad = float(T)(0.01)) where {T} = ((ϵmax + ϵmin) / 2, (ϵmax - ϵmin) / (2 - pad))

bandrangeKPM(h::Hamiltonian) = bandrangeKPM(matrixKPM(h, ArnoldiMethodPackage()))

function bandrangeKPM(h::AbstractMatrix{T}) where {T}
@warn "Computing spectrum bounds... Consider using the `bandrange` option for faster performance."
checkloaded(:ArnoldiMethod)
R = real(T)
decompl, _ = Main.ArnoldiMethod.partialschur(h, nev=1, tol=1e-4, which = Main.ArnoldiMethod.LR());
Expand Down Expand Up @@ -245,13 +267,11 @@ Same as above with momenta `μ` as input.
Equivalent to `dosKPM(bloch(h); kw...)` for finite hamiltonians (zero dimensional).
"""
dosKPM(h::AbstractMatrix; resolution = 2, kw...) =
dosKPM(h; resolution = 2, kw...) =
dosKPM(momentaKPM(h; kw...), resolution = resolution)

dosKPM(μ::MomentaKPM; resolution = 2) = real.(densityKPM(μ; resolution = resolution))

dosKPM(h::Hamiltonian; kw...) = dosKPM(matrixKPM(h); kw...)

"""
densityKPM(h::AbstractMatrix, A; resolution = 2, kw...)
Expand All @@ -271,17 +291,15 @@ Same as above with the KPM momenta as input (see `momentaKPM`).
Equivalent to `densityKPM(bloch(h), bloch(A); kw...)` for finite Hamiltonians (zero dimensional).
"""
densityKPM(h::AbstractMatrix, A; resolution = 2, kw...) =
densityKPM(h, A; resolution = 2, kw...) =
densityKPM(momentaKPM(h, A; kw...); resolution = resolution)

densityKPM(h::Hamiltonian, A::Hamiltonian; kw...) = densityKPM(matrixKPM(h), matrixKPM(A); kw...)

function densityKPM(momenta::MomentaKPM{T}; resolution = 2) where {T}
checkloaded(:FFTW)
(center, halfwidth) = momenta.bandbracket
numpoints = round(Int, length(momenta.μlist) * resolution)
numpoints = round(Int, length(momenta.mulist) * resolution)
ρlist = zeros(T, numpoints)
copyto!(ρlist, momenta.μlist)
copyto!(ρlist, momenta.mulist)
Main.FFTW.r2r!(ρlist, Main.FFTW.REDFT01, 1) # DCT-III in FFTW
xk = [cos(π * (k + 0.5) / numpoints) for k in 0:numpoints - 1]
@. ρlist = center + halfwidth * ρlist / (π * sqrt(1.0 - xk^2))
Expand All @@ -307,18 +325,16 @@ Same as above with the KPM momenta as input (see `momentaKPM`).
Equivalent to `averageKPM(bloch(h), bloch(A); kw...)` for finite Hamiltonians (zero
dimensional).
"""
averageKPM(h::AbstractMatrix, A; kBT = 0, Ef = 0, kw...) = averageKPM(momentaKPM(h, A; kw...); kBT = kBT, Ef = Ef)

averageKPM(h::Hamiltonian, A::Hamiltonian; kw...) = averageKPM(matrixKPM(h), matrixKPM(A); kw...)
averageKPM(h, A; kBT = 0, Ef = 0, kw...) = averageKPM(momentaKPM(h, A; kw...); kBT = kBT, Ef = Ef)

function averageKPM(momenta::MomentaKPM{T}; kBT = 0.0, Ef = 0.0) where {T}
(center, halfwidth) = momenta.bandbracket
order = length(momenta.μlist) - 1
order = length(momenta.mulist) - 1
# if !iszero(kBT)
# @warn "Finite temperature requires numerical evaluation of the integrals"
# checkloaded(:QuadGK)
# end
average = sum(n -> momenta.μlist[n + 1] * fermicheby(n, Ef, kBT, center, halfwidth), 0:order)
average = sum(n -> momenta.mulist[n + 1] * fermicheby(n, Ef, kBT, center, halfwidth), 0:order)
return average
end

Expand Down
14 changes: 14 additions & 0 deletions src/diagonalizer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,20 @@ end

similarmatrix(h, ::ArpackPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})

## ArnoldiMethod ##
struct ArnoldiMethodPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
kw::K
end

ArnoldiMethodPackage(; kw...) = (checkloaded(:ArnoldiMethod); ArnoldiMethodPackage(values(kw)))

function diagonalize(matrix, method::ArnoldiMethodPackage)
ϵ, ψ = Main.ArnoldiMethod.partialschur(matrix; (method.kw)...)
return ϵ, ψ
end

similarmatrix(h, ::ArnoldiMethodPackage) = similarmatrix(h, SparseMatrixCSC{blockeltype(h)})

## IterativeSolvers ##

struct KrylovKitPackage{K<:NamedTuple} <: AbstractDiagonalizeMethod
Expand Down
2 changes: 2 additions & 0 deletions src/hamiltonian.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -837,6 +837,8 @@ function similarmatrix(h, ::Type{A´} = matrixtype(h)) where {A´<:AbstractMatri
return _similarmatrix(parent(h), matrixtype(h), A´)
end

similarmatrix(h, ::Missing) = similarmatrix(h)

# We only provide the type combinastions that make sense
_similarmatrix(h, ::Type{A}, ::Type{A´}) where {A´,A<:A´} =
similar(h.harmonics[1].h)
Expand Down
1 change: 1 addition & 0 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,5 @@ using Quantica
include("test_hamiltonian.jl")
include("test_mesh.jl")
include("test_bandstructure.jl")
include("test_KPM.jl")
end
Loading

0 comments on commit 38c69dd

Please sign in to comment.