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

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