Add FermiSurface and DensityOfStates.

src/TightBindingApproximation.jl

@@ -21,7 +21,7 @@ import QuantumLattices: contentnames
 
 export Bosonic, Fermionic, Phononic, TBAKind
 export AbstractTBA, TBA, TBAMatrix, TBAMatrixRepresentation, commutator
-export BerryCurvature, EnergyBands, InelasticNeutronScatteringSpectra
+export BerryCurvature, DensityOfStates, EnergyBands, FermiSurface, InelasticNeutronScatteringSpectra
 export SampleNode, deviation, optimize!
 
 const tbatimer = TimerOutput()
@@ -473,6 +473,93 @@ end
     pack[2].data[1:2]
 end
 
+function spectralfunction(tbakind::TBAKind, ω::Real, values::Vector{<:Real}, vectors::Matrix{<:Number}, orbitals::Union{Colon, Vector{Int}}=:; σ::Real)
+    result = zero(ω)
+    start = isa(tbakind, TBAKind{:TBA}) ? 1 : length(values)÷2
+    for i = start:length(values)
+        factor = mapreduce(abs2, +, vectors[orbitals, i])
+        result += factor*exp(-(ω-values[i])^2/2/σ^2)
+    end
+    return result/√(2pi)/σ
+end
+"""
+    FermiSurface{B<:Union{BrillouinZone, ReciprocalZone}, L<:Tuple{Vararg{Union{Colon, Vector{Int}}}}, O} <: Action
+
+Fermi surface of a free fermionic system.
+"""
+struct FermiSurface{B<:Union{BrillouinZone, ReciprocalZone}, L<:Tuple{Vararg{Union{Colon, Vector{Int}}}}, O} <: Action
+    reciprocalspace::B
+    μ::Float64
+    orbitals::L
+    options::O
+end
+@inline FermiSurface(reciprocalspace::Union{BrillouinZone, ReciprocalZone}, μ::Real=0.0, orbitals::Union{Colon, Vector{Int}}...=:; options...) = FermiSurface(reciprocalspace, μ, orbitals, options)
+function initialize(fs::FermiSurface, ::AbstractTBA)
+    @assert length(fs.reciprocalspace.reciprocals)==2 "initialize error: only two dimensional reciprocal spaces are supported."
+    ny, nx = map(length, shape(fs.reciprocalspace))
+    z = zeros(Float64, ny, nx, length(fs.orbitals))
+    return (fs.reciprocalspace, z)
+end
+function run!(tba::Algorithm{<:AbstractTBA{<:Fermionic{:TBA}}}, fs::Assignment{<:FermiSurface})
+    count = 1
+    σ = get(fs.action.options, :fwhm, 0.1)/2/√(2*log(2))
+    eigenvalues, eigenvectors = eigen(tba, fs.action.reciprocalspace)
+    ny, nx = map(length, shape(fs.action.reciprocalspace))
+    for i=1:nx, j=1:ny
+        for (k, orbitals) in enumerate(fs.action.orbitals)
+            fs.data[2][j, i, k] += spectralfunction(kind(tba.frontend), fs.action.μ, eigenvalues[count], eigenvectors[count], orbitals; σ=σ)
+        end
+        count += 1
+    end
+end
+@recipe function plot(pack::Tuple{Algorithm{<:AbstractTBA}, Assignment{<:FermiSurface}})
+    if size(pack[2].data[2])[3]==1
+        title --> nameof(pack...)
+        titlefontsize --> 10
+        pack[2].data[1], pack[2].data[2][:, :, 1]
+    else
+        subtitles --> [@sprintf("orbitals: %s", (orbitals==:) ? "all" : join(orbitals, ", ")) for orbitals in pack[2].action.orbitals]
+        subtitlefontsize --> 8
+        plot_title --> nameof(pack[1], pack[2])
+        plot_titlefontsize --> 10
+        pack[2].data
+    end
+end
+
+"""
+    DensityOfStates{B<:BrillouinZone, L<:Tuple{Vararg{Union{Colon, Vector{Int}}}}, O} <: Action
+
+Density of states of a tight-binding system.
+"""
+struct DensityOfStates{B<:BrillouinZone, L<:Tuple{Vararg{Union{Colon, Vector{Int}}}}, O} <: Action
+    brillouinzone::B
+    orbitals::L
+    options::O
+end
+@inline DensityOfStates(brillouinzone::BrillouinZone, orbitals::Union{Colon, Vector{Int}}...=:; options...) = DensityOfStates(brillouinzone, orbitals, options)
+@inline function initialize(dos::DensityOfStates, ::AbstractTBA)
+    ne = get(dos.options, :ne, 100)
+    x = zeros(Float64, ne)
+    z = zeros(Float64, ne, length(dos.orbitals))
+    return (x, z)
+end
+function run!(tba::Algorithm{<:AbstractTBA{<:Fermionic{:TBA}}}, dos::Assignment{<:DensityOfStates})
+    σ = get(dos.action.options, :fwhm, 0.1)/2/√(2*log(2))
+    eigenvalues, eigenvectors = eigen(tba, dos.action.brillouinzone)
+    emin = get(dos.action.options, :emin, mapreduce(minimum, min, eigenvalues))
+    emax = get(dos.action.options, :emin, mapreduce(maximum, max, eigenvalues))
+    ne = get(dos.action.options, :ne, 100)
+    nk = length(dos.action.brillouinzone)
+    for (i, ω) in enumerate(LinRange(emin, emax, ne))
+        dos.data[1][i] = ω
+        for (j, orbitals) in enumerate(dos.action.orbitals)
+            for (values, vectors) in zip(eigenvalues, eigenvectors)
+                dos.data[2][i, j] += spectralfunction(kind(tba.frontend), ω, values, vectors, orbitals; σ=σ)/nk
+            end
+        end
+    end
+end
+
 """
     InelasticNeutronScatteringSpectra{P<:ReciprocalSpace, E<:AbstractVector, O} <: Action
 

test/TightBindingApproximation.jl

@@ -1,10 +1,10 @@
 using LinearAlgebra: Diagonal, Eigen, Hermitian, eigen, eigvals, eigvecs, ishermitian
 using Plots: plot, plot!, savefig
 using QuantumLattices: dimension, kind, matrix, update!
-using QuantumLattices: Coupling, Hilbert, Metric, OperatorUnitToTuple
+using QuantumLattices: Coupling, Hilbert, MatrixCoupling, Metric, OperatorUnitToTuple
 using QuantumLattices: Algorithm, Parameters
 using QuantumLattices: Elastic, FID, Fock, Hooke, Hopping, Kinetic, Onsite, Pairing, Phonon
-using QuantumLattices: BrillouinZone, Lattice, ReciprocalPath, ReciprocalZone, azimuth, rcoordinate, reciprocals, @rectangle_str
+using QuantumLattices: BrillouinZone, Lattice, ReciprocalPath, ReciprocalZone, azimuth, rcoordinate, reciprocals, @rectangle_str, @σ_str
 using QuantumLattices: contentnames
 using TightBindingApproximation
 
@@ -100,35 +100,44 @@ end
     @test isapprox(op.minimizer, [-1.0, 0.0], atol=10^-10)
 end
 
-@time @testset "plot" begin
+@time @testset "EnergyBands and BerryCurvature" begin
     unitcell = Lattice([0.0, 0.0]; name=:Square, vectors=[[1.0, 0.0], [0.0, 1.0]])
     hilbert = Hilbert(site=>Fock{:f}(1, 1) for site=1:length(unitcell))
     t = Hopping(:t, 1.0, 1)
     μ = Onsite(:μ, 3.5)
     Δ = Pairing(:Δ, Complex(0.5), 1, Coupling(:, FID, :, :, (1, 1)); amplitude=bond->exp(im*azimuth(rcoordinate(bond))))
     sc = Algorithm(Symbol("p+ip"), TBA(unitcell, hilbert, (t, μ, Δ)))
     @test eigen(sc) == Eigen(eigvals(sc), eigvecs(sc))
+
     path = ReciprocalPath(reciprocals(unitcell), rectangle"Γ-X-M-Γ", length=100)
     @test eigen(sc, path) == (eigvals(sc, path), eigvecs(sc, path))
-    energybands = sc(:EB, EnergyBands(path))
-    plt = plot(energybands)
-    display(plt)
-    savefig(plt, "eb.png")
+    savefig(plot(sc(:EB, EnergyBands(path))), "eb.png")
 
     brillouin = BrillouinZone(reciprocals(unitcell), 100)
-    berry = sc(:BerryCurvature, BerryCurvature(brillouin, [1, 2]));
-    plt = plot(berry)
-    display(plt)
-    savefig(plt, "bc.png")
+    savefig(plot(sc(:BerryCurvature, BerryCurvature(brillouin, [1, 2]))), "bc.png")
 
     reciprocalzone = ReciprocalZone(reciprocals(unitcell), [-2.0=>2.0, -2.0=>2.0]; length=201, ends=(true, true))
-    berry = sc(:BerryCurvatureExtended, BerryCurvature(reciprocalzone, [1, 2]))
-    plt = plot(berry)
-    display(plt)
-    savefig(plt, "bcextended.png")
+    savefig(plot(sc(:BerryCurvatureExtended, BerryCurvature(reciprocalzone, [1, 2]))), "bcextended.png")
+end
+
+@time @testset "FermiSurface and DensityOfStates" begin
+    unitcell = Lattice([0.0, 0.0]; name=:Square, vectors=[[1.0, 0.0], [0.0, 1.0]])
+    hilbert = Hilbert(site=>Fock{:f}(1, 2) for site=1:length(unitcell))
+    t = Hopping(:t, 1.0, 1)
+    h = Onsite(:h, 0.1, MatrixCoupling(:, FID, :, σ"z", :))
+    tba = Algorithm(:tba, TBA(unitcell, hilbert, (t, h)))
+
+    brillouin = BrillouinZone(reciprocals(unitcell), 200)
+    savefig(plot(tba(:FermiSurface, FermiSurface(brillouin, 0.0))), "fs-all.png")
+    savefig(plot(tba(Symbol("FermiSurface-SpinDependent"), FermiSurface(brillouin, 0.0, [1], [2]))), "fs-spin.png")
+    savefig(plot(tba(:DensityOfStates, DensityOfStates(brillouin, :, [1], [2]))), "dos.png")
+
+    reciprocalzone = ReciprocalZone(reciprocals(unitcell), [-2.0=>2.0, -2.0=>2.0]; length=401, ends=(true, true))
+    savefig(plot(tba(:FermiSurfaceExtended, FermiSurface(reciprocalzone, 0.0))), "fs-extended-all.png")
+    savefig(plot(tba(Symbol("FermiSurfaceExtended-SpinDependent"), FermiSurface(reciprocalzone, 0.0, [1], [2]))), "fs-extended-spin.png")
 end
 
-@time @testset "phonon" begin
+@time @testset "InelasticNeutronScatteringSpectra" begin
     unitcell = Lattice([0.0, 0.0]; name=:Square, vectors=[[1.0, 0.0], [0.0, 1.0]])
     hilbert = Hilbert(site=>Phonon(2) for site=1:length(unitcell))
     T = Kinetic(:T, 0.5)
@@ -138,14 +147,11 @@ end
     path = ReciprocalPath(reciprocals(unitcell), rectangle"Γ-X-M-Γ", length=100)
 
     energybands = phonon(:EB, EnergyBands(path))
-    plt = plot(energybands)
-    display(plt)
-    savefig(plt, "phonon.png")
+    savefig(plot(energybands), "phonon.png")
 
     inelastic = phonon(:INSS, InelasticNeutronScatteringSpectra(path, range(0.0, 2.5, length=501); fwhm=0.05, scale=log))
     plt = plot()
     plot!(plt, inelastic)
     plot!(plt, energybands, color=:white, linestyle=:dash)
-    display(plt)
     savefig("inelastic.png")
 end