From 1fa298a92eca94a8db2a9be76861971ab5167454 Mon Sep 17 00:00:00 2001 From: waltergu Date: Wed, 31 May 2023 16:58:47 +0800 Subject: [PATCH] Add `FermiSurface` and `DensityOfStates`. --- src/TightBindingApproximation.jl | 89 ++++++++++++++++++++++++++++++- test/TightBindingApproximation.jl | 46 +++++++++------- 2 files changed, 114 insertions(+), 21 deletions(-) diff --git a/src/TightBindingApproximation.jl b/src/TightBindingApproximation.jl index f6307f7..e277fa6 100644 --- a/src/TightBindingApproximation.jl +++ b/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 diff --git a/test/TightBindingApproximation.jl b/test/TightBindingApproximation.jl index 9a4230b..d0835e8 100644 --- a/test/TightBindingApproximation.jl +++ b/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,7 +100,7 @@ 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) @@ -108,27 +108,36 @@ end Δ = 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