Merge branch 'master' into compathelper/new_version/2022-03-09-00-33-…

…28-885-02818818054

.gitignore

@@ -29,4 +29,7 @@ Manifest.toml
 
 *.DS_Store
 
-.ipynb_checkpoints
+.ipynb_checkpoints
+*.png
+*.pdf
+run

Project.toml

@@ -5,11 +5,13 @@ version = "0.1.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 StaticArrays = "1"
+AbstractTrees = "0.3"
 julia = "1.6"
 
 [extras]

example/FebasedSC.jl

@@ -0,0 +1,205 @@
+using TightBinding
+using LinearAlgebra
+using Plots
+la = set_Lattice(2,[[1,0],[0,1]])
+add_atoms!(la,[0,0])
+add_atoms!(la,[0,0])
+add_atoms!(la,[0,0])
+add_atoms!(la,[0,0])
+add_atoms!(la,[0,0])
+
+tmat = [
+-0.7    0 -0.4  0.2 -0.1
+-0.8    0    0    0    0
+ 0.8 -1.5    0    0 -0.3
+   0  1.7    0    0 -0.1
+-3.0    0    0 -0.2    0
+-2.1  1.5    0    0    0
+ 1.3    0  0.2 -0.2    0
+ 1.7    0    0  0.2    0
+-2.5  1.4    0    0    0
+-2.1  3.3    0 -0.3  0.7
+ 1.7  0.2    0  0.2    0
+ 2.5    0    0  0.3    0
+ 1.6  1.2 -0.3 -0.3 -0.3
+   0    0    0 -0.1    0
+ 3.1 -0.7 -0.2    0    0
+]
+tmat = 0.1.*tmat
+imap = zeros(Int64,5,5)
+count = 0
+for μ=1:5
+    for ν=μ:5
+        global count += 1
+        imap[μ,ν] = count
+    end
+end
+Is = [1,-1,-1,1,1,1,1,-1,-1,1,-1,-1,1,1,1]
+σds = [1,-1,1,1,-1,1,-1,-1,1,1,1,-1,1,-1,1]
+tmat_σy = tmat[:,:]
+tmat_σy[imap[1,2],:] = -tmat[imap[1,3],:]
+tmat_σy[imap[1,3],:] = -tmat[imap[1,2],:]
+tmat_σy[imap[1,4],:] = -tmat[imap[1,4],:]
+tmat_σy[imap[2,2],:] = tmat[imap[3,3],:]
+tmat_σy[imap[2,4],:] = tmat[imap[3,4],:]
+tmat_σy[imap[2,5],:] = -tmat[imap[3,5],:]
+tmat_σy[imap[3,3],:] = tmat[imap[2,2],:]
+tmat_σy[imap[3,4],:] = tmat[imap[2,4],:]
+tmat_σy[imap[3,5],:] = -tmat[imap[2,5],:]
+tmat_σy[imap[4,5],:] = -tmat[imap[4,5],:]
+
+hoppingmatrix = zeros(Float64,5,5,5,5)
+hops = [-2,-1,0,1,2]
+hopelements = [[1,0],[1,1],[2,0],[2,1],[2,2]]
+
+for μ = 1:5
+    for ν=μ:5
+        for ii=1:5
+            ihop = hopelements[ii][1]
+            jhop = hopelements[ii][2]
+            #[a,b],[a,-b],[-a,-b],[-a,b],[b,a],[b,-a],[-b,a],[-b,-a]
+
+            #[a,b]
+            i = ihop +3
+            j = jhop +3
+            hoppingmatrix[μ,ν,i,j]=tmat[imap[μ,ν],ii]
+            #[a,-b] = σy*[a,b] [1,1] -> [1,-1]
+            if jhop != 0
+                i = ihop +3
+                j = -jhop +3
+                hoppingmatrix[μ,ν,i,j]=tmat_σy[imap[μ,ν],ii]
+            end
+
+            if μ != ν
+                #[-a,-b] = I*[a,b] [1,1] -> [-1,-1],[1,0]->[-1,0]
+                i = -ihop +3
+                j = -jhop +3
+                hoppingmatrix[μ,ν,i,j]=Is[imap[μ,ν]]*tmat[imap[μ,ν],ii]
+                #[-a,b] = I*[a,-b] = I*σy*[a,b]  #[2,0]->[-2,0]
+                if jhop != 0
+                    i = -ihop +3
+                    j = jhop +3
+                    hoppingmatrix[μ,ν,i,j]=Is[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
+                end
+            end
+            #[b,a],[b,-a],[-b,a],[-b,-a]
+            if jhop != ihop
+                #[b,a] = σd*[a,b]
+                i = jhop +3
+                j = ihop +3
+                hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*tmat[imap[μ,ν],ii]
+                #[-b,a] = σd*σy*[a,b]
+                if jhop != 0
+                    i = -jhop +3
+                    j = ihop +3
+                    hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
+                end
+
+                if μ != ν
+                    #[-b,-a] = σd*[-a,-b] = σd*I*[a,b]
+                    i = -jhop +3
+                    j = -ihop +3
+                    hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*Is[imap[μ,ν]]*tmat[imap[μ,ν],ii]
+                    #[b,-a] = σd*[-a,b] = σd*I*[a,-b] = σd*I*σy*[a,b]  #[2,0]->[-2,0]
+                    if jhop != 0
+                        i = jhop +3
+                        j = -ihop +3
+                        hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*Is[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
+                    end
+                end
+            end
+        end
+
+
+    end
+end
+
+for μ=1:5
+    for ν=μ:5
+        for i = 1:5
+            ih = hops[i]
+            for j = 1:5
+                jh = hops[j]
+                if hoppingmatrix[μ,ν,i,j] != 0.0                
+                    add_hoppings!(la,hoppingmatrix[μ,ν,i,j],μ,ν,[ih,jh])
+                end
+            end
+        end
+    end
+end
+
+onsite = [10.75,10.96,10.96,11.12,10.62]
+set_onsite!(la,onsite)
+
+ham = hamiltonian_k(la) #Construct the Hamiltonian
+function dispersion_Fe(k)
+    eigen_graphene = eigen(ham(k))
+    energy = eigen_graphene.values
+    U_matrix = eigen_graphene.vectors
+    return energy, U_matrix
+end
+
+
+"""
+This function returns the diagonalized green's function G_diag of Fe-based superconductor, together with the U matrix which can reproduce the original Green's function by:
+    G = U*G_diag*inv(U)
+"""
+function green_Fe(k, μ, ω)
+    energies, U = dispersion_Fe(k)
+    Green = zeros(ComplexF64, length(energies), length(energies))
+    for i=1:length(energies)
+        Green[i,i]= -1/( im*ω - energies[i] + μ) 
+    end
+    return Green, U
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    k = [π/4,π/4]
+    energies, U = dispersion_Fe(k)
+    #U = transpose(U)
+    Ematrix = zeros(Float64, length(energies), length(energies))
+    for i=1:length(energies)
+        Ematrix[i,i]=energies[i]
+    end
+    dim = la.dim
+    n = la.numatoms
+    print("Error of matrix diagonalization is ", maximum(abs.(U*Ematrix*inv(U)-ham(k))),"\n")
+    nk = 100
+    energies = zeros(Float64,n,nk,nk)
+    kx = range(-2π, stop = 2π, length = nk)
+    ky = range(-2π, stop = 2π, length = nk)
+    for i1=1:nk
+        for i2=1:nk
+            energy = dispersion(dim,n,ham,[kx[i1],ky[i2]])
+            for j=1:n
+                energies[j,i1,i2] = energy[j]
+            end
+
+        end
+    end
+    pls = heatmap(kx, ky, energies[1,:,:])
+    #pls = surface!(kx, ky, energies[2,:,:], legend=false, fillalpha=0.5)
+    #pls = plot_fermisurface_2D(la, Eshift = -0.2, nk = 1000)
+    savefig("Fe5_lowest_band.png")
+
+"""
+Test band structure of Fe-based 5 orbital superconductor.
+Should be consistent with T. Nomura, J. Phys. Soc. Jpn. 78, 034716 (2009)
+"""
+    set_μ!(la,10.96) #set the chemical potential
+    klines = set_Klines()
+    kmin = [0,0]
+    kmax = [π,0]
+    add_Kpoints!(klines,kmin,kmax,"(0,0)","(pi,0)",nk=nk)
+
+    kmin = [π,0]
+    kmax = [π,π]
+    add_Kpoints!(klines,kmin,kmax,"(pi,0)","(pi,pi)",nk=nk)
+
+    kmin = [π,π]
+    kmax = [0,0]
+    add_Kpoints!(klines,kmin,kmax,"(pi,pi)","(0,0)",nk=nk)
+
+    pls = calc_band_plot(klines,la)
+    savefig("Fe5band.png")
+end

example/graphene.jl

@@ -0,0 +1,88 @@
+using TightBinding
+using Plots
+using LinearAlgebra
+#Primitive vectors
+a1 = [sqrt(3) / 2, 1 / 2]
+a2 = [0, 1]
+#set lattice
+la = set_Lattice(2, [a1, a2])
+#add atoms
+add_atoms!(la, [1 / 3, 1 / 3])
+add_atoms!(la, [2 / 3, 2 / 3])
+show_neighbors(la)
+t = 1.0
+add_hoppings!(la, -t, 1, 2, [1 / 3, 1 / 3])
+add_hoppings!(la, -t, 1, 2, [-2 / 3, 1 / 3])
+add_hoppings!(la, -t, 1, 2, [1 / 3, -2 / 3])
+# plot_lattice_2d(la)
+nk = 100 #numer ob meshes. nk^d meshes are used. d is a dimension.
+# plot_DOS(la, nk)
+#show the band structure
+klines = set_Klines()
+kmin = [0, 0]
+kmax = [2π / sqrt(3), 0]
+add_Kpoints!(klines, kmin, kmax, "G", "K")
+
+kmin = [2π / sqrt(3), 0]
+kmax = [2π / sqrt(3), 2π / 3]
+add_Kpoints!(klines, kmin, kmax, "K", "M")
+
+kmin = [2π / sqrt(3), 2π / 3]
+kmax = [0, 0]
+add_Kpoints!(klines, kmin, kmax, "M", "G")
+# calc_band_plot(klines,la)
+
+ham = hamiltonian_k(la) #Construct the Hamiltonian
+
+
+function dispersion_graphene(k)
+    eigen_graphene = eigen(ham(k))
+    energy = eigen_graphene.values
+    U_matrix = eigen_graphene.vectors
+    return energy, U_matrix
+end
+
+
+"""
+This function returns the diagonalized green's function G_diag of graphene, together with the U matrix which can reproduce the original Green's function by:
+    G = U*G_diag*inv(U)
+"""
+function green_graphene(k, μ, ω)
+    energies, U = dispersion_graphene(k)
+    Green = zeros(ComplexF64, length(energies), length(energies))
+    for i = 1:length(energies)
+        Green[i, i] = -1 / (im * ω - energies[i] + μ)
+    end
+    return Green, U
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    k = [π / 4, π / 4]
+    energies, U = dispersion_graphene(k)
+    #U = transpose(U)
+    Ematrix = zeros(Float64, length(energies), length(energies))
+    for i = 1:length(energies)
+        Ematrix[i, i] = energies[i]
+    end
+    dim = la.dim
+    n = la.numatoms
+    print("Error of matrix diagonalization is ", maximum(abs.(U * Ematrix * inv(U) - ham(k))), "\n")
+    nk = 100
+    energies = zeros(Float64, n, nk, nk)
+    kx = range(-2π, stop=2π, length=nk)
+    ky = range(-2π, stop=2π, length=nk)
+    for i1 = 1:nk
+        for i2 = 1:nk
+            energy = dispersion(dim, n, ham, [kx[i1], ky[i2]])
+            for j = 1:n
+                energies[j, i1, i2] = energy[j]
+            end
+
+        end
+    end
+    pls = heatmap(kx, ky, energies[1, :, :])
+    #pls = surface(kx, ky, energies[1,:,:], legend=false, fillalpha=0.5)
+    #pls = surface!(kx, ky, energies[2,:,:], legend=false, fillalpha=0.5)
+    #pls = plot_fermisurface_2D(la, Eshift = -0.2, nk = 1000)
+    savefig("run/graphene_band.png")
+end

example/graphene_treetest.jl

@@ -0,0 +1,42 @@
+using SpaceGrid
+using SpaceGrid.AbstractTrees, SpaceGrid.GridTree, SpaceGrid.BaseMesh
+using Plots
+using LinearAlgebra
+
+include("graphene.jl")
+
+function density(k; band=1)
+    T = 0.1
+    μ = -2.5
+
+    dim = la.dim
+    n = la.numatoms
+    ϵ = dispersion(dim, n, ham, k)[band] - μ
+
+    # return 1 / (exp((ϵ) / T) + 1.0)
+    return (π * T)^2 / ((π * T)^2 + ϵ^2)
+    # return (exp((ϵ) / T) / T)/(exp((ϵ) / T) + 1.0)^2
+end
+
+latvec = [2 0; 1 sqrt(3)] .* (2π)
+# println(SpaceGrid.GridTree._calc_subpoints(3, [3,3], latvec, 2))
+tg = treegridfromdensity(k -> density(k), latvec; atol=1 / 2^18, maxdepth=9, mindepth=1, N=2)
+
+X, Y = zeros(Float64, size(tg)), zeros(Float64, size(tg))
+for (pi, p) in enumerate(tg)
+    X[pi], Y[pi] = p[1], p[2]
+    # println(p)
+end
+
+println("size:$(size(tg)),\t length:$(length(tg)),\t efficiency:$(efficiency(tg))")
+
+smap = SymMap(tg, k -> density(k); atol=1e-6)
+# println(smap.map)
+# println(smap.inv_map)
+println("compress:$(smap.reduced_length/length(smap.map))")
+
+k = 4π
+# p = plot(legend = false, size = (1024, 1024), xlim = (-1, 1), ylim = (-1, 1))
+p = plot(legend=false, size=(1024, 1024), xlim=(-k, k), ylim=(-k, k))
+scatter!(p, X, Y, marker=:cross, markersize=2)
+savefig(p, "run/tg.pdf")