Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

#

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Jump to bottom

add Fe-based 5 orbital superconductor and graphene #5

Merged
merged 1 commit into from
Mar 9, 2022
Merged

add Fe-based 5 orbital superconductor and graphene #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
205 changes: 205 additions & 0 deletions example/FebasedSC.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
88 changes: 88 additions & 0 deletions example/graphene.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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("graphene_band.png")
end