使用动态数组造成数组越界的问题
当我使用动态数组的时候，也就是下面的tH，首先我在ELEMENTS中声明了这个数组，然后我定义了赋值规则，tH的大小应该是(N,N,2,2)，后面我在子程序GETH中完成对tH的具体赋值，这里产生了数组越界的错误。我的问题首先在于，这里对动态数组的大小命令也就是“allocate(tH(N, N, 2,2))”应该在哪个子程序中呢，或者是不需要声明数组的大小，如何判断是否需要声明大小呢。
其次就在于这个数组越界的错误如何产生的
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*, )
complex(8), parameter :: Tqtr(2 , 2) = reshape(w*[(1.d0, 0.d0), exp(CI*Phi) , exp(-CI*Phi), (1.d0, 0.d0)], )
complex(8), parameter :: Tqtl(2 , 2) = reshape(w*[(1.d0, 0.d0), exp(-CI*Phi) , exp(CI*Phi), (1.d0, 0.d0)], )
! geometry of graphene
real(8), parameter :: a1(2) = a*, a2(2) = a*[-1.d0/2.d0, sqrt(3.d0)/2.d0]
real(8), parameter :: b1(2) = ((4.d0*PI)/(sqrt(3.d0)*a))*,&
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)], )
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)], )
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, n2ket index
integer :: i, j ! Hamiltonian matrix index
b1m = (8.d0*PI*sin(theta/2.d0)/(sqrt(3.d0)*a))*
b2m = (8.d0*PI*sin(theta/2.d0)/(sqrt(3.d0)*a))*
qb = (8*PI*sin(theta/2.d0)/(3.d0*a))*
qtr = (8*PI*sin(theta/2.d0)/(3.d0*a))* ! 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
页:
[1]