[求助] 使用动态数组造成数组越界的问题

module TBG
    implicit none

    complex(8), parameter :: CI = (0.0, 1.0)
    real(8), parameter :: PI = acos(-1.0) , d = 1.42d-010, a = sqrt(3.d0)*d, eV = 1.602176634D-19, hbar = 1.054572663D-034
    real(8), parameter ::  Phi = 2.d0*PI/3.d0, vF = 5.944d-010*eV, w = 0.11D0*eV !(vF = vF/hbar)

    complex(8), parameter :: Tqb(2 , 2) = reshape(w*[1.d0 , 1.d0 , 1.d0 , 1.d0], [2 , 2])
    complex(8), parameter :: Tqtr(2 , 2) = reshape(w*[(1.d0, 0.d0), exp(CI*Phi) , exp(-CI*Phi), (1.d0, 0.d0)], [2, 2])
    complex(8), parameter :: Tqtl(2 , 2) = reshape(w*[(1.d0, 0.d0), exp(-CI*Phi) , exp(CI*Phi), (1.d0, 0.d0)], [2, 2])

    ! geometry of graphene
    real(8), parameter :: a1(2) = a*[1.d0/2.d0, sqrt(3.d0)/2.d0], a2(2) = a*[-1.d0/2.d0, sqrt(3.d0)/2.d0]
    real(8), parameter :: b1(2) = ((4.d0*PI)/(sqrt(3.d0)*a))*[sqrt(3.d0)/2.d0, 1.d0/2.d0],&
                          b2(2) = ((4.d0*PI)/(sqrt(3.d0)*a))*[-sqrt(3.d0)/2.d0, 1.d0/2.d0]

    ! truncation and order of Hamiltonian
    integer, parameter :: tr = 1 , N =  2*((2*tr+1)**2)

contains

! Singel Layer Hamiltonian of Graphene
subroutine HSLG(k, theta, qb, b1m, b2m, layer, m1, m2, M)
    implicit none

    integer, intent(in) :: m1, m2, layer
    real(8), intent(in) :: k(2), theta, qb(2), b1m(2), b2m(2)
    real(8) :: k1, k2
    complex(8), intent(out) :: M(2,2) 

    real(8) :: theta_k

    if (layer == 1) then ! bottom layer
        
        k1 = k(1) - qb(1) - m1*b1m(1) - m2*b2m(1)
        k2 = k(2) - qb(2) - m1*b1m(2) - m2*b2m(2)
        theta_k = atan2(k2, k1)

        M = vF*sqrt(k1**2 + k2**2)*reshape([(0.d0, 0.d0), exp(CI*(theta_k + theta)), &
        exp(-CI*(theta_k + theta)) , (0.d0, 0.d0)], [2, 2])
    
    elseif (layer == 2) then ! top layer

        k1 = k(1) - m1*b1m(1) - m2*b2m(1)
        k2 = k(2) - m1*b1m(2) - m2*b2m(2)
        theta_k = atan2(k2, k1)

        M = vF*sqrt(k1**2 + k2**2)*reshape([(0.d0, 0.d0), exp(CI*(theta_k - theta)),&
         exp(-CI*(theta_k - theta)) , (0.d0, 0.d0)], [2, 2])
    endif

end subroutine HSLG

subroutine INDEX(layer, n1, n2, i) ! layer = 1(bottom), 2(top)
    implicit none

    integer, intent(in) :: layer, n1, n2
    integer, intent(out) :: i

    i = (layer - 1)*((2*tr+1)**2)+(n1+tr)*(2*tr+1)+(n2+tr)+1

end subroutine INDEX

subroutine ELEMENTS(tH, alpha, n1, n2, beta, m1, m2, M)
    implicit none

    complex(8), allocatable, intent(inout) :: tH(: , :, :, :)
    integer, intent(in) :: alpha, n1, n2, beta, m1, m2 
    complex(8), intent(in) :: M(2, 2)
    integer :: i, j

    allocate(tH(N, N, 2,2))
  
    call INDEX(alpha, m1, m2, i) ! bra 
    call INDEX(beta, n1, n2, j) ! ket
    if (alpha == beta .and. m1 == n1 .and. m2 == n2) then
        tH(i, j, :, :) = M 
    else
        tH(i, j, :, :) = M 
        tH(j, i, :, :) = transpose(conjg(M))
    endif
endsubroutine ELEMENTS

subroutine GETH(k, theta, H)

    implicit none
    
    real(8), intent(in) :: k(2), theta
    complex(8), allocatable, intent(out) :: H(:, :) 
    complex(8), allocatable :: tH(:, :, :, :)
    complex(8) :: M(2, 2)

    real(8) :: b1m(2), b2m(2)
    real(8) :: qb(2), qtr(2), qtl(2)

    integer :: alpha, beta, n1, n2 ! n1, n2  ket index
    integer :: i, j    ! Hamiltonian matrix index

    b1m = (8.d0*PI*sin(theta/2.d0)/(sqrt(3.d0)*a))*[1.d0/2.d0 , -sqrt(3.d0)/2.d0]
    b2m = (8.d0*PI*sin(theta/2.d0)/(sqrt(3.d0)*a))*[1.d0/2.d0 , sqrt(3.d0)/2.d0]
    qb = (8*PI*sin(theta/2.d0)/(3.d0*a))*[0.d0 , -1.d0]
    qtr = (8*PI*sin(theta/2.d0)/(3.d0*a))*[sqrt(3.d0)/2.d0 , 1.d0/2.d0] ! q_tr
    qtl = (8*PI*sin(theta/2.d0)/(3.d0*a))*[-sqrt(3.d0)/2.d0 , 1.d0/2.d0] !q_tl

    tH = 0.d0

    do alpha = 1, 2
        do beta = 1, 2
            do n1 = -tr, tr
                do n2 = -tr, tr
                    call INDEX(alpha, n1, n2, i)

                if  (alpha ==1 .and. beta == 1) then       
                    call ELEMENTS(tH, alpha, n1, n2, beta, n1, n2, M)
                elseif (alpha ==2 .and. beta == 2) then
                    call HSLG(k, theta/2, qb, b1m, b2m, 2, n1, n2, M)
                    call ELEMENTS(tH, alpha, n1, n2, beta, n1, n2, M) 
                elseif (alpha == 1 .and. beta == 2) then ! Tqb: n1 = m1, n2 = m2
                        call ELEMENTS(tH, alpha, n1, n2, beta, n1, n2, Tqb)
                    if (abs(n2 - 1) <= tr) then ! T2: m1 = n1, m2 = n2 - 1
                        call ELEMENTS(tH, alpha, n1, n2, beta, n1, n2 - 1, Tqtr)
                    elseif (abs(n1 + 1) <= tr) then  ! T3: m1 = n1 + 1, m2 = n2
                        call ELEMENTS(tH, alpha, n1, n2, beta, n1 + 1, n2, Tqtl)
                    end if
                end if
                end do
            end do
        end do
    end do

    H = 0.d0
    do i = 1, N
        do j = 1, N
            H(2*i-1:2*i, 2*j-1:2*j) = tH(i, j, : ,:)
        enddo
    enddo

    deallocate(tH)

end subroutine GETH

end module TBG


program BM_MODEL
    use TBG

    implicit none
    ! parameter
    integer, parameter :: numk = 100

    integer :: i, j
    real(8) :: l, k(2)
    real(8) :: kx(3*numk), ky(3*numk), dk(3*numk)

    ! matrix
    complex(8), allocatable :: H(:, :)
    real(8) :: theta, En(2*N)

    theta = (1.05d0)*PI/180.d0

    l = 8*PI*sin(theta/2.d0)/(3.d0*a)
    do i = 1, 3*numk
        if (i == 1) then
            kx(i) = 0.d0
            ky(i) = -1.d0*l
            dk(i) = 0.d0
        else
        if ( i <= numk ) then
            kx(i) = (sqrt(3.d0)*i/(2*numk))*l
            ky(i) = (real(i)/(2.d0*numk) - 1)*l
        elseif ( i<= 2*numk) then
            kx(i) = sqrt(3.d0)*l*(1.d0 - real(i)/(2*numk))
            ky(i) = -l/2.d0
        elseif (i <= 3*numk) then
            kx(i) = 0.d0
            ky(i) = (l/2.d0)*(real(i)/numk - 3)       
        endif
        dk(i) = dk(i - 1) + sqrt((kx(i) - kx(i-1))**2 +(ky(i) - ky(i-1))**2)
    endif
    enddo

    open(1, file = 'eb.dat' , status = 'replace')

    do i = 1, 3*numk
            k(1) = kx(i)
            k(2) = ky(i)
            call GETH(k, theta, H)
            call myzheev(2*N, En, H)
            write(1,'(200es32.16,1x)')dk(i), (En(j)/eV, j = 1, 2*N)
    enddo

    deallocate(H)

    close(1)
   
endprogram BM_MODEL

subroutine myzheev(N,W,H)
    implicit none
    integer,intent(in) :: N
    real(8),intent(out) :: W(N)
    complex(8),intent(inout) :: H(N,N)
    integer :: INFO=0
    real(8) :: RWORK(3*N-2)
    complex(8) :: WORK(3*N)

    call ZHEEV('V', 'U', N, H, N, W, WORK, 3*N, RWORK, INFO)

endsubroutine myzheev