|
|
本帖最后由 Kieran 于 2025-1-14 00:47 编辑
大家好，
这是我写的一段Fortran代码，其中使用了ZHEEV函数对角化矩阵，以计算本征值。
[Fortran] 纯文本查看 复制代码 SUBROUTINE DOTPRODUCT(v1,v2,v3)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
REAL (KIND=dp) :: v1(3)
REAL (KIND=dp) :: v2(3)
REAL (KIND=dp) :: v3(3)
v3(1) = v1(1)*v2(1)
v3(2) = v1(2)*v2(2)
v3(3) = v1(3)*v2(3)
RETURN
END SUBROUTINE DOTPRODUCT
SUBROUTINE CROSSPRODUCT(v1,v2,v3)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
REAL (KIND=dp) :: v1(3)
REAL (KIND=dp) :: v2(3)
REAL (KIND=dp) :: v3(3)
v3(1) = v1(2)*v2(3)-v1(3)*v2(2)
v3(2) = v1(3)*v2(1)-v1(1)*v2(3)
v3(3) = v1(1)*v2(2)-v1(2)*v2(1)
RETURN
END SUBROUTINE CROSSPRODUCT
SUBROUTINE LATTICEINK(lr,lk)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
REAL (KIND=dp) :: lr(3,2) !matrix from the main code
REAL (KIND=dp) :: lk(3,2) !matrix back to the main code
REAL (KIND=dp) :: em(3,3) !meta results
REAL (KIND=dp) :: cr(3,3) !meta results
REAL (KIND=dp) :: vl(3) !Volume of the unit cell in the real space
REAL (KIND=dp) :: pi !constant pi
em = 0.0d0
em(:,1:2) = lr
em(3,3) = 1.5d1
pi = DACOS(-1.0d0)
CALL CROSSPRODUCT(em(:,2),em(:,3),cr(:,1))
CALL CROSSPRODUCT(em(:,3),em(:,1),cr(:,2))
CALL CROSSPRODUCT(em(:,1),em(:,2),cr(:,3))
CALL DOTPRODUCT(em(:,1),cr(:,1),vl)
vl(1) = DSQRT(vl(1)**2+vl(2)**2+vl(3)**2)
CALL DOTPRODUCT(em(:,2),cr(:,2),vl)
vl(1) = DSQRT(vl(1)**2+vl(2)**2+vl(3)**2)
CALL DOTPRODUCT(em(:,3),cr(:,3),vl)
vl(1) = DSQRT(vl(1)**2+vl(2)**2+vl(3)**2)
lk(:,1) = 2.0d0*pi*cr(:,1)/vl(1)
lk(:,2) = 2.0d0*pi*cr(:,2)/vl(1)
RETURN
END SUBROUTINE LATTICEINK
PROGRAM GRAPHENE_TIGHTBINDING
USE LAPACK
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,14)
INTEGER :: i, j
REAL (KIND=dp) :: dc !Distant between carbon atom in pristine graphene unit cell
REAL (KIND=dp) :: lr(3,2) !Lattice constant in real space
REAL (KIND=dp) :: lk(3,2) !Lattice constant in k space
REAL (KIND=dp) :: cr(3,3) !matrix to store the meta data
REAL (KIND=dp) :: pi
REAL (KIND=dp) :: kh(3,3) !High symmetry k point
INTEGER :: km(3) !k point sampling
REAL (KIND=dp), ALLOCATABLE :: kc(:,:) !k point coordinate
REAL (KIND=dp) :: mk(3,2)
REAL (KIND=dp) :: ep(2) !onsite
REAL (KIND=dp) :: hp !hopping
COMPLEX (KIND=dp) :: ha(2,2) !Hamiltonian
REAL (KIND=dp), ALLOCATABLE :: bd(:,:) !band
INTEGER :: nu, lda, lwork, info
REAL (KIND=dp), ALLOCATABLE :: ei(:), rwork(:)
COMPLEX (KIND=dp), ALLOCATABLE :: work(:)
!Parameters
dc = 1.42d0
pi = DACOS(-1.0d0)
lr = 0.0d0
lr(1,1) = DCOS(pi/6.0d0)*dc*2.0d0
lr(1,2) = DCOS(pi*5.0d0/6.0d0)*dc
lr(2,2) = DSIN(pi/6.0d0)*dc+dc
kh = 0.0d0
kh(1,2) = 1.0d0/3.0d0
kh(2,2) = 1.0d0/3.0d0
kh(1,3) = 5.0d-1
km(1) = 100
km(2) = 100
km(3) = 100
ep = 0.0d0
hp = 2.8d0
nu = 2
lda = nu
lwork = -1
!
OPEN (UNIT=3, FILE='bandstructure.dat', STATUS='UNKNOWN')
!Computing the lattice constant in k space
CALL LATTICEINK(lr,lk)
!
!Computing the k point coordinate
ALLOCATE (kc(4,km(1)+km(2)+km(3)))
kc = 0.0d0
mk(1,1) = kh(1,1)*lk(1,1)+kh(2,1)*lk(1,2)
mk(2,1) = kh(1,1)*lk(2,1)+kh(2,1)*lk(2,2)
mk(3,1) = kh(1,1)*lk(3,1)+kh(2,1)*lk(3,2)
mk(1,2) = kh(1,2)*lk(1,1)+kh(2,2)*lk(1,2)
mk(2,2) = kh(1,2)*lk(2,1)+kh(2,2)*lk(2,2)
mk(3,2) = kh(1,2)*lk(3,1)+kh(2,2)*lk(3,2)
DO i = 1, km(1), 1
kc(1,i) = mk(1,1)+DBLE(i-1)*(mk(1,2)-mk(1,1))/DBLE(km(1))
kc(2,i) = mk(2,1)+DBLE(i-1)*(mk(2,2)-mk(2,1))/DBLE(km(1))
kc(3,i) = mk(3,1)+DBLE(i-1)*(mk(3,2)-mk(3,1))/DBLE(km(1))
END DO
mk(1,1) = kh(1,2)*lk(1,1)+kh(2,2)*lk(1,2)
mk(2,1) = kh(1,2)*lk(2,1)+kh(2,2)*lk(2,2)
mk(3,1) = kh(1,2)*lk(3,1)+kh(2,2)*lk(3,2)
mk(1,2) = kh(1,3)*lk(1,1)+kh(2,3)*lk(1,2)
mk(2,2) = kh(1,3)*lk(2,1)+kh(2,3)*lk(2,2)
mk(3,2) = kh(1,3)*lk(3,1)+kh(2,3)*lk(3,2)
DO i = km(1)+1, km(1)+km(2), 1
kc(1,i) = mk(1,1)+DBLE(i-km(1)-1)*(mk(1,2)-mk(1,1))/DBLE(km(2))
kc(2,i) = mk(2,1)+DBLE(i-km(1)-1)*(mk(2,2)-mk(2,1))/DBLE(km(2))
kc(3,i) = mk(3,1)+DBLE(i-km(1)-1)*(mk(3,2)-mk(3,1))/DBLE(km(2))
END DO
mk(1,1) = kh(1,3)*lk(1,1)+kh(2,3)*lk(1,2)
mk(2,1) = kh(1,3)*lk(2,1)+kh(2,3)*lk(2,2)
mk(3,1) = kh(1,3)*lk(3,1)+kh(2,3)*lk(3,2)
mk(1,2) = kh(1,1)*lk(1,1)+kh(2,1)*lk(1,2)
mk(2,2) = kh(1,1)*lk(2,1)+kh(2,1)*lk(2,2)
mk(3,2) = kh(1,1)*lk(3,1)+kh(2,1)*lk(3,2)
DO i = km(1)+km(2)+1, km(1)+km(2)+km(3), 1
kc(1,i) = mk(1,1)+DBLE(i-km(1)-km(2)-1)*(mk(1,2)-mk(1,1))/DBLE(km(3))
kc(2,i) = mk(2,1)+DBLE(i-km(1)-km(2)-1)*(mk(2,2)-mk(2,1))/DBLE(km(3))
kc(3,i) = mk(3,1)+DBLE(i-km(1)-km(2)-1)*(mk(3,2)-mk(3,1))/DBLE(km(3))
END DO
!
!Building Hamiltonian and computing the band structure
ALLOCATE (bd(nu,km(1)+km(2)+km(3)))
ALLOCATE (ei(nu))
ALLOCATE (work(lwork))
ALLOCATE (rwork(3*nu-2))
ha = DCMPLX(0.0d0, 0.0d0)
ha(1,1) = DCMPLX(ep(1), 0.0d0)
ha(1,2) = hp*(1.0d0+EXP(DCMPLX(0.0d0,-1.0d0)*(kc(1,km(1)+1)*lr(1,1)&
+kc(2,km(1)+1)*lr(2,1)+kc(3,km(1)+1)*lr(3,1)))+EXP(DCMPLX(0.0d0,-1.0d0)&
*(kc(1,km(1)+1)*lr(1,2)+kc(2,km(1)+1)*lr(2,2)+kc(3,km(1)+1)*lr(3,2))))
ha(2,2) = DCMPLX(ep(2), 0.0d0)
ha(2,1) = DCONJG(ha(1,2))
print*,'outsideloop'
print*,'ha(1,1)',ha(1,1)
print*,'ha(1,2)',ha(1,2)
print*,'ha(2,1)',ha(2,1)
print*,'ha(2,2)',ha(2,2)
print*,'nu',nu
print*,'lda',lda
print*,'ei',ei
print*,'work',work
print*,'lwork',lwork
print*,'rwork',rwork
print*,'info',info
CALL ZHEEV('V','U',nu,ha,lda,ei,work,lwork,rwork,info)
print*,'ei,kc',ei,kc(:,km(1)+1)
lwork = INT(REAL(work(1)))
DEALLOCATE (work)
ALLOCATE (work(lwork))
DO i = 1, km(1)+km(2)+km(3), 1
rwork = 0.0d0
ei = 0.0d0
ha = DCMPLX(0.0d0, 0.0d0)
ha(1,1) = DCMPLX(ep(1), 0.0d0)
ha(1,2) = hp*(1.0d0+EXP(DCMPLX(0.0d0,-1.0d0)*(kc(1,i)*lr(1,1)&
+kc(2,i)*lr(2,1)+kc(3,i)*lr(3,1)))+EXP(DCMPLX(0.0d0,-1.0d0)&
*(kc(1,i)*lr(1,2)+kc(2,i)*lr(2,2)+kc(3,i)*lr(3,2))))
ha(2,2) = DCMPLX(ep(2), 0.0d0)
ha(2,1) = DCONJG(ha(1,2))
if (i == km(1)+1) then
print*,'within the loop at K'
print*,'ha(1,1)',ha(1,1)
print*,'ha(1,2)',ha(1,2)
print*,'ha(2,1)',ha(2,1)
print*,'ha(2,2)',ha(2,2)
print*,'nu',nu
print*,'lda',lda
print*,'ei',ei
print*,'work',work
print*,'lwork',lwork
print*,'rwork',rwork
print*,'info',info
end if
CALL ZHEEV('V','U',nu,ha,lda,ei,work,lwork,rwork,info)
if (i == km(1)+1) then
print*,'within the loop at K',ei,kc(:,i)
end if
bd(:,i) = ei
END DO
!
!Plotting band structure
kc(4,1) = 0.0d0
DO i = 2, km(1)+km(2)+km(3), 1
kc(4,i) = kc(4,i-1)+DSQRT((kc(1,i)-kc(1,i-1))**2+(kc(2,i)-kc(2,i-1))**2+(kc(3,i)-kc(3,i-1))**2)
END DO
DO i = 1, nu, 1
DO j = 1, km(1)+km(2)+km(3), 1
WRITE (UNIT=3, FMT='(F15.10,3X,F15.10)') kc(4,j), bd(i,j)
END DO
WRITE (UNIT=3, FMT=*)
END DO
!
DEALLOCATE (bd)
DEALLOCATE (ei)
DEALLOCATE (kc)
DEALLOCATE (rwork)
DEALLOCATE (work)
CLOSE (UNIT=3)
END PROGRAM GRAPHENE_TIGHTBINDING
这个程序中，我两次输出本征值数据，一次是在循环外，第172行，print*,'ei,kc',ei,kc(:,km(1)+1)；另一次是在循环内，第202行，print*,'within the loop at K',ei,kc(:,i)。这两次计算过程中ha复数矩阵都是相同的。但输出的本征值ei却不同。
输出结果如下：
outsideloop
ha(1,1) (0.0000000000000000,0.0000000000000000)
ha(1,2) (0.0000000000000000,-4.8497422611928558)
ha(2,1) (0.0000000000000000,4.8497422611928558)
ha(2,2) (0.0000000000000000,0.0000000000000000)
nu 2
lda 2
ei 0.0000000000000000 0.0000000000000000
work
lwork -1
rwork 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
info 0
ei,kc 0.0000000000000000 0.0000000000000000 0.85154899729305999 1.4749261284459125 0.0000000000000000 0.0000000000000000
within the loop at K
ha(1,1) (0.0000000000000000,0.0000000000000000)
ha(1,2) (0.0000000000000000,-4.8497422611928558)
ha(2,1) (0.0000000000000000,4.8497422611928558)
ha(2,2) (0.0000000000000000,0.0000000000000000)
nu 2
lda 2
ei 0.0000000000000000 0.0000000000000000
work (66.000000000000000,0.0000000000000000) (0.0000000000000000,0.0000000000000000) (32.000000000000000,0.0000000000000000) (0.0000000000000000,0.0000000000000000) (0.0000000000000000,0.0000000000000000) (0.0000000000000000,0.0000000000000000)
......
lwork 66
rwork 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
info 0
within the loop at K -4.8497422611928558 4.8497422611928558 0.85154899729305999 1.4749261284459125 0.0000000000000000 0.0000000000000000
可以看到循环外的本征都是0.0；而循环内的本征却是-4.8497422611928558和4.8497422611928558。两次的输入矩阵都是完全相同的。
请问是哪里不合理呢？望指教。
谢谢。
|
|
F 币
贡献
权杖
窥视卡
发表于 2025-1-14 00:33:44
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
楼主