最近刚刚学习了julia, 手头上也正好在重复一篇文章，就正好拿新学习的内容一边温习一边做研究。
导入函数库
# Import external package that used in program
import Pkg
# Pkg.add("PyPlot")
# Pkg.add("LinearAlgebra")
# Pkg.add("CPUTime")
#Pkg.add("SparseArrays")
#Pkg.add("DelimitedFiles")
函数定义
# Construct function
# Phase factor in Landu guage
function phi(x)
return exp(-im*B*x)
end
# --------------------------------------
function boundary(bry)
for m1 in 1:yn
for m2 in 1:xn
i = (m1 - 1)*xn + m2
bry[1,i] = i + 1
if m2 == xn
bry[1,i] = bry[1,i] - xn
end
bry[2,i] = i - 1
if m2 == 1
bry[2,i] = bry[2,i] + xn
end
bry[3,i] = i + xn
if m1 == yn
bry[3,i] = bry[3,i] - len2
end
bry[4,i] = i - xn
if m1 == 1
bry[4,i] = bry[4,i] + len2
end
end
end
end
# --------------------------------------
function hopping()
boundary(bry) # Start the lattice indeces
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
ham[i,i] = m0 - mu + u0*((m1-yp)^2+(m2-xp)^2)
ham[len2+i,len2+i] = - m0 - mu + u0*((m1-yp)^2+(m2-xp)^2)
ham[len2*2+i,len2*2+i] = m0 - mu + u0*((m1-yp)^2+(m2-xp)^2)
ham[len2*3+i,len2*3+i] = - m0 - mu + u0*((m1-yp)^2+(m2-xp)^2)
ham[len2*4+i,len2*4+i] = - m0 + mu - u0*((m1-yp)^2+(m2-xp)^2)
ham[len2*5+i,len2*5+i] = m0 + mu - u0*((m1-yp)^2+(m2-xp)^2)
ham[len2*6+i,len2*6+i] = - m0 + mu - u0*((m1-yp)^2+(m2-xp)^2)
ham[len2*7+i,len2*7+i] = m0 + mu - u0*((m1-yp)^2+(m2-xp)^2)
end
end
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
ham[i,bry[1,i]] = -tx/2*phi(m1)
ham[i,bry[2,i]] = -tx/2*conj(phi(m1))
ham[i,bry[3,i]] = -ty/2
ham[i,bry[4,i]] = -ty/2
if m1 == 1
ham[i,bry[4,i]] = -ty/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[i,bry[3,i]] = -ty/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[i,bry[2,i]] = -tx/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[i,bry[1,i]] = -tx/2*phi(m1)*bcx
end
#--------------
ham[len2 + i,len2 + bry[1,i]] = tx/2*phi(m1)
ham[len2 + i,len2 + bry[2,i]] = tx/2*conj(phi(m1))
ham[len2 + i,len2 + bry[3,i]] = ty/2
ham[len2 + i,len2 + bry[4,i]] = ty/2
if m1 == 1
ham[len2 + i,len2 + bry[4,i]] = ty/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[len2 + i,len2 + bry[3,i]] = ty/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[len2 + i,len2 + bry[2,i]] = tx/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[len2 + i,len2 + bry[1,i]] = tx/2*phi(m1)*bcx
end
#----------------------------------------------------
ham[len2*2 + i,len2*2 + bry[1,i]] = -tx/2*phi(m1)
ham[len2*2 + i,len2*2 + bry[2,i]] = -tx/2*conj(phi(m1))
ham[len2*2 + i,len2*2 + bry[3,i]] = -ty/2
ham[len2*2 + i,len2*2 + bry[4,i]] = -ty/2
if m1 == 1
ham[len2*2 + i,len2*2 + bry[4,i]] = -ty/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[len2*2 + i,len2*2 + bry[3,i]] = -ty/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[len2*2 + i,len2*2 + bry[2,i]] = -tx/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[len2*2 + i,len2*2 + bry[1,i]] = -tx/2*phi(m1)*bcx
end
#--------------------------------------------------------
ham[len2*3+i,len2*3+bry[1,i]] = tx/2*phi(m1)
ham[len2*3+i,len2*3+bry[2,i]] = tx/2*conj(phi(m1))
ham[len2*3+i,len2*3+bry[3,i]] = ty/2
ham[len2*3+i,len2*3+bry[4,i]] = ty/2
if m1 == 1
ham[len2*3 + i,len2*3 + bry[4,i]] = ty/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[len2*3 + i,len2*3 + bry[3,i]] = ty/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[len2*3 + i,len2*3 + bry[2,i]] = tx/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[len2*3 + i,len2*3 + bry[1,i]] = tx/2*phi(m1)*bcx
end
#-----------------
ham[len2*4+i,len2*4+bry[1,i]] = tx/2*conj(phi(m1))
ham[len2*4+i,len2*4+bry[2,i]] = tx/2*phi(m1)
ham[len2*4+i,len2*4+bry[3,i]] = ty/2
ham[len2*4+i,len2*4+bry[4,i]] = ty/2
if m1 == 1
ham[len2*4 + i,len2*4 + bry[4,i]] = ty/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*4 + i,len2*4 + bry[3,i]] = ty/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*4 + i,len2*4 + bry[2,i]] = tx/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*4 + i,len2*4 + bry[1,i]] = tx/2*conj(phi(m1))*bcx
end
#-------------------
ham[len2*5+i,len2*5+bry[1,i]] = -tx/2*conj(phi(m1))
ham[len2*5+i,len2*5+bry[2,i]] = -tx/2*phi(m1)
ham[len2*5+i,len2*5+bry[3,i]] = -ty/2
ham[len2*5+i,len2*5+bry[4,i]] = -ty/2
if m1 == 1
ham[len2*5 + i,len2*5 + bry[4,i]] = -ty/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*5 + i,len2*5 + bry[3,i]] = -ty/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*5 + i,len2*5 + bry[2,i]] = -tx/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*5 + i,len2*5 + bry[1,i]] = -tx/2*conj(phi(m1))*bcx
end
#--------------------
ham[len2*6+i,len2*6+bry[1,i]] = tx/2*conj(phi(m1))
ham[len2*6+i,len2*6+bry[2,i]] = tx/2*phi(m1)
ham[len2*6+i,len2*6+bry[3,i]] = ty/2
ham[len2*6+i,len2*6+bry[4,i]] = ty/2
if m1 == 1
ham[len2*6 + i,len2*6 + bry[4,i]] = ty/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*6 + i,len2*6 + bry[3,i]] = ty/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*6 + i,len2*6 + bry[2,i]] = tx/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*6 + i,len2*6 + bry[1,i]] = tx/2*conj(phi(m1))*bcx
end
#----------------------
ham[len2*7+i,len2*7+bry[1,i]] = -tx/2*conj(phi(m1))
ham[len2*7+i,len2*7+bry[2,i]] = -tx/2*phi(m1)
ham[len2*7+i,len2*7+bry[3,i]] = -ty/2
ham[len2*7+i,len2*7+bry[4,i]] = -ty/2
if m1 == 1
ham[len2*7 + i,len2*7 + bry[4,i]] = -ty/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*7 + i,len2*7 + bry[3,i]] = -ty/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*7 + i,len2*7 + bry[2,i]] = -tx/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*7 + i,len2*7 + bry[1,i]] = -tx/2*conj(phi(m1))*bcx
end
end
end
end
# ===============================================================================
function couple()
boundary(bry)
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
ham[i,len2 + bry[1,i]] = -im*ax/2*phi(m1)
ham[i,len2 + bry[2,i]] = im*ax/2*conj(phi(m1))
ham[i,len2 + bry[3,i]] = -ay/2
ham[i,len2 + bry[4,i]] = ay/2
if m1 == 1
ham[i,len2 + bry[4,i]] = ay/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[i,len2 + bry[3,i]] = -ay/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[i,len2 + bry[2,i]] = im*ax/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[i,len2 + bry[1,i]] = -im*ax/2*phi(m1)*bcx
end
#---------------------------------------------
ham[len2 + i, bry[1,i]] = -im*ax/2*phi(m1)
ham[len2 + i,bry[2,i]] = im*ax/2*conj(phi(m1))
ham[len2 + i,bry[3,i]] = ay/2
ham[len2 + i,bry[4,i]] = -ay/2
if m1 == 1
ham[len2 + i,bry[4,i]] = -ay/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[len2 + i,bry[3,i]] = ay/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[len2 + i,bry[2,i]] = im*ax/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[len2 + i,bry[1,i]] = -im*ax/2*phi(m1)*bcx
end
#---------------------------------------------
ham[len2*2 + i,len2*3 + bry[1,i]] = im*ax/2*phi(m1)
ham[len2*2 + i,len2*3 + bry[2,i]] = -im*ax/2*conj(phi(m1))
ham[len2*2 + i,len2*3 + bry[3,i]] = -ay/2
ham[len2*2 + i,len2*3 + bry[4,i]] = ay/2
if m1 == 1
ham[len2*2 + i,len2*3 + bry[4,i]] = ay/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[len2*2 + i,len2*3 + bry[3,i]] = -ay/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[len2*2 + i,len2*3 + bry[2,i]] = -im*ax/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[len2*2 + i,len2*3 + bry[1,i]] = im*ax/2*phi(m1)*bcx
end
#---------------------------------------------
ham[len2*3 + i,len2*2 + bry[1,i]] = im*ax/2*phi(m1)
ham[len2*3 + i,len2*2 + bry[2,i]] = -im*ax/2*conj(phi(m1))
ham[len2*3 + i,len2*2 + bry[3,i]] = ay/2
ham[len2*3 + i,len2*2 + bry[4,i]] = -ay/2
if m1 == 1
ham[len2*3 + i,len2*2 + bry[4,i]] = -ay/2*phi(m2*yn)*bcy
end
if m1 == yn
ham[len2*3 + i,len2*2 + bry[3,i]] = ay/2*conj(phi(m2*yn))*bcy
end
if m2 == 1
ham[len2*3 + i,len2*2 + bry[2,i]] = -im*ax/2*conj(phi(m1))*bcx
end
if m2 == xn
ham[len2*3 + i,len2*2 + bry[1,i]] = im*ax/2*phi(m1)*bcx
end
#---------------------------------------------
ham[len2*4 + i,len2*5 + bry[1,i]] = -im*ax/2*conj(phi(m1))
ham[len2*4 + i,len2*5 + bry[2,i]] = im*ax/2*phi(m1)
ham[len2*4 + i,len2*5 + bry[3,i]] = ay/2
ham[len2*4 + i,len2*5 + bry[4,i]] = -ay/2
if m1 == 1
ham[len2*4 + i,len2*5 + bry[4,i]] = -ay/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*4 + i,len2*5 + bry[3,i]] = ay/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*4 + i,len2*5 + bry[2,i]] = im*ax/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*4 + i,len2*5 + bry[1,i]] = -im*ax/2*conj(phi(m1))*bcx
end
#---------------------------------------------
ham[len2*5 + i,len2*4 + bry[1,i]] = -im*ax/2*conj(phi(m1))
ham[len2*5 + i,len2*4 + bry[2,i]] = im*ax/2*phi(m1)
ham[len2*5 + i,len2*4 + bry[3,i]] = -ay/2
ham[len2*5 + i,len2*4 + bry[4,i]] = ay/2
if m1 == 1
ham[len2*5 + i,len2*4 + bry[4,i]] = ay/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*5 + i,len2*4 + bry[3,i]] = -ay/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*5 + i,len2*4 + bry[2,i]] = im*ax/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*5 + i,len2*4 + bry[1,i]] = -im*ax/2*conj(phi(m1))*bcx
end
#---------------------------------------------
ham[len2*6 + i,len2*7 + bry[1,i]] = im*ax/2*conj(phi(m1))
ham[len2*6 + i,len2*7 + bry[2,i]] = -im*ax/2*phi(m1)
ham[len2*6 + i,len2*7 + bry[3,i]] = ay/2
ham[len2*6 + i,len2*7 + bry[4,i]] = -ay/2
if m1 == 1
ham[len2*6 + i,len2*7 + bry[4,i]] = -ay/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*6 + i,len2*7 + bry[3,i]] = ay/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*6 + i,len2*7 + bry[2,i]] = -im*ax/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*6 + i,len2*7 + bry[1,i]] = im*ax/2*conj(phi(m1))*bcx
end
#---------------------------------------------
ham[len2*7 + i,len2*6 + bry[1,i]] = im*ax/2*conj(phi(m1))
ham[len2*7 + i,len2*6 + bry[2,i]] = -im*ax/2*phi(m1)
ham[len2*7 + i,len2*6 + bry[3,i]] = -ay/2
ham[len2*7 + i,len2*6 + bry[4,i]] = ay/2
if m1 == 1
ham[len2*7 + i,len2*6 + bry[4,i]] = ay/2*conj(phi(m2*yn))*bcy
end
if m1 == yn
ham[len2*7 + i,len2*6 + bry[3,i]] = -ay/2*phi(m2*yn)*bcy
end
if m2 == 1
ham[len2*7 + i,len2*6 + bry[2,i]] = -im*ax/2*phi(m1)*bcx
end
if m2 == xn
ham[len2*7 + i,len2*6 + bry[1,i]] = im*ax/2*conj(phi(m1))*bcx
end
end
end
end
# ====================================================================================
function scpair()
boundary(bry) # Setting boundary conditions
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
ham[i,len2*6 + i] = -del0/2
ham[len2 + i,len2*7 + i] = -del0/2
ham[len2*2 + i,len2*4 + i] = del0/2
ham[len2*3 + i,len2*5 + i] = del0/2
ham[len2*6 + i,i] = -del0/2
ham[len2*7 + i,len2 + i] = -del0/2
ham[len2*4 + i,len2*2 + i] = del0/2
ham[len2*5 + i,len2*3 + i] = del0/2
end
end
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
ham[i,len2*6 + bry[1,i]] = -delx/2
ham[i,len2*6 + bry[2,i]] = -delx/2
ham[i,len2*6 + bry[3,i]] = -dely/2
ham[i,len2*6 + bry[4,i]] = -dely/2
if m1 == 1
ham[i,len2*6 + bry[4,i]] = -dely/2*bcy
end
if m1 == yn
ham[i,len2*6 + bry[3,i]] = -dely/2*bcy
end
if m2 == 1
ham[i,len2*6 + bry[2,i]] = -delx/2*bcx
end
if m2 == xn
ham[i,len2*6 + bry[1,i]] = -delx/2*bcx
end
#----------------------------------
ham[len2 + i,len2*7 + bry[1,i]] = -delx/2
ham[len2 + i,len2*7 + bry[2,i]] = -delx/2
ham[len2 + i,len2*7 + bry[3,i]] = -dely/2
ham[len2 + i,len2*7 + bry[4,i]] = -dely/2
if m1 == 1
ham[len2 + i,len2*7 + bry[4,i]] = -dely/2*bcy
end
if m1 == yn
ham[len2 + i,len2*7 + bry[3,i]] = -dely/2*bcy
end
if m2 == 1
ham[len2 + i,len2*7 + bry[2,i]] = -delx/2*bcx
end
if m2 == xn
ham[len2 + i,len2*7 + bry[1,i]] = -delx/2*bcx
end
#---------------------------------------------
ham[len2*2 + i,len2*4 + bry[1,i]] = delx/2
ham[len2*2 + i,len2*4 + bry[2,i]] = delx/2
ham[len2*2 + i,len2*4 + bry[3,i]] = dely/2
ham[len2*2 + i,len2*4 + bry[4,i]] = dely/2
if m1 == 1
ham[len2*2 + i,len2*4 + bry[4,i]] = dely/2*bcy
end
if m1 == yn
ham[len2*2 + i,len2*4 + bry[3,i]] = dely/2*bcy
end
if m2 == 1
ham[len2*2 + i,len2*4 + bry[2,i]] = delx/2*bcx
end
if m2 == xn
ham[len2*2 + i,len2*4 + bry[1,i]] = delx/2*bcx
end
#----------------------------------------
ham[len2*3 + i,len2*5 + bry[1,i]] = delx/2
ham[len2*3 + i,len2*5 + bry[2,i]] = delx/2
ham[len2*3 + i,len2*5 + bry[3,i]] = dely/2
ham[len2*3 + i,len2*5 + bry[4,i]] = dely/2
if m1 == 1
ham[len2*3 + i,len2*5 + bry[4,i]] = dely/2*bcy
end
if m1 == yn
ham[len2*3 + i,len2*5 + bry[3,i]] = dely/2*bcy
end
if m2 == 1
ham[len2*3 + i,len2*5 + bry[2,i]] = delx/2*bcx
end
if m2 == xn
ham[len2*3 + i,len2*5 + bry[1,i]] = delx/2*bcx
end
#------------------------------------------
ham[len2*6 + i,bry[1,i]] = -delx/2
ham[len2*6 + i,bry[2,i]] = -delx/2
ham[len2*6 + i,bry[3,i]] = -dely/2
ham[len2*6 + i,bry[4,i]] = -dely/2
if m1 == 1
ham[len2*6 + i,bry[4,i]] = -dely/2*bcy
end
if m1 == yn
ham[len2*6 + i,bry[3,i]] = -dely/2*bcy
end
if m2 == 1
ham[len2*6 + i,bry[2,i]] = -delx/2*bcx
end
if m2 == xn
ham[len2*6 + i,bry[1,i]] = -delx/2*bcx
end
#---------------------------------------
ham[len2*7 + i,len2 + bry[1,i]] = -delx/2
ham[len2*7 + i,len2 + bry[2,i]] = -delx/2
ham[len2*7 + i,len2 + bry[3,i]] = -dely/2
ham[len2*7 + i,len2 + bry[4,i]] = -dely/2
if m1 == 1
ham[len2*7 + i,len2 + bry[4,i]] = -dely/2*bcy
end
if m1 == yn
ham[len2*7 + i,len2 + bry[3,i]] = -dely/2*bcy
end
if m2 == 1
ham[len2*7 + i,len2 + bry[2,i]] = -delx/2*bcx
end
if m2 == xn
ham[len2*7 + i,len2 + bry[1,i]] = -delx/2*bcx
end
#---------------------------------------
ham[len2*4 + i,len2*2 + bry[1,i]] = delx/2
ham[len2*4 + i,len2*2 + bry[2,i]] = delx/2
ham[len2*4 + i,len2*2 + bry[3,i]] = dely/2
ham[len2*4 + i,len2*2 + bry[4,i]] = dely/2
if m1 == 1
ham[len2*4 + i,len2*2 + bry[4,i]] = dely/2*bcy
end
if m1 == yn
ham[len2*4 + i,len2*2 + bry[3,i]] = dely/2*bcy
end
if m2 == 1
ham[len2*4 + i,len2*2 + bry[2,i]] = delx/2*bcx
end
if m2 == xn
ham[len2*4 + i,len2*2 + bry[1,i]] = delx/2*bcx
end
#-----------------------------------------
ham[len2*5 + i,len2*3 + bry[1,i]] = delx/2
ham[len2*5 + i,len2*3 + bry[2,i]] = delx/2
ham[len2*5 + i,len2*3 + bry[3,i]] = dely/2
ham[len2*5 + i,len2*3 + bry[4,i]] = dely/2
if m1 == 1
ham[len2*5 + i,len2*3 + bry[4,i]] = dely/2*bcy
end
if m1 == yn
ham[len2*5 + i,len2*3 + bry[3,i]] = dely/2*bcy
end
if m2 == 1
ham[len2*5 + i,len2*3 + bry[2,i]] = delx/2*bcx
end
if m2 == xn
ham[len2*5 + i,len2*3 + bry[1,i]] = delx/2*bcx
end
end
end
end
# ===================================================
function delta(x)
gamma = 0.005
return 1/pi*gamma/(x^2+gamma^2)
end
# ===================================================
function ldos(val,vec)
s = zeros(Float32,len2)
fileio = open("ldos.dat","w")
for m1 in 1:yn
for m2=1:yn
i = (m1-1)*xn + m2
for kk = 1:N
s[i] = s[i] + delta(val[i])*(abs(vec[i,kk])^2)+delta(-val[kk])*(abs(vec[i+len2*4,kk])^2);
end
writedlm(fileio,[m1 m2 s[m1]])
end
end
close(fileio)
return s
end
# ==========================================================
function opc(x,y,val,vec,input)
beta = 10^5
s = 0
if input == 1
for kk in 1:N
s = s + (vec[x,kk]*conj(vec[len2*6+y,kk])+vec[y,kk]*conj(vec[len2*6+x,kk]))*tanh(val[kk]/2*beta)
end
elseif input == 2
for kk in 1:N
s = s + (vec[len2+x,kk]*conj(vec[len2*7+y,kk])+vec[len2+y,kk]*conj(vec[len2*7+x,kk]))*tanh(val[kk]/2*beta)
end
end
return vp/4*s
end
# ===============================================================
function order(val,vec)
boundary(bry)
absop = zeros(Float32,len2)
fileio = open("op.dat","w")
for i in 1:len2
d1[i,bry[1,i]] = opc(i,bry[1,i],val,vec,1)
d1[i,bry[2,i]] = opc(i,bry[2,i],val,vec,1)
d1[i,bry[3,i]] = opc(i,bry[3,i],val,vec,1)
d1[i,bry[4,i]] = opc(i,bry[4,i],val,vec,1)
#-------------------------------------------
d2[i,bry[1,i]] = -d1[bry[1,i],i]
d2[i,bry[2,i]] = -d1[bry[2,i],i]
d2[i,bry[3,i]] = -d1[bry[3,i],i]
d2[i,bry[4,i]] = -d1[bry[4,i],i]
#-------------------------------------------
d3[i,bry[1,i]] = opc(i,bry[1,i],val,vec,2)
d3[i,bry[2,i]] = opc(i,bry[2,i],val,vec,2)
d3[i,bry[3,i]] = opc(i,bry[3,i],val,vec,2)
d3[i,bry[4,i]] = opc(i,bry[4,i],val,vec,2)
#------------------------------------------
d4[i,bry[1,i]] = -d3[bry[1,i],i]
d4[i,bry[2,i]] = -d3[bry[2,i],i]
d4[i,bry[3,i]] = -d3[bry[3,i],i]
d4[i,bry[4,i]] = -d3[bry[4,i],i]
end
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
a1 = d1[i,bry[1,i]]*phi(m1)
a2 = d1[i,bry[2,i]]*conj(phi(m1))
a3 = d1[i,bry[3,i]]
a4 = d1[i,bry[4,i]]
# -------------------------
b1 = d3[i,bry[1,i]]*phi(m1)
b2 = d3[i,bry[2,i]]*conj(phi(m1))
b3 = d3[i,bry[3,i]]
b4 = d3[i,bry[4,i]]
#------------------------------
if m1 == yn
a3 = a3*conj(phi(m2*yn))
b3 = b3*conj(phi(m2*yn))
end
if m1 == 1
a4 = a4*phi(m2*yn)
b4 = b4*phi(m2*yn)
end
sop[m1,m2] = (a1 + a2 - a3 - a4 + b1 + b2 - b3 - b4)/8.0
absop[i] = abs(sop[m1,m2])^2
writedlm(fileio,[m1 m2 absop[i]])
end
end
close(fileio)
return absop
end
# ================================================================
function checkElementary()
for m1 in 1:N
for m2 in 1:N
if ham[m1,m2] != conj(ham[m2,m1])
println(m1,m2)
return false
end
end
end
end
# =======================================================
function pairCon()
for i in 1:len2
ham[i,len2*6 + bry[1,i]] = d1[i,bry[1,i]]
ham[i,len2*6 + bry[2,i]] = d1[i,bry[2,i]]
ham[i,len2*6 + bry[3,i]] = d1[i,bry[3,i]]
ham[i,len2*6 + bry[4,i]] = d1[i,bry[4,i]]
#----------------------------------------
ham[len2*2 + i,len2*4 + bry[1,i]] = d2[i,bry[1,i]]
ham[len2*2 + i,len2*4 + bry[2,i]] = d2[i,bry[2,i]]
ham[len2*2 + i,len2*4 + bry[3,i]] = d2[i,bry[3,i]]
ham[len2*2 + i,len2*4 + bry[4,i]] = d2[i,bry[4,i]]
#----------------------------------------
ham[len2 + i,len2*7 + bry[1,i]] = d3[i,bry[1,i]]
ham[len2 + i,len2*7 + bry[2,i]] = d3[i,bry[2,i]]
ham[len2 + i,len2*7 + bry[3,i]] = d3[i,bry[3,i]]
ham[len2 + i,len2*7 + bry[4,i]] = d3[i,bry[4,i]]
#----------------------------------------
ham[len2*3 + i,len2*5 + bry[1,i]] = d4[i,bry[1,i]]
ham[len2*3 + i,len2*5 + bry[2,i]] = d4[i,bry[2,i]]
ham[len2*3 + i,len2*5 + bry[3,i]] = d4[i,bry[3,i]]
ham[len2*3 + i,len2*5 + bry[4,i]] = d4[i,bry[4,i]]
end
#------------------------------------------------------
for m1 in 1:len2
for m2 in 1:len2
ham[len2*6 + m1,m2] = conj(ham[m2,len2*6 + m1])
ham[len2*7 + m1,len2 + m2] = conj(ham[len2 + m2,len2*7 + m1])
ham[len2*4 + m1,len2*2 + m2] = conj(ham[len2*2 + m2,len2*4 + m1])
ham[len2*5 + m1,len2*3 + m2] = conj(ham[len2*3 + m2,len2*5 + m1])
end
end
end
# =========================================================
function matrixSet(input)
hopping()
couple()
if input == 0
scpair()
else
pairCon()
end
matrixCheck()
end
# ==================================================
function fermi(x)
beta = 10^5
return 1.0/(exp(x*beta) + 1.0)
end
# ======================================================
function MagNum(val,vec)
fileio = open("magnum.dat","w")
s = zeros(Float32,xn,yn)
for m1 in 1:yn
for m2 in 1:xn
i = (m1-1)*xn + m2
s1, s2 = 0, 0
for m3 in 1:N
s1 = s1 + (abs(vec[i,m3])^2 + abs(vec[len2 + i,m3])^2)*fermi(val[m3])
s2 = s2 + (abs(vec[len2*6 + i,m3])^2 + abs(vec[len2*7 + i,m3])^2)*(1.0-fermi(val[m3]))
end
s[m1,m2] = 0.5*(s1-s2)
writedlm(fileio,[m1 m2 0.5*(s1-s2)])
end
end
return s
close(fileio)
end
# ===================================================================
function ParticleNumber(val,vec)
fileio = open("particlenumber.dat","w")
s = zeros(Float32,xn,yn)
for m1 in 1:yn
for m2 in 1:yn
i = (m1-1)*xn + m2
s1, s2 = 0,0
for m3 in 1:N
s1 = s1 + (abs(vec[i,m3])^2 + abs(vec[len2 + i,m3])^2)*fermi(val[m3])
s2 = s2 + (abs(vec[len2*6 + i,m3])^2 + abs(vec[len2*7 + i,m3])^2)*(1.0-fermi(val[m3]))
end
s[m1,m2] = s1 + s2
writedlm(fileio,[m1 m2 s1+s2])
end
end
return s
close(fileio)
end
# ===========================================================================
function loop(val,vec)
del_loop = zeros(ComplexF32,xn,yn)
del_err = zeros(ComplexF32,xn,yn)
ref = 3
num = 0
refArray = Int64[]
fileio = open("check.dat","w")
order(val,vec)
while ref > 0
ref = 0
num += 1
for m1 in 1:yn
for m2 in 1:xn
del_err[m1,m2] = sop[m1,m2] - del_loop[m1,m2]
del_loop[m1,m2] = sop[m1,m2]
if abs(real(del_err[m1,m2])) > err
ref += 1
end
if abs(imag(del_err[m1,m2])) > err
ref += 1
end
end
end
writedlm(fileio,ref)
matrixSet(10)
if ham' == ham
eigval,eigvec = eigen(ham)
fileio = open("eigenval.dat","w")
for m1 in 1:N
writedlm(fileio,eigval[m1])
end
close(fileio)
order(eigval,eigvec)
ldos(eigval,eigvec)
ParticleNumber(eigval,eigvec)
MagNum(eigval,eigvec)
end
end
end
# ==================================================
function matrixCheck()
for m1 in 1:N
for m2 in 1:N
if ham[m1,m2] != conj(ham[m2,m1])
println("(",m1,",",m2,")",ham[m1,m2],ham[m2,m1])
break
end
end
end
end
# =======================================================
function imageShow(val,vec,scale)
figure(figsize=(16,5))
subplot(121)
PyPlot.plot(val,"b.") # visual the eigenval
xlim(length(val)/2-40,length(val)/2+40)
ylim(-delx,delx)
subplot(122)
ld = ldos(eigval,eigvec)
x = Int32[]
y = Int32[]
for m1 in 1:yn
for m2 in 1:xn
append!(x,m2)
append!(y,m1)
end
end
scatter(x,y,ld*scale,c=ld,edgecolors="b",cmap="Reds")
colorbar()
PyPlot.savefig("result.png",bbox_inches="tight",dpi=500)
end
主函数
# Run program
using CPUTime
using PyPlot,LinearAlgebra
using DelimitedFiles
# parameter setting
xn = 30
yn = 30
N = xn*yn*8
len1 = xn
len2 = xn*yn
# Boundary condition setting
bcx = 0 # x boundary (zero is open boundary,one is periodic boundary)
bcy = 0 # y boundary (zero is open boundary,one is periodic boundary)
T = 1e-4
Kb = 1
beta = 1/(Kb*T)
err = 10^(-6)
# System parameter
m0 = 1.5
tx = 1.0
ty = 1.0
mu = 0.0
V = 0 # point potential
ax = 1.0
ay = 1.0
# Spuerconduct pair
del0 = 0
delx = 0.5
dely = -0.5
vp = 5.0
# magnetic paramater
phi0 = pi
#B = 2.0*phi0/(xn*yn)
B = 0
# -----------------------------------
u0 = 0
if u0 != 0
vp = 2.0
xp = ceil(xn/2)
yp = ceil(yn/2)
B = 0
else
xp = 0
yp = 0
end
# ---------------------------------------------
bry = zeros(Int32,4,len2) # store boundary indices
ham = zeros(ComplexF32,N,N) # Hamiltonian initional
op = zeros(Float32,len2)
# ---------------------------------------
# self consistance variables
d1 = zeros(ComplexF32,len2,len2)
d2 = zeros(ComplexF32,len2,len2)
d3 = zeros(ComplexF32,len2,len2)
d4 = zeros(ComplexF32,len2,len2)
sop = zeros(ComplexF32,xn,yn)
# -----------------------------
matrixSet(0)
if ham' == ham # Check the Hamiltonian is Hermitian matrix
#@time @CPUtime eigval,eigvec = eigen(ham)
eigval,eigvec = eigen(ham)
fileio = open("eigenval.dat","w")
for m1 in 1:N
writedlm(fileio,eigval[m1])
end
close(fileio)
end
结果展示
#loop(eigval,eigvec)
imageShow(eigval,eigvec,100)
边界态计算
import os
import numpy as np
from math import *
import matplotlib.pyplot as plt
#-----------------------------------
def Pauli():
hn = 2
i0 = np.zeros([hn,hn],complex)
x = np.zeros([hn,hn],complex)
y = np.zeros([hn,hn],complex)
z = np.zeros([hn,hn],complex)
i0[0,0] = 1
i0[1,1] = 1
x[0,1] = 1
x[1,0] = 1
y[0,1] = -1j
y[1,0] = 1j
z[0,0] = 1
z[1,1] = -1
g1 = np.kron(np.kron(z,i0),z)
g2 = np.kron(np.kron(x,z),i0)
g3 = np.kron(np.kron(y,i0),z)
g4 = np.kron(np.kron(i0,y),y)
g5 = np.kron(np.kron(i0,x),z)
return g1,g2,g3,g4,g5
#----------------------------------
def openy(h0, kx, yn):
hn = 8
N = yn*hn
m0 = 1.
tx = 2.0
ty = 2.0
ax = 2.0
ay = 2.0
dx = 0.
dy = -dx
d0 = 0.4
g1,g2,g3,g4,g5 = Pauli()
ham = np.zeros([N,N],dtype=complex)
for k in range(0,yn):
if k==0:
for i1 in range(0,hn):
for i2 in range(0,hn):
ham[i1,i2] = (m0-tx*cos(kx))*g1[i1,i2] + ax*sin(kx)*g2[i1,i2] + (d0 + dx*cos(kx))*g4[i1,i2]+ h0*g5[i1,i2]
ham[i1,i2 + hn] = (-ty*g1[i1,i2] - 1j*ay*g3[i1,i2])/2 + dy/2.0*g4[i1,i2]
elif k == yn - 1:
for i1 in range(0,hn):
for i2 in range(0,hn):
ham[k*hn + i1,k*hn + i2] = (m0-tx*cos(kx))*g1[i1,i2] + ax*sin(kx)*g2[i1,i2] + (d0 + dx*cos(kx))*g4[i1,i2]+ h0*g5[i1,i2]
ham[k*hn + i1,k*hn + i2 - hn] = -ty*g1[i1,i2]/2 + 1j*ay*g3[i1,i2]/2 + dy/2.0*g4[i1,i2]
else:
for i1 in range(0,hn):
for i2 in range(0,hn):
ham[k*hn + i1,k*hn + i2] = (m0-tx*cos(kx))*g1[i1,i2] + ax*sin(kx)*g2[i1,i2] + (d0 + dx*cos(kx))*g4[i1,i2]+ h0*g5[i1,i2]
ham[k*hn + i1,k*hn + i2 + hn] = (-ty*g1[i1,i2] - 1j*ay*g3[i1,i2] )/2 + dy/2.0*g4[i1,i2]
ham[k*hn + i1,k*hn + i2 - hn] = -ty*g1[i1,i2]/2 + 1j*ay*g3[i1,i2]/2 + dy/2.0*g4[i1,i2]
return ham
#----------------------------------
def openx(h0, ky, yn):
hn = 8
N = yn*hn
m0 = 1.
tx = 2.0
ty = 2.0
ax = 2.0
ay = 2.0
dx = 0.
dy = -dx
d0 = 0.4
g1,g2,g3,g4,g5 = Pauli()
ham = np.zeros([N,N],dtype = complex)
for k in range(0,yn):
if k == 0:
for i1 in range(0,hn):
for i2 in range(0,hn):
ham[i1,i2] = (m0 - ty*cos(ky))*g1[i1,i2] + ay*sin(ky)*g3[i1,i2] + (d0 + dy*cos(ky))*g4[i1,i2]+ h0*g5[i1,i2]
ham[i1,i2 + hn] = (-tx*g1[i1,i2] - 1j*ax*g2[i1,i2])/2 + dx/2.0*g4[i1,i2]
elif k == yn - 1:
for i1 in range(0,hn):
for i2 in range(0,hn):
ham[k*hn + i1,k*hn + i2] = (m0 - ty*cos(ky))*g1[i1,i2] + ay*sin(ky)*g3[i1,i2] + (d0 + dy*cos(ky))*g4[i1,i2]+ h0*g5[i1,i2]
ham[k*hn + i1,k*hn + i2 - hn] = (-tx*g1[i1,i2] + 1j*ax*g2[i1,i2])/2 + dx/2.0*g4[i1,i2]
else:
for i1 in range(0,hn):
for i2 in range(0,hn):
ham[k*hn + i1,k*hn + i2] = (m0 - ty*cos(ky))*g1[i1,i2] + ay*sin(ky)*g3[i1,i2] + (d0 + dy*cos(ky))*g4[i1,i2]+ h0*g5[i1,i2]
ham[k*hn + i1,k*hn + i2 + hn] = (-tx*g1[i1,i2] - 1j*ax*g2[i1,i2])/2 + dx/2.0*g4[i1,i2]
ham[k*hn + i1,k*hn + i2 - hn] = (-tx*g1[i1,i2] + 1j*ax*g2[i1,i2])/2 + dx/2.0*g4[i1,i2]
return ham
#------------------------------------------
def cylinder():
h0 = 0.
hn = 8
yn = 50
N = hn*yn
kn = 50
ham = np.zeros([N,N],dtype=complex)
klist = []
vallist1 = []
vallist2 = []
for i1 in range(-kn,kn):
kx = i1*pi/kn
klist.append(kx/pi)
ham = openx(h0,kx,yn)
val,vec = np.linalg.eigh(ham)
vallist1.append(val)
ham = openy(h0,kx,yn)
val,vec = np.linalg.eigh(ham)
vallist2.append(val)
return klist,vallist1,vallist2
#-------------------------------------------
def main():
# Wilsonloop()
f1name = "cylinder.png"
r1,r2,r3 = cylinder()
fig = plt.figure(figsize=(14, 7))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.plot(r1,r2,color = 'blue')
ax2.plot(r1,r3,color = 'blue')
plt.show()
plt.savefig(f1name, dpi = 600,bbox_inches = 'tight')
# print(r1)
#---------------------------------
def main():
g1,g2,g3,g4,g5 = Pauli()
print(g1)
#----------------------------------
if __name__=="__main__":
os.chdir(os.getcwd())
ham = main()
公众号
相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
