Skip to content

Commit

Permalink
Add checkoptions for EnergyBands, BerryCurvature, `FermiSurface…
Browse files Browse the repository at this point in the history
…`, `DensityOfStates` and `InelasticNeutronScatteringSpectra`.
  • Loading branch information
waltergu committed Dec 3, 2024
1 parent 745269b commit e2f8b50
Show file tree
Hide file tree
Showing 8 changed files with 98 additions and 49 deletions.
1 change: 0 additions & 1 deletion .github/workflows/CI.yml
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@ jobs:
fail-fast: false
matrix:
version:
- '1.9'
- '1.10'
- '1.11'
os:
Expand Down
2 changes: 1 addition & 1 deletion Manifest.toml
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

julia_version = "1.11.1"
manifest_format = "2.0"
project_hash = "691253e15f09d19c3f3c1735ca8559810a4ed375"
project_hash = "d1819e50fb7295b196f9e8cb59c99c06975dcc53"

[[deps.Adapt]]
deps = ["LinearAlgebra", "Requires"]
Expand Down
6 changes: 3 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "TightBindingApproximation"
uuid = "d3576f78-cb9d-4774-9dc5-b247cf392ff4"
authors = ["waltergu <[email protected]>"]
version = "0.2.0"
version = "0.2.1"

[deps]
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
Expand All @@ -16,8 +16,8 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
[compat]
DelimitedFiles = "1.9 - 1"
Optim = "1.7"
QuantumLattices = "0.10.1"
QuantumLattices = "0.10.2"
RecipesBase = "1.2"
StaticArrays = "1.4 - 1"
TimerOutputs = "0.5"
julia = "1.9 - 1.11"
julia = "1.10 - 1.11"
2 changes: 1 addition & 1 deletion docs/src/examples/Phonon.md
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ The inelastic neutron scattering spectra of phonons can also be computed:
# scale: the scale of the intensity
spectra = phonon(
:INSS,
InelasticNeutronScatteringSpectra(path, range(0.0, 2.5, length=501); fwhm=0.05, scale=log)
InelasticNeutronScatteringSpectra(path, range(0.0, 2.5, length=501); fwhm=0.05, rescale=log)
)
plt = plot()
plot!(plt, spectra)
Expand Down
4 changes: 2 additions & 2 deletions docs/src/examples/Superconductor.md
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ t = Hopping(:t, 1.0, 1)
## p+ip pairing term
Δ = Pairing(
:Δ, Complex(0.5), 1, Coupling(:, 𝕗, :, :, (1, 1));
:Δ, Complex(0.5), 1, Coupling(𝕗, :, :, :, (1, 1));
amplitude=bond->exp(im*azimuth(rcoordinate(bond)))
)
Expand Down Expand Up @@ -135,7 +135,7 @@ hilbert = Hilbert(site=>Fock{:f}(1, 1) for site=1:length(unitcell))
t = Hopping(:t, symbols("t", real=true), 1)
μ = Onsite(:μ, symbols("μ", real=true))
Δ = Pairing(
:Δ, symbols("Δ", real=true), 1, Coupling(:, 𝕗, :, :, (1, 1));
:Δ, symbols("Δ", real=true), 1, Coupling(𝕗, :, :, :, (1, 1));
amplitude=bond->exp(im*azimuth(rcoordinate(bond)))
)
Expand Down
88 changes: 69 additions & 19 deletions src/Core.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@ using LinearAlgebra: Diagonal, Eigen, cholesky, dot, inv, norm, logdet, normaliz
using Printf: @printf, @sprintf
using QuantumLattices: atol, lazy, plain, rtol
using QuantumLattices: AbstractLattice, Action, Algorithm, Assignment, BrillouinZone, Boundary, CoordinatedIndex, Elastic, FockIndex, Fock, Formula, Frontend, Generator, Hilbert, Hooke, Hopping, ID, Index, Internal, Kinetic, LinearTransformation, Matrixization, Metric, Neighbors, Onsite, Operator, OperatorPack, OperatorSet, OperatorSum, OperatorUnitToTuple, Pairing, Phonon, PhononIndex, ReciprocalPath, ReciprocalSpace, ReciprocalZone, Term
using QuantumLattices: , bonds, dtype, expand, icoordinate, idtype, isannihilation, iscreation, optype, parametertype, rank, rcoordinate, shape, shrink, statistics, tostr, volume
using QuantumLattices: , bonds, checkoptions, dtype, expand, icoordinate, idtype, isannihilation, iscreation, optype, parametertype, rank, rcoordinate, shape, shrink, statistics, tostr, volume
using RecipesBase: RecipesBase, @recipe, @series, @layout
using TimerOutputs: TimerOutput, @timeit_debug

import LinearAlgebra: eigen, eigvals, eigvecs, ishermitian, Hermitian
import QuantumLattices: Parameters, Table, add!, dimension, getcontent, initialize, kind, matrix, parameternames, run!, update!
import QuantumLattices: Parameters, Table, add!, dimension, getcontent, initialize, kind, matrix, options, parameternames, run!, update!

const tbatimer = TimerOutput()

Expand Down Expand Up @@ -407,7 +407,7 @@ struct SimpleTBA{
H::Union{Formula, OperatorSet{<:Quadratic}, Generator{<:OperatorSet{<:Quadratic}}},
commutator::Union{AbstractMatrix, Nothing}
) where {K<:TBAKind}
@assert check_commutator(commutator) "SimpleTBA error: unsupported input commutator."
checkcommutator(commutator)
new{K, typeof(lattice), typeof(H), typeof(commutator)}(lattice, H, commutator)
end
end
Expand All @@ -417,12 +417,12 @@ end
end
return tba
end
@inline check_commutator(::Nothing) =true
function check_commutator(commutator::AbstractMatrix)
@inline checkcommutator(::Nothing) = nothing
function checkcommutator(commutator::AbstractMatrix)
values = eigvals(commutator)
num₁ = count(isapprox(+1, atol=atol, rtol=rtol), values)
num₂ = count(isapprox(-1, atol=atol, rtol=rtol), values)
return num₁==num₂==length(values)//2
@assert num₁==num₂==length(values)//2 "checkcommutator error: unsupported input commutator."
end

"""
Expand Down Expand Up @@ -456,7 +456,7 @@ struct CompositeTBA{
transformation::Quadraticization,
commutator::Union{AbstractMatrix, Nothing}
) where {K<:TBAKind}
@assert check_commutator(commutator) "SimpleTBA error: unsupported input commutator."
checkcommutator(commutator)
H = transformation(system)
new{K, typeof(lattice), typeof(system), typeof(transformation), typeof(H), typeof(commutator)}(lattice, system, transformation, H, commutator)
end
Expand Down Expand Up @@ -503,6 +503,19 @@ end
@inline wrapper(x) = (x,)
@inline wrapper(xs::Tuple) = xs

"""
const basicoptions = Dict(
:gauge => "gauge used to perform the Fourier transformation",
:infinitesimal => "infinitesimal added to the diagonal of the matrix representation of the tight-binding Hamiltonian"
)
Basic options of tight-binding actions.
"""
const basicoptions = Dict(
:gauge => "gauge used to perform the Fourier transformation",
:infinitesimal => "infinitesimal added to the diagonal of the matrix representation of the tight-binding Hamiltonian"
)

"""
EnergyBands{P<:ReciprocalSpace, L<:Union{Colon, AbstractVector{Int}}, O} <: Action
Expand All @@ -513,15 +526,21 @@ struct EnergyBands{P<:ReciprocalSpace, L<:Union{Colon, AbstractVector{Int}}, O}
bands::L
options::O
end
@inline EnergyBands(reciprocalspace::ReciprocalSpace, bands::Union{Colon, AbstractVector{Int}}=Colon(); options...) = EnergyBands(reciprocalspace, bands, options)
@inline options(::Type{<:EnergyBands}) = merge(basicoptions, Dict(
:tol => "maximum tolerance of the imaginary part of eigen energies"
))
@inline function EnergyBands(reciprocalspace::ReciprocalSpace, bands::Union{Colon, AbstractVector{Int}}=Colon(); options...)
checkoptions(EnergyBands; options...)
return EnergyBands(reciprocalspace, bands, options)
end
@inline initialize(eb::EnergyBands{P, Colon}, tba::TBA) where {P<:ReciprocalSpace} = (eb.reciprocalspace, zeros(Float64, length(eb.reciprocalspace), dimension(tba)))
@inline initialize(eb::EnergyBands{P, <:AbstractVector{Int}}, ::TBA) where {P<:ReciprocalSpace} = (eb.reciprocalspace, zeros(Float64, length(eb.reciprocalspace), length(eb.bands)))
@inline Base.nameof(tba::Algorithm{<:TBA}, eb::Assignment{<:EnergyBands}) = @sprintf "%s_%s" repr(tba; context=:select=>(names(eb.action.reciprocalspace))) eb.id
function run!(tba::Algorithm{<:TBA}, eb::Assignment{<:EnergyBands})
atol = get(eb.action.options, :atol, 10^-12)
tol = get(eb.action.options, :tol, 10^-12)
for (i, k) in enumerate(eb.action.reciprocalspace)
eigenvalues = eigvals(tba, k; eb.action.options...)[eb.action.bands]
norm(imag(eigenvalues))>atol && @warn("run! warning: imaginary eigen energies at $k with the norm of all imaginary parts being $(norm(imag(eigenvalues))).")
norm(imag(eigenvalues))>tol && @warn("run! warning: imaginary eigen energies at $k with the norm of all imaginary parts being $(norm(imag(eigenvalues))).")
eb.data[2][i, :] = real(eigenvalues)
end
end
Expand Down Expand Up @@ -588,9 +607,19 @@ struct BerryCurvature{B<:ReciprocalSpace, M<:BerryCurvatureMethod, O} <: Action
method::M
options::O
end
@inline BerryCurvature(reciprocalspace::ReciprocalSpace, method::BerryCurvatureMethod; gauge=:rcoordinate, options...) = BerryCurvature(reciprocalspace, method, (gauge=gauge, options...))
@inline BerryCurvature(reciprocalspace::ReciprocalSpace, μ::Real, d::Real=0.1, kx::T=nothing, ky::T=nothing; gauge=:rcoordinate, options...) where {T<:Union{Nothing, Vector{Float64}}} = BerryCurvature(reciprocalspace, Kubo(μ, d, kx, ky), (gauge=gauge, options...))
@inline BerryCurvature(reciprocalspace::Union{BrillouinZone, ReciprocalZone}, bands::AbstractVector{Int}, abelian::Bool=true; gauge=:rcoordinate, options...) = BerryCurvature(reciprocalspace, Fukui(bands; abelian=abelian), (gauge=gauge, options...))
@inline options(::Type{<:BerryCurvature}) = basicoptions
@inline function BerryCurvature(reciprocalspace::ReciprocalSpace, method::BerryCurvatureMethod; gauge=:rcoordinate, options...)
checkoptions(BerryCurvature; options...)
return BerryCurvature(reciprocalspace, method, (gauge=gauge, options...))
end
@inline function BerryCurvature(reciprocalspace::ReciprocalSpace, μ::Real, d::Real=0.1, kx::T=nothing, ky::T=nothing; gauge=:rcoordinate, options...) where {T<:Union{Nothing, Vector{Float64}}}
checkoptions(BerryCurvature; options...)
return BerryCurvature(reciprocalspace, Kubo(μ, d, kx, ky), (gauge=gauge, options...))
end
@inline function BerryCurvature(reciprocalspace::Union{BrillouinZone, ReciprocalZone}, bands::AbstractVector{Int}, abelian::Bool=true; gauge=:rcoordinate, options...)
checkoptions(BerryCurvature; options...)
return BerryCurvature(reciprocalspace, Fukui(bands; abelian=abelian), (gauge=gauge, options...))
end

# For the Berry curvature and Chern number (Berry phase ÷ 2π) on the first Brillouin zone
@inline function initialize(bc::BerryCurvature{<:BrillouinZone, <:Fukui}, ::TBA)
Expand Down Expand Up @@ -783,7 +812,11 @@ struct FermiSurface{B<:Union{BrillouinZone, ReciprocalZone}, A<:Union{Colon, Abs
orbitals::L
options::O
end
@inline options(::Type{<:FermiSurface}) = merge(basicoptions, Dict(
:fwhm => "full width at half maximum for the Gaussian broadening"
))
@inline function FermiSurface(reciprocalspace::Union{BrillouinZone, ReciprocalZone}, μ::Real=0.0, bands::Union{Colon, AbstractVector{Int}}=:, orbitals::Union{Colon, AbstractVector{Int}}...=:; options...)
checkoptions(FermiSurface; options...)
return FermiSurface(reciprocalspace, convert(Float64, μ), bands, orbitals, options)
end
function initialize(fs::FermiSurface, ::TBA)
Expand Down Expand Up @@ -831,7 +864,16 @@ struct DensityOfStates{B<:BrillouinZone, A<:Union{Colon, AbstractVector{Int}}, L
orbitals::L
options::O
end
@inline DensityOfStates(brillouinzone::BrillouinZone, bands::Union{Colon, AbstractVector{Int}}=:, orbitals::Union{Colon, AbstractVector{Int}}...=:; options...) = DensityOfStates(brillouinzone, bands, orbitals, options)
@inline options(::Type{<:DensityOfStates}) = merge(basicoptions, Dict(
:fwhm => "full width at half maximum for the Gaussian broadening",
:ne => "number of energy sample points",
:emin => "minimum value of the energy window",
:emax => "maximum value of the energy window"
))
@inline function DensityOfStates(brillouinzone::BrillouinZone, bands::Union{Colon, AbstractVector{Int}}=:, orbitals::Union{Colon, AbstractVector{Int}}...=:; options...)
checkoptions(DensityOfStates; options...)
return DensityOfStates(brillouinzone, bands, orbitals, options)
end
@inline function initialize(dos::DensityOfStates, ::TBA)
ne = get(dos.options, :ne, 100)
x = zeros(Float64, ne)
Expand Down Expand Up @@ -870,7 +912,15 @@ struct InelasticNeutronScatteringSpectra{P<:ReciprocalSpace, E<:AbstractVector,
new{typeof(reciprocalspace), typeof(energies), typeof(options)}(reciprocalspace, energies, options)
end
end
@inline InelasticNeutronScatteringSpectra(reciprocalspace::ReciprocalSpace, energies::AbstractVector; options...) = InelasticNeutronScatteringSpectra(reciprocalspace, energies, options)
@inline options(::Type{<:InelasticNeutronScatteringSpectra}) = merge(basicoptions, Dict(
:fwhm => "full width at half maximum for the Gaussian broadening",
:check => "whether the polarization consistency of phonons will be checked",
:rescale => "function used to rescale the intensity of the spectrum at each energy-momentum point"
))
@inline function InelasticNeutronScatteringSpectra(reciprocalspace::ReciprocalSpace, energies::AbstractVector; options...)
checkoptions(InelasticNeutronScatteringSpectra; options...)
return InelasticNeutronScatteringSpectra(reciprocalspace, energies, options)
end
@inline initialize(inss::InelasticNeutronScatteringSpectra, ::TBA) = (inss.reciprocalspace, inss.energies, zeros(Float64, length(inss.energies), length(inss.reciprocalspace)))

# Inelastic neutron scattering spectra for phonons.
Expand All @@ -881,7 +931,7 @@ function run!(tba::Algorithm{<:CompositeTBA{Phononic, <:AbstractLattice}}, inss:
sequences = Dict(site=>[tba.frontend.transformation.table[Index(site, PhononIndex{:u}(Char(Int('x')+i-1)))] for i=1:phonon.ndirection] for (site, phonon) in pairs(tba.frontend.system.hilbert))
eigenvalues, eigenvectors = eigen(tba, inss.action.reciprocalspace; inss.action.options...)
for (i, (momentum, values, vectors)) in enumerate(zip(inss.action.reciprocalspace, eigenvalues, eigenvectors))
check && @timeit_debug tba.timer "check" check_polarizations(@views(vectors[(dim÷2+1):dim, 1:(dim÷2)]), @views(vectors[(dim÷2+1):dim, dim:-1:(dim÷2+1)]), momentum./pi)
check && @timeit_debug tba.timer "check" checkpolarizations(@views(vectors[(dim÷2+1):dim, 1:(dim÷2)]), @views(vectors[(dim÷2+1):dim, dim:-1:(dim÷2+1)]), momentum./pi)
@timeit_debug tba.timer "spectra" begin
for j = 1:dim
factor = 0
Expand All @@ -896,11 +946,11 @@ function run!(tba::Algorithm{<:CompositeTBA{Phononic, <:AbstractLattice}}, inss:
end
end
end
inss.data[3][:, :] = get(inss.action.options, :scale, identity).(inss.data[3].+1)
inss.data[3][:, :] = get(inss.action.options, :rescale, identity).(inss.data[3].+1)
end
@inline function check_polarizations(qs₁::AbstractMatrix, qs₂::AbstractMatrix, momentum)
@inline function checkpolarizations(qs₁::AbstractMatrix, qs₂::AbstractMatrix, momentum)
inner = mapreduce((e₁, e₂)->norm(conj(e₁)*e₂), +, qs₁, qs₂)/norm(qs₁)/norm(qs₂)
isapprox(inner, 1; atol=100*atol, rtol=100*rtol) || begin
@warn("check_polarizations: small inner product $inner at π*$momentum, indication of degeneracy, otherwise inconsistent polarization vectors.")
@warn("checkpolarizations: small inner product $inner at π*$momentum, indication of degeneracy, otherwise inconsistent polarization vectors.")
end
end
30 changes: 15 additions & 15 deletions test/Manifest.toml
Original file line number Diff line number Diff line change
Expand Up @@ -210,9 +210,9 @@ version = "1.3.7"

[[deps.FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
git-tree-sha1 = "fa8e19f44de37e225aa0f1695bc223b05ed51fb4"
git-tree-sha1 = "786e968a8d2fb167f2e4880baba62e0e26bd8e4e"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.13.3+0"
version = "2.13.3+1"

[[deps.FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl"]
Expand Down Expand Up @@ -399,10 +399,10 @@ uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.11.0+0"

[[deps.Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "6f73d1dd803986947b2c750138528a999a6c7733"
deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "ff3b4b9d35de638936a525ecd36e86a8bb919d11"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.6.0+0"
version = "1.7.0+0"

[[deps.Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl"]
Expand Down Expand Up @@ -871,9 +871,9 @@ version = "0.4.1"

[[deps.Unitful]]
deps = ["Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "d95fe458f26209c66a187b1114df96fd70839efd"
git-tree-sha1 = "01915bfcd62be15329c9a07235447a89d588327c"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.21.0"
version = "1.21.1"

[deps.Unitful.extensions]
ConstructionBaseUnitfulExt = "ConstructionBase"
Expand Down Expand Up @@ -944,9 +944,9 @@ version = "1.2.4+0"

[[deps.Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495"
git-tree-sha1 = "9dafcee1d24c4f024e7edc92603cedba72118283"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.8.6+0"
version = "1.8.6+1"

[[deps.Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl"]
Expand All @@ -968,9 +968,9 @@ version = "1.1.4+1"

[[deps.Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
git-tree-sha1 = "d2d1a5c49fae4ba39983f63de6afcbea47194e85"
git-tree-sha1 = "d7155fea91a4123ef59f42c4afb5ab3b4ca95058"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.6+0"
version = "1.3.6+1"

[[deps.Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
Expand Down Expand Up @@ -1010,9 +1010,9 @@ version = "0.1.1+1"

[[deps.Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "bcd466676fef0878338c61e655629fa7bbc69d8e"
git-tree-sha1 = "1a74296303b6524a0472a8cb12d3d87a78eb3612"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.17.0+0"
version = "1.17.0+1"

[[deps.Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
Expand Down Expand Up @@ -1099,9 +1099,9 @@ version = "0.56.3+0"

[[deps.gperf_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3516a5630f741c9eecb3720b1ec9d8edc3ecc033"
git-tree-sha1 = "0ba42241cb6809f1a278d0bcb976e0483c3f1f2d"
uuid = "1a1c6b14-54f6-533d-8383-74cd7377aa70"
version = "3.1.1+0"
version = "3.1.1+1"

[[deps.libaom_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl"]
Expand Down
Loading

0 comments on commit e2f8b50

Please sign in to comment.