前面有两篇博客主要关注的不同的基矢下如何构建BdG形式的哈密顿量,既然只不过是基矢的变化,那么两个哈密顿量肯定就是等价的,这里就给一个具体的实例,来看看怎么把两个不同基矢下的哈密顿量进行相互的改写,并利用程序验证其正确性.
{:.info}

模型

这里选用参考(5)文献中的哈密顿量

$\begin{equation} \begin{aligned} H(\mathbf{k})&=2\lambda_x\sin k_x\sigma_xs_z\tau_z+2\lambda_y\sin k_y\sigma_y\tau_z+(\xi_{\bf k}\sigma_z-\mu)\tau_z+\Delta_0\tau_x+{\bf h\cdot s}\\ \xi_{\bf k}&=\epsilon_0-2t_x\cos k_x-2t_y\cos k_y \end{aligned} \end{equation}\label{h1}$

这里的基矢选择为$\Psi^\dagger=(c^\dagger_{a\uparrow\mathbf{k}},c^\dagger_{b\uparrow\mathbf{k}},c^\dagger_{a\downarrow\mathbf{k}},c^\dagger_{b\downarrow\mathbf{k}},c_{a\downarrow\mathbf{-k}},c_{b\downarrow\mathbf{-k}},-c_{a\uparrow\mathbf{-k}},-c_{b\uparrow\mathbf{-k}})=(C_\mathbf{k}^\dagger,-is_y\sigma_0C_\mathbf{-k})$

在超导态基矢选择对构建BdG哈密顿量的影响这篇博客中我也提到过,这种基矢的选取虽然从形式上分析是比较方便的,但是在写程序的时候,不仅仅需要考虑哈密顿量中的符号,基矢中的符号也同时需要考虑,根据我自己的习惯,这时候写程序就会变的有些麻烦,还是想选用$\Psi^\dagger=(c^\dagger_{a\uparrow\mathbf{k}},c^\dagger_{b\uparrow\mathbf{k}},c^\dagger_{a\downarrow\mathbf{k}},c^\dagger_{b\downarrow\mathbf{k}},c_{a\downarrow\mathbf{-k}},c_{b\downarrow\mathbf{-k}},c_{a\uparrow\mathbf{-k}},c_{b\uparrow\mathbf{-k}})$这种基矢,这样的话就只需要考虑哈密顿量的具体形式即可,不需要额外再去考虑基矢中的符号.

同样利用超导态基矢选择对构建BdG哈密顿量的影响这篇文章中的方法,来把BdG形式的构建过程算一下,具体我就不演示这个哈密顿量是怎么做的了,可以参考前面的文章,做法都是完全相同的,因为正常态的基矢选择都是相同的,唯一需要考虑的就是超导的时候,空穴部分的基矢是如何选取的,完成之后,哈密顿量可以改写为

$\begin{equation} \begin{aligned} H^{\mathrm{BdG}}(\mathbf{k})&=(\xi_{\bf k}\sigma_z-\mu +{\bf h\cdot s})\tau_z+2\lambda_x\sin k_x\sigma_xs_z+2\lambda_y\sin k_y\sigma_y\tau_z+\Delta_0s_y\tau_y\\ \xi_{\bf k}&=\epsilon_0-2t_x\cos k_x-2t_y\cos k_y \end{aligned} \end{equation}\label{h2}$

简记$\Gamma_1=\sigma_zs_0\tau_z,\Gamma_2=\sigma_0s_x\tau_z,\Gamma_3=\sigma_xs_z\tau_0,\Gamma_4=\sigma_ys_0\tau_z,\Gamma_5=\sigma_0s_y\tau_y$

写程序计算

Cylinder Geomotery

!   https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.124.227001
!   将基矢写作通常的情况
    module pub
    implicit none
    integer xn,N,hn,kn
    parameter(xn = 50,kn = 50,hn = 8,N = xn*hn)
    real,parameter::pi = 3.1415926535
    complex,parameter::im = (0.0,1.0)  !虚数单位
    complex::Ham(N,N) = 0
!=================================
    real m0   !Driac mass
    real tx,ty
    real lamx,lamy
    real del0
    real h  ! Zeeman field
    complex g1(hn,hn),g2(hn,hn),g3(hn,hn),g4(hn,hn),g5(hn,hn),g6(hn,hn)
!================MKL===============
    integer::lda = N
    integer,parameter::lwmax=2*N+N**2
    real,allocatable::w(:)
    complex,allocatable::work(:)
    real,allocatable::rwork(:)
    integer,allocatable::iwork(:)
    integer lwork   ! at least 2*N+N**2
    integer lrwork    ! at least 1 + 5*N +2*N**2
    integer liwork   ! at least 3 +5*N
    integer info
    end module pub
!================= PROGRAM START ============================
    program ex01
    use pub
    integer m,l,i
    real ky
    !=====================
    allocate(w(N))
    allocate(work(lwmax))
    allocate(rwork(1+5*N+2*N**2))
    allocate(iwork(3+5*N))
    !--------------------------
    m0 = 1.0
    tx = 1.0
    ty = 1.0
    lamx = 1.0
    lamy = 1.0
    del0 = 0.4
    h = 0.4
    call band()
    stop 
    end program
!========================================================================
    subroutine band()
    use pub
    real ky
    integer ik
    open(20,file="x-open.dat")
    open(21,file="y-open.dat")
    do ik = -kn,kn
        ky = ik*pi/kn
        call openx(ky)
        write(20,999)ky/pi,(w(i),i = 1,N)
        call openy(ky)
        write(21,999)ky/pi,(w(i),i = 1,N)
    end do
    close(20)
    close(21)
999 format(4001f11.6)
    end subroutine band
!======================== Pauli Matrix driect product============================
    subroutine Pauli()
    use pub
    !----mass term
    g1(1,1) = 1
    g1(2,2) = -1
    g1(3,3) = 1
    g1(4,4) = -1
    g1(5,5) = -1
    g1(6,6) = 1
    g1(7,7) = -1
    g1(8,8) = 1
    !----Zeeman field
    g2(1,3) = 1
    g2(2,4) = 1
    g2(3,1) = 1
    g2(4,2) = 1
    g2(5,7) = -1
    g2(6,8) = -1
    g2(7,5) = -1
    g2(8,6) = -1
    !--- sin(kx)
    g3(1,2) = 1
    g3(2,1) = 1
    g3(3,4) = -1
    g3(4,3) = -1
    g3(5,6) = 1
    g3(6,5) = 1
    g3(7,8) = -1
    g3(8,7) = -1
    !----sin(ky)
    g4(1,2) = -im
    g4(2,1) = im
    g4(3,4) = -im
    g4(4,3) = im
    g4(5,6) = im
    g4(6,5) = -im
    g4(7,8) = im
    g4(8,7) = -im
    !----Delta(k)
    g5(1,7) = -1
    g5(2,8) = -1
    g5(3,5) = 1
    g5(4,6) = 1
    g5(5,3) = 1
    g5(6,4) = 1
    g5(7,1) = -1
    g5(8,2) = -1
    !---mu
    g6(1,1) = 1
    g6(2,2) = 1
    g6(3,3) = 1
    g6(4,4) = 1
    g6(5,5) = -1
    g6(6,6) = -1
    g6(7,7) = -1
    g6(8,8) = -1
    return
    end subroutine Pauli
!==========================================================
    subroutine openx(ky)
    use pub
    real ky
    integer m,l,k
    call Pauli()
    Ham = 0
    !=================== Non-diag Term========================
    do k = 0,xn-1
        if (k == 0) then
            do m = 1,hn
                do l = 1,hn
                    Ham(m,l) = 2*lamy*sin(ky)*g4(m,l) + (m0 - 2*ty*cos(ky))*g1(m,l) + del0*g5(m,l) + h*g2(m,l) + mu*g6(m,l)
                    Ham(m,l + 8) = -im*lamx*g3(m,l) - tx*g1(m,l)
                end do
            end do
        elseif ( k == xn-1 ) then
            do m = 1,hn
                do l = 1,hn
                    Ham(k*hn + m,k*hn + l) = 2*lamy*sin(ky)*g4(m,l) + (m0 - 2*ty*cos(ky))*g1(m,l) + del0*g5(m,l) + h*g2(m,l) + mu*g6(m,l)
                    Ham(k*hn + m,k*hn + l - hn) = im*lamx*g3(m,l) - tx*g1(m,l)
                end do
            end do
        else
            do m = 1,hn
                do l = 1,hn
                    Ham(k*hn + m,k*hn + l) = 2*lamy*sin(ky)*g4(m,l) + (m0 - 2*ty*cos(ky))*g1(m,l) + del0*g5(m,l) + h*g2(m,l) + mu*g6(m,l)
                    Ham(k*hn + m,k*hn + l + hn) = -im*lamx*g3(m,l) - tx*g1(m,l)
                    Ham(k*hn + m,k*hn + l - hn) = im*lamx*g3(m,l) - tx*g1(m,l)
                end do
            end do
        end if
    end do
    call ishermitian()
    call eigsol()
    end subroutine openx
!==========================================================
    subroutine openy(kx)
    use pub
    real kx
    integer m,l,k
    call Pauli()
    Ham = 0
    !=================== Non-diag Term========================
    do k = 0,xn-1
        if (k == 0) then
            do m = 1,hn
                do l = 1,hn
                    Ham(m,l) = 2*lamx*sin(kx)*g3(m,l) + (m0 - 2*tx*cos(kx))*g1(m,l) + del0*g5(m,l) + h*g2(m,l) + mu*g6(m,l)
                    Ham(m,l + 8) = -im*lamy*g4(m,l) - ty*g1(m,l)
                end do
            end do
        elseif ( k == xn-1 ) then
            do m = 1,hn
                do l = 1,hn
                    Ham(k*hn + m,k*hn + l) = 2*lamx*sin(kx)*g3(m,l) + (m0 - 2*tx*cos(kx))*g1(m,l) + del0*g5(m,l) + h*g2(m,l) + mu*g6(m,l)
                    Ham(k*hn + m,k*hn + l - hn) = im*lamy*g4(m,l) - ty*g1(m,l)
                end do
            end do
        else
            do m = 1,hn
                do l = 1,hn
                    Ham(k*hn + m,k*hn + l) = 2*lamx*sin(kx)*g3(m,l) + (m0 - 2*tx*cos(kx))*g1(m,l) + del0*g5(m,l) + h*g2(m,l) + mu*g6(m,l)
                    Ham(k*hn + m,k*hn + l + hn) = -im*lamy*g4(m,l) - tx*g1(m,l)
                    Ham(k*hn + m,k*hn + l - hn) = im*lamy*g4(m,l) - tx*g1(m,l)
                end do
            end do
        end if
    end do
    call ishermitian()
    call eigsol()
    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')
                ccc = ccc +1
                write(160,*)i,j
                write(160,*)Ham(i,j)
                write(160,*)Ham(j,i)
                write(160,*)"===================="
                stop
            end if
        end do
    end do
    close(160)
    return
    end subroutine ishermitian
!================= 矩阵本征值求解 ==============
    subroutine eigSol()
    use pub
    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.txt",status="unknown")
        write(110,*)'The algorithm failed to compute eigenvalues.'
        close(110)
    end if
    open(120,file="eigval.dat",status="unknown")
    do m = 1,N
        write(120,*)m,w(m)
    end do
    close(120)
    return
    end subroutine eigSol

这个程序是用来计算边界态的,即一个方向是开边界,另外一个方向是周期边界条件,我已验证过,相同的参数下和文章中结果是相同.

Real-space

!   https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.124.227001
!   将基矢写作通常的情况
    module pub
    implicit none
    integer xn,N,yn,hn,kn,len2
    parameter(xn = 36,yn = xn,kn = 50,hn = 8,len2 = xn*yn,N = len2*hn)
    real eps,pi
    parameter (pi = 3.1415926535, eps = 1e-5)
    complex,parameter::im = (0.0,1.0)  !虚数单位
    complex::Ham(N,N) = 0
    !-----------------------------------------
    real m0,tx,ty,lamx,lamy,del0,h
    integer bry(4,len2)  ! boundary index
    !------------------------------------------
    integer::lda = N
    integer,parameter::lwmax=2*N+N**2
    real,allocatable::w(:)
    complex,allocatable::work(:)
    real,allocatable::rwork(:)
    integer,allocatable::iwork(:)
    integer lwork   ! at least 2*N+N**2
    integer lrwork    ! at least 1 + 5*N +2*N**2
    integer liwork   ! at least 3 +5*N
    integer info
    end module pub
!================= PROGRAM START ============================
    program ex01
    use pub
    integer m,l,i
    !=====================
    allocate(w(N))
    allocate(work(lwmax))
    allocate(rwork(1+5*N+2*N**2))
    allocate(iwork(3+5*N))
    !--------------------------
    m0 = 1.0
    tx = 1.0
    ty = 1.0
    lamx = 1.0
    lamy = 1.0
    del0 = 0.4
    h = 0.8
    call matset2()
    call ldos()
    stop 
    end program
!===========================Local Density of State=============================
    subroutine ldos()
    use pub
    integer ix,iy
    real,external::wave
    open(12,file="wavenorm.dat") 
    do iy = 1,yn
        do ix = 1,xn
            write(12,*)ix,iy,wave(ix,iy)
        end do
    end do
    close(12)
    return
    end subroutine ldos
!=======================================================================
    real function wave(x,y)
    use pub
    integer x,y,ind
    ind = (y-1)*xn + x
    wave = 0
    do i2 = 0,7
        do zo = N/2 - 3,N/2 + 4
            wave = wave + abs(Ham(ind + len2*i2,zo))**2
        end do
    end do
    end function wave
!==========================================================
    subroutine boundary()
    use pub
    integer ix,iy,i
    bry = 0
    do iy = 1,yn
        do ix = 1,xn
            i = (iy - 1)*xn + ix
            bry(1,i) = i + 1
            if(ix == xn)bry(1,i) = bry(1,i) - xn
            bry(2,i) = i - 1
            if(ix == 1)bry(2,i) = bry(2,i) + xn
            bry(3,i) = i + yn
            if(iy == yn)bry(3,i) = bry(3,i) - len2
            bry(4,i) = i - yn
            if(iy == 1)bry(4,i) = bry(4,i) + len2
        end do
    end do
    return
    end subroutine boundary
!=====================================================================
    subroutine matset2()
    use pub
    integer ix,iy,i
    call boundary()
    Ham = 0
    do iy = 1,yn
        do ix = 1,xn
            i = (iy - 1)*xn + ix
            !-------------------------------
            !  Mass term
            !(1,1)
            ham(i,i) = m0 + mu
            if(ix.ne.xn)ham(i,bry(1,i)) = tx
            if(ix.ne.1)ham(i,bry(2,i)) = tx
            if(iy.ne.yn)ham(i,bry(3,i)) = ty
            if(iy.ne.1)ham(i,bry(4,i)) = ty
            !(2,2)
            ham(len2 + i,len2 + i) = -m0 + mu
            if(ix.ne.xn)ham(len2 + i,len2 + bry(1,i)) = -tx
            if(ix.ne.1)ham(len2 + i,len2 + bry(2,i)) = -tx
            if(iy.ne.yn)ham(len2 + i,len2 + bry(3,i)) = -ty
            if(iy.ne.1)ham(len2 + i,len2 + bry(4,i)) = -ty
            !(3,3)
            ham(len2*2 + i,len2*2 + i) = m0 + mu
            if(ix.ne.xn)ham(len2*2 + i,len2*2 + bry(1,i)) = tx
            if(ix.ne.1)ham(len2*2 + i,len2*2 + bry(2,i)) = tx
            if(iy.ne.yn)ham(len2*2 + i,len2*2 + bry(3,i)) = ty
            if(iy.ne.1)ham(len2*2 + i,len2*2 + bry(4,i)) = ty
            !(4,4)
            ham(len2*3 + i,len2*3 + i) = -m0 + mu
            if(ix.ne.xn)ham(len2*3 + i,len2*3 + bry(1,i)) = -tx
            if(ix.ne.1)ham(len2*3 + i,len2*3 + bry(2,i)) = -tx
            if(iy.ne.yn)ham(len2*3 + i,len2*3 + bry(3,i)) = -ty
            if(iy.ne.1)ham(len2*3 + i,len2*3 + bry(4,i)) = -ty
            !(5,5)
            ham(len2*4 + i,len2*4 + i) = -m0 - mu
            if(ix.ne.xn)ham(len2*4 + i,len2*4 + bry(1,i)) = -tx
            if(ix.ne.1)ham(len2*4 + i,len2*4 + bry(2,i)) = -tx
            if(iy.ne.yn)ham(len2*4 + i,len2*4 + bry(3,i)) = -ty
            if(iy.ne.1)ham(len2*4 + i,len2*4 + bry(4,i)) = -ty
            !(6,6)
            ham(len2*5 + i,len2*5 + i) = m0 - mu
            if(ix.ne.xn)ham(len2*5 + i,len2*5 + bry(1,i)) = tx
            if(ix.ne.1)ham(len2*5 + i,len2*5 + bry(2,i)) = tx
            if(iy.ne.yn)ham(len2*5 + i,len2*5 + bry(3,i)) = ty
            if(iy.ne.1)ham(len2*5 + i,len2*5 + bry(4,i)) = ty
            !(7,7)
            ham(len2*6 + i,len2*6 + i) = -m0 - mu
            if(ix.ne.xn)ham(len2*6 + i,len2*6 + bry(1,i)) = -tx
            if(ix.ne.1)ham(len2*6 + i,len2*6 + bry(2,i)) = -tx
            if(iy.ne.yn)ham(len2*6 + i,len2*6 + bry(3,i)) = -ty
            if(iy.ne.1)ham(len2*6 + i,len2*6 + bry(4,i)) = -ty
            !(8,8)
            ham(len2*7 + i,len2*7 + i) = m0 - mu
            if(ix.ne.xn)ham(len2*7 + i,len2*7 + bry(1,i)) = tx
            if(ix.ne.1)ham(len2*7 + i,len2*7 + bry(2,i)) = tx
            if(iy.ne.yn)ham(len2*7 + i,len2*7 + bry(3,i)) = ty
            if(iy.ne.1)ham(len2*7 + i,len2*7 + bry(4,i)) = ty
            !----------------------------------------------
            !   SOC
            !(1,2)
            if(ix.ne.xn)ham(i,len2 + bry(1,i)) = lamx/im
            if(ix.ne.1)ham(i,len2 + bry(2,i)) = -lamx/im
            if(iy.ne.yn)ham(i,len2 + bry(3,i)) = -im*lamx/im
            if(iy.ne.1)ham(i,len2 + bry(4,i)) = im*lamx/im
            !(2,1)
            if(ix.ne.xn)ham(len2 + i,bry(1,i)) = lamx/im
            if(ix.ne.1)ham(len2 + i,bry(2,i)) = -lamx/im
            if(iy.ne.yn)ham(len2 + i,bry(3,i)) = im*lamx/im
            if(iy.ne.1)ham(len2 + i,bry(4,i)) = -im*lamx/im
            !(3,4)
            if(ix.ne.xn)ham(len2*2 + i,len2*3 + bry(1,i)) = -lamx/im
            if(ix.ne.1)ham(len2*2 + i,len2*3 + bry(2,i)) = lamx/im
            if(iy.ne.yn)ham(len2*2 + i,len2*3 + bry(3,i)) = -im*lamx/im
            if(iy.ne.1)ham(len2*2 + i,len2*3 + bry(4,i)) = im*lamx/im
            !(4,3)
            if(ix.ne.xn)ham(len2*3 + i,len2*2 + bry(1,i)) = -lamx/im
            if(ix.ne.1)ham(len2*3 + i,len2*2 + bry(2,i)) = lamx/im
            if(iy.ne.yn)ham(len2*3 + i,len2*2 + bry(3,i)) = im*lamx/im
            if(iy.ne.1)ham(len2*3 + i,len2*2 + bry(4,i)) = -im*lamx/im
            !(5,6)
            if(ix.ne.xn)ham(len2*4 + i,len2*5 + bry(1,i)) = lamx/im
            if(ix.ne.1)ham(len2*4 + i,len2*5 + bry(2,i)) = -lamx/im
            if(iy.ne.yn)ham(len2*4 + i,len2*5 + bry(3,i)) = im*lamx/im
            if(iy.ne.1)ham(len2*4 + i,len2*5 + bry(4,i)) = -im*lamx/im
            !(6,5)
            if(ix.ne.xn)ham(len2*5 + i,len2*4 + bry(1,i)) = lamx/im
            if(ix.ne.1)ham(len2*5 + i,len2*4 + bry(2,i)) = -lamx/im
            if(iy.ne.yn)ham(len2*5 + i,len2*4 + bry(3,i)) = -im*lamx/im
            if(iy.ne.1)ham(len2*5 + i,len2*4 + bry(4,i)) = im*lamx/im
            !(7,8)
            if(ix.ne.xn)ham(len2*6 + i,len2*7 + bry(1,i)) = -lamx/im
            if(ix.ne.1)ham(len2*6 + i,len2*7 + bry(2,i)) = lamx/im
            if(iy.ne.yn)ham(len2*6 + i,len2*7 + bry(3,i)) = im*lamx/im
            if(iy.ne.1)ham(len2*6 + i,len2*7 + bry(4,i)) = -im*lamx/im
            !(8,7)
            if(ix.ne.xn)ham(len2*7 + i,len2*6 + bry(1,i)) = -lamx/im
            if(ix.ne.1)ham(len2*7 + i,len2*6 + bry(2,i)) = lamx/im
            if(iy.ne.yn)ham(len2*7 + i,len2*6 + bry(3,i)) = -im*lamx/im
            if(iy.ne.1)ham(len2*7 + i,len2*6 + bry(4,i)) = im*lamx/im
            !---------------------------------------------------
            !  SC
            !(1,7)
            ham(i,len2*6 + i) = -del0
            !(2,8)
            ham(len2 + i,len2*7 + i) = -del0
            !(3,5)
            ham(len2*2 + i,len2*4 + i) = del0
            !(4,6)
            ham(len2*3 + i,len2*5 + i) = del0
            !(5,3)
            ham(len2*4 + i,len2*2 + i) = del0
            !(6,4)
            ham(len2*5 + i,len2*3 + i) = del0
            !(7,1)
            ham(len2*6 + i,i) = -del0
            !(8,2)
            ham(len2*7 + i,len2 + i) = -del0
            !---------------------------------------------------
            ! Zeeman field
            !(1,3)
            ham(i,len2*2 + i) = h
            !(2,4)
            ham(len2 + i,len2*3 + i) = h
            !(3,1)
            ham(len2*2 + i,i) = h
            !(4,2)
            ham(len2*3 + i,len2 + i) = h
            !(5,7)
            ham(len2*4 + i,len2*6 + i) = -h
            !(6,8)
            ham(len2*5 + i,len2*7 + i) = -h
            !(7,5)
            ham(len2*6 + i,len2*4 + i) = -h
            !(8,6)
            ham(len2*7 + i,len2*5 + i) = -h
        end do
    end do
    !-------------------------
    call ishermitian()
    call eigsol()
    end subroutine matset2
!============================================================
    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')
                ccc = ccc +1
                write(160,*)i,j
                write(160,*)Ham(i,j)
                write(160,*)Ham(j,i)
                write(160,*)"===================="
                write(*,*)"Ham isn't Hermitian"
                stop
            end if
        end do
    end do
    close(160)
    return
    end subroutine ishermitian
!================= 矩阵本征值求解 ==============
    subroutine eigSol()
    use pub
    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",status="unknown")
    do m = 1,N
        write(120,*)m,w(m)
    end do
    close(120)
    return
    end subroutine eigSol