利用格林函数求解边界态-V2

这篇博客中利用另外一种格林函数的方法来计算拓扑绝缘体的边界态.
{:.info}

在前面利用格林函数求解边界态这篇博客中利用迭代格林函数的方法计算了拓扑绝缘体的边界态,但是结果中边界态的确很明显,但是体态的性质却不会表现的很清晰,这里就像利用矩阵求逆的方式,直接计算哈密顿量对应的谱函数,从而来计算边界态和体态的性质,但是这个方法的缺点也是比较明显的,那就是耗时比较长.

模型

仍然以拓扑绝缘体的哈密顿量为出发点

$$H(\mathbf{k})=(m_0-t_x\cos k_x-t_y\cos k_y)\sigma_z+\lambda_x\sin k_x\sigma_xs_z+\lambda_y\sin k_y\sigma_y\label{ham}$$

对这个哈密顿量沿着$x$方向取开边界,沿$y$方向仍然是周期边界条件,这样就可以在一个圆柱结构上计算边界态了,对应的谱函数为

$$A(k,\omega)=-2\text{Im}[\omega+i\epsilon-H(\mathcal{k})]^{-1}$$

这里的$H(\mathbf{k})$是开边界之后的哈密顿量.

代码

    module pub
    implicit none
    integer yn,kn,ne,N,enn,hn
    real eta,dk,de
    parameter(yn = 50,hn =  4,kn = 50, ne = 50,N = yn*4,eta = 0.01,dk = 0.01,de = dk)
    real,parameter::pi = 3.1415926535
    complex,parameter::im = (0.,1.0)  
    complex Ham(N,N),one(N,N) 
!=================================
    real m0,mu   
    real tx,ty
    real ax,ay,gamma
    complex g1(hn,hn) 
    complex g2(hn,hn) 
    complex g3(hn,hn) 
!================cheev===============
    integer::lda = N
    integer,parameter::lwmax = 2*N + N**2
    real,allocatable::w(:)
    complex,allocatable::work(:)
    real,allocatable::rwork(:)
    integer,allocatable::iwork(:)
    integer lwork   
    integer lrwork    
    integer liwork   
    integer info
    end module pub
!================= PROGRAM START ============================
    program sol
    use pub
    integer i1
!====空间申请==================
    allocate(w(N))
    allocate(work(lwmax))
    allocate(rwork(1+5*N+2*N**2))
    allocate(iwork(3+5*N))
!======parameter value setting =====
    m0 = 1.5
    mu = 0
    tx = 1.0
    ty = 1.0
    ax = 1.0
    ay = 1.0
    do i1 = 1,N
        one(i1,i1) = 1  !单位矩阵
    end do
    call main()
    stop 
    end program sol
!=============================================================================================
    subroutine main()
    !  Calculate edge spectrum function
    use pub
    integer m1,m2,i1
    real kx,ky,omega,re2
    complex h1(N,N),h2(N,N)
    complex re1
    open(30,file="openy-bhz.dat")
    !-------------------------------------------------
    !   y-direction is open
    do omega = -3.0,3.0,de
        do kx = -pi,pi,dk
            re1 = 0
            call openy(kx)
            h1 = omega*one - ham + im*eta*one
            call inv(h1,h2)
            do i1 = 1,N
                re1 = re1 + h2(i1,i1)
            end do
            re2 = -2*aimag(re1)/pi/N
            write(30,999)kx/pi,omega,re2
        end do
        write(30,*)" "
    end do
    close(30)
    !---------------------------------------------------------------------
    !  x-direction is open
    open(31,file="openx-bhz.dat")
    do omega = -3.0,3.0,de
        do ky = -pi,pi,dk
            re1 = 0
            call openx(ky)
            h1 = omega*one - ham + im*eta*one
            call inv(h1,h2)
            do i1 = 1,N
                re1 = re1 + h2(i1,i1)
            end do
            re2 = -2*aimag(re1)/pi/N
            write(31,999)ky/pi,omega,re2
        end do
        write(31,*)" "
    end do
    close(31)
999 format(3f11.5)
    return
    end subroutine main
!======================== Pauli Matrix driect product============================
    subroutine Pauli()
    use pub
    !   TI
!=== Kinetic energy
    g1(1,1) = 1
    g1(2,2) = -1
    g1(3,3) = 1
    g1(4,4) = -1
!====== SOC-x
    g2(1,2) = 1
    g2(2,1) = 1
    g2(3,4) = -1
    g2(4,3) = -1
!====== SOC-y
    g3(1,2) = -im
    g3(2,1) = im
    g3(3,4) = -im
    g3(4,3) = im
    return
    end subroutine Pauli
!==========================================================
    subroutine openx(ky)
    use pub
    real ky
    integer m,l,k
    call Pauli()
    Ham = 0
!========== Positive energy ========
    do k = 0,yn-1
        if (k == 0) then ! Only right block in first line
            do m = 1,hn
                do l = 1,hn
                    Ham(m,l) = (m0-ty*cos(ky))*g1(m,l) + ay*sin(ky)*g3(m,l) 

                    Ham(m,l + hn) = (-tx*g1(m,l) - im*ax*g2(m,l))/2
                end do
            end do
        elseif ( k==yn-1 ) then ! Only left block in last line
            do m = 1,hn
                do l = 1,hn
                    Ham(k*hn + m,k*hn + l) = (m0-ty*cos(ky))*g1(m,l) + ay*sin(ky)*g3(m,l)

                    Ham(k*hn + m,k*hn + l - hn) = -tx*g1(m,l)/2 + im*ax*g2(m,l)/2 
                end do
            end do
        else
            do m = 1,hn
                do l = 1,hn ! k start from 1,matrix block from 2th row
                    Ham(k*hn + m,k*hn + l) = (m0 - ty*cos(ky))*g1(m,l) + ay*sin(ky)*g3(m,l)

                    Ham(k*hn + m,k*hn + l + hn) = (-tx*g1(m,l) - im*ax*g2(m,l))/2 
                    Ham(k*hn + m,k*hn + l - hn) = -tx*g1(m,l)/2 + im*ax*g2(m,l)/2 
                end do
            end do
        end if
    end do
    !------------------------
    call isHermitian()
    ! call eigsol()
    return
    end subroutine openx
!==========================================================
    subroutine openy(kx)
    use pub
    real kx
    integer m,l,k
    call Pauli()
    Ham = 0
!========== Positive energy ========
    do k = 0,yn-1
        if (k == 0) then ! Only right block in first line
            do m = 1,hn
                do l = 1,hn
                    Ham(m,l) = (m0-tx*cos(kx))*g1(m,l) + ax*sin(kx)*g2(m,l)

                    Ham(m,l + hn) = (-ty*g1(m,l) - im*ay*g3(m,l))/2
                end do
            end do
        elseif ( k==yn-1 ) then ! Only left block in last line
            do m = 1,hn
                do l = 1,hn
                    Ham(k*hn + m,k*hn + l) = (m0-tx*cos(kx))*g1(m,l) + ax*sin(kx)*g2(m,l)

                    Ham(k*hn + m,k*hn + l - hn) = -ty*g1(m,l)/2 + im*ay*g3(m,l)/2 
                end do
            end do
        else
            do m = 1,hn
                do l = 1,hn ! k start from 1,matrix block from 2th row
                    Ham(k*hn + m,k*hn + l) = (m0-tx*cos(kx))*g1(m,l) + ax*sin(kx)*g2(m,l)

                    Ham(k*hn + m,k*hn + l + hn) = (-ty*g1(m,l) - im*ay*g3(m,l) )/2 
                    Ham(k*hn + m,k*hn + l - hn) = -ty*g1(m,l)/2 + im*ay*g3(m,l)/2 
                end do
            end do
        end if
    end do
    !---------------------------------
    call isHermitian()
    ! call eigsol()
    return
    end subroutine openy
!============================================================
    subroutine isHermitian()
    use pub
    integer i,j
    do i = 1,N
        do j = 1,N
            if (Ham(i,j) .ne. conjg(Ham(j,i)))then
                open(160,file = 'hermitian.dat')
                write(160,*)i,j
                write(160,*)Ham(i,j)
                write(160,*)Ham(j,i)
                write(160,*)"===================="
                write(*,*)"Hamiltonian is not hermitian"
                stop
            end if
        end do
    end do
    close(160)
    return
    end subroutine isHermitian
!================= 矩阵本征值求解 ==============
    subroutine eigSol()
    use pub
    integer m
    lwork = -1
    liwork = -1
    lrwork = -1
    call cheevd('V','Upper',N,Ham,lda,w,work,lwork &
        ,rwork,lrwork,iwork,liwork,info)
    lwork = min(2*N+N**2, int( work( 1 ) ) )
    lrwork = min(1+5*N+2*N**2, int( rwork( 1 ) ) )
    liwork = min(3+5*N, iwork( 1 ) )
    call cheevd('V','Upper',N,Ham,lda,w,work,lwork &
        ,rwork,lrwork,iwork,liwork,info)
    if( info .GT. 0 ) then
        open(110,file="mes.dat",status="unknown")
        write(110,*)'The algorithm failed to compute eigenvalues.'
        close(110)
    end if
    ! open(120,file="eigval.dat")
    ! do m = 1,N
    !     write(120,*)w(m)
    ! end do
    ! close(120)
    return
    end subroutine eigSol
!=================================================================
    subroutine inv(matin,matout)
    use pub
    complex,intent(in) :: matin(N,N)
    complex:: matout(size(matin,1),size(matin,2))
    real:: work2(size(matin,1))            ! work array for LAPACK
    integer::info2,ipiv(size(matin,1))     ! pivot indices
    ! Store A in Ainv to prevent it from being overwritten by LAPACK
    matout = matin
    ! SGETRF computes an LU factorization of a general M-by-N matrix A
    ! using partial pivoting with row interchanges.
    call CGETRF(N,N,matout,N,ipiv,info2)
    if (info2.ne.0)then
      write(*,*)'Matrix is numerically singular!'
      stop
    end if
    ! SGETRI computes the inverse of a matrix using the LU factorization
    ! computed by SGETRF.
    call CGETRI(N,matout,N,ipiv,work2,N,info2)
    if (info2.ne.0)then
        write(*,*)'Matrix inversion failed!'
        stop
    end if
    return
    end subroutine inv

编译命令

ifort -mkl file.f90

绘图

因为想要看到清晰的结果,在离散撒点的时候,间隔必须很小,那么数据文件就会很大,这里推荐gnuplot绘图,

set encoding iso_8859_1
#set terminal  postscript enhanced color
#set output 'arc_r.eps'
#set terminal pngcairo truecolor enhanced  font ",50" size 1920, 1680
set terminal png truecolor enhanced font ",50" size 1920, 1680
set output 'density.png'
#set palette defined ( -10 "#194eff", 0 "white", 10 "red" )
set palette defined ( -10 "blue", 0 "white", 10 "red" )
#set palette rgbformulae 33,13,10
unset ztics
unset key
set pm3d
set border lw 6
set size ratio 1
set view map
set xtics
set ytics
#set xlabel "K_1 (1/{\305})"
set xlabel "X_1"
#set ylabel "K_2 (1/{\305})"
set ylabel "Y"
set ylabel offset 1, 0
set colorbox
set xrange [-1:1]
set yrange [-3:3]
set pm3d interpolate 4,4
#splot 'wavenorm.dat' u 1:2:3 w pm3d
#splot 'wavenorm.dat' u 1:2:3 w pm3d
splot 'openy-bhz.dat' u 1:2:3 w pm3d

绘图

gnuplot file.gnu

最后的结果为

相比较于利用格林函数求解边界态中的结果,可以清晰的看到边界态和体态的结果.如果接的图中的颜色并不能很好的展现,可以在gnuplot绘图中修改颜色方案

set palette defined ( -10 "white", 0 "blue", 10 "red" )

在新的配色方案中,边界态就可以很清晰的表现出来

代码下载

点击这里下载代码.

公众号

相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
{:.info}

Email yxliphy@gmail.com