Fortran Coder

查看: 4373|回复: 1
打印 上一主题 下一主题

[通用算法] 石墨烯边缘态能带

[复制链接]

2

帖子

1

主题

0

精华

新人

F 币
27 元
贡献
9 点
跳转到指定楼层
楼主
发表于 2021-10-19 13:52:33 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
计算石墨烯边缘态能带的时候,我用32个点可以拟合能带图形,但是在64个点能带的波动很大,是哪里出问题了呢?文章的doi:10.1103/PhysRevLett.102.096801
代码如下
[Fortran] 纯文本查看 复制代码
include 'cxml_include.f90'
	  parameter (ny=32,nx=32,nxy=nx)
        complex ham(nxy,nxy),eivec(nxy,nxy)
        real eival(nxy),pi
	
	pi=2.0*acos(0.0)
	t1=-1.0
      dmu=0.0
	dmu1=dmu+0.03
	dmu2=dmu-0.03   
      sq3=sqrt(3.0)

        open(14,file='32.dat',status='unknown')
	do 151 idy=0,ny 
	dky=float(idy)*2.0*pi/(ny*sq3) 
		  do 152 i=1,nxy
        do 152 j=1,nxy
         ham(i,j)=0.0
152     enddo
             	  
	do 153 ix=1,nx     
         ipx=ix+1     
	   imx=ix-1

	if(mod(ix,2).eq.0) then
           ham(ix,ix)=-dmu2    
	endif

	if(mod(ix,2).ne.0) then
           ham(ix,ix)=-dmu1      
	endif

CCCCCCCCCC n.n. hopping--------------------------------------------   
	if(mod(ix,2).eq.0) then

	      if(ipx.le.nx) then
	      ham(ix,ipx)=-t1
	      endif

	      if(imx.ge.1) then
	      ham(ix,imx)=-2.0*t1*cos(sq3*dky/2.0)
	      endif

	endif

c---------------------------------------------
	if(mod(ix,2).ne.0) then

	      if(imx.ge.1) then
	      ham(ix,imx)=-t1
	      endif

	      if(ipx.le.nx) then
	      ham(ix,ipx)=-2.0*t1*cos(sq3*dky/2.0)
	      endif

	endif
153	continue !scan ix

c-------diagonization
        neigen=nxy
        call eigen(neigen,ham,eival,eivec)

c      do ixe=1,nxy4
	write(14,990),dky*sq3,eival(5),eival(6),eival(7)
     &	,eival(8),eival(9),eival(10),eival(11),eival(12)
     &	,eival(13),eival(14),eival(15),eival(16),eival(17)
     &    ,eival(18),eival(19),eival(20),eival(21),eival(22)
     &    ,eival(23),eival(24),eival(25),eival(26),eival(27)
     &    ,eival(28)
c      enddo

990   format(1x,26f12.6)



151   continue !scan dky

	close(14)


888   write(88,*),886
	end

c--------------------------------------------------------------------
      SUBROUTINE eigen (n,a,val,vec)
      parameter (ny=32,nx=32,n3=nx)                
      PARAMETER (NMAX=n3,LDA=NMAX)
c      PARAMETER (LWORK=NMAX*NMAX+2*NMAX,LIWORK=3+5*NMAX,
c     +  LRWORK=2*NMAX*NMAX+5*NMAX+1)
      PARAMETER (LWORK=2*NMAX*NMAX+3*NMAX,LIWORK=6+8*NMAX,
     +  LRWORK=4*NMAX*NMAX+5*NMAX+2)

      CHARACTER JOB,UPLO
      complex vec(n,n)
      real val(n)

      complex A(LDA,NMAX),WORK(LWORK)
      real     RWORK(LRWORK)
      INTEGER IWORK(LIWORK),INFO

c	write(*,*)'before'
      N=n3
      UPLO='U'
      JOB='V'
      IF (UPLO.EQ.'U') THEN
      END IF
      CALL cheevd (JOB,UPLO,N,A,LDA,val,WORK,LWORK,RWORK,LRWORK,IWORK,
     +LIWORK,INFO)


      IF (INFO.GT.0)THEN
      WRITE (*,*) 'Failure to converge'
      ENDIF


      do 3 i=1,n
      do 3 j=1,n
    3  vec(i,j)=A(i,j)
      return
      end

附件也有代码文件




32.f

2.56 KB, 下载次数: 0

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2

帖子

1

主题

0

精华

新人

F 币
27 元
贡献
9 点
沙发
 楼主| 发表于 2021-10-19 18:19:16 | 只看该作者
已经解决了,代码没有问题,是自己操作失误了。
(不知道怎么删帖,就说明一下)
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-1 07:30

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表