module aa
   integer n,m
   parameter (n=100,m=n/2)
   real t,eps
   INTEGER LDA,LDAINV,LDEVEC
   parameter (LDA=n,LDAINV=n,LDEVEC=n,t=1.0d0,eps=0.001d0)
   complex cj
      parameter(cj=(0.,1.))
   end module

 !  *************************************************************
      PROGRAM MAIN

! ***********************************

 
   use imsl
   use aa
     implicit none

 
      INTEGER NT
      INTEGER NOUT
      COMPLEX(8),ALLOCATABLE,DIMENSION (:) :: H(:,:),H0(:,:),Ht(:,:),G(:,:),EVEC(:,:),state1(:,:),state2(:,:)
   COMPLEX(8),ALLOCATABLE,DIMENSION (:) :: S(:,:),A(:,:),AINV(:,:)
      REAL(8),ALLOCATABLE,DIMENSION (:) :: energy1(:),energy2(:),EVAL(:),rou1(:),rou2(:),rou(:)
   integer i,j,k,l,p,q
   real(8) z2,omega2,eta,pi,sin,cos,faik,theta,eta1,eta2,time,temp
   real(8) delta1,delta2,delta,mu,mu1,mu2,pr
   complex(8) pr2
   real(8) tanh,tanh2,rout,tanh3,eta3

      ALLOCATE(energy1(n))
   ALLOCATE(energy2(n))
      ALLOCATE(H(n,n))
   ALLOCATE(H0(n,n))
   ALLOCATE(Ht(n,n))
   ALLOCATE(S(n,n))
      ALLOCATE(EVEC(n,n))
   ALLOCATE(state1(n,n))
   ALLOCATE(state2(n,n))
      ALLOCATE(EVAL(n))
   ALLOCATE(G(n,n))
   ALLOCATE(A(n,n))
   ALLOCATE(AINV(n,n))
   ALLOCATE(rou(n/2))
   ALLOCATE(rou1(n/2))
   ALLOCATE(rou2(n/2))

   pi=4.0d0*datan(1.0d0)


   open(2,file='energy2')

 eta1=0.0d0
 eta2=0.0d0
 mu1=0.0d0
 mu2=0.0d0
 delta1=0.1d0
 delta2=0.0d0


 delta=delta1
 eta=eta1
 mu=mu1

 H=0.0d0
 Ht=0.0d0
 energy1=0.0
 energy2=0.0
 state1=0.0
 state2=0.0

!    do delta=-1.0d0,1.0d0,0.1d0
 do i=1,m-1,1

    H(2*i-1,2*i-1)=mu   !alphaj_dagger
 H(2*i,2*i)=-mu   !alphaj

 H(2*i-1,2*i+1)=-t
    H(2*i+1,2*i-1)=-t
 H(2*i+2,2*i)=t
    H(2*i,2*i+2)=t

 H(2*i,2*i+1)=delta
 H(2*i+1,2*i)=delta
 H(2*i-1,2*i+2)=-delta
 H(2*i+2,2*i-1)=-delta
 enddo

 i=m

    H(2*i-1,2*i-1)=mu   !alphaj_dagger
 H(2*i,2*i)=-mu   !alphaj


 H=H*0.25d0

 H0=H


    CALL DEVCHF(n, H, LDA, EVAL, EVEC, LDEVEC)
 CALL ARRANGE(n,EVAL,EVEC)

 energy1=EVAL
    state1=EVEC
 delta=delta2
 mu=mu2
 eta=eta2

 do i=1,m-1,1

    H(2*i-1,2*i-1)=mu   !alphaj_dagger
 H(2*i,2*i)=-mu   !alphaj

 H(2*i-1,2*i+1)=-t
    H(2*i+1,2*i-1)=-t
 H(2*i+2,2*i)=t
    H(2*i,2*i+2)=t

 H(2*i,2*i+1)=delta
 H(2*i+1,2*i)=delta
 H(2*i-1,2*i+2)=-delta
 H(2*i+2,2*i-1)=-delta
 enddo

 i=m

    H(2*i-1,2*i-1)=mu   !alphaj_dagger
 H(2*i,2*i)=-mu   !alphaj

 H=H*0.25d0

    CALL DEVCHF(n, H, LDA, EVAL, EVEC, LDEVEC)
 CALL ARRANGE(n,EVAL,EVEC)


 energy2=EVAL
    state2=EVEC

    time=25.0d0
 Ht=0.0d0
 S=0.0d0

 do i=1,n,1
 do j=1,n,1
 do k=1,n,1

    S(i,j)=S(i,j)+CONJG(state1(k,i))*state2(k,j)

 enddo
 enddo
 enddo


 do i=1,n,1
    write(*,*) i
 do j=1,n,1
 do k=1,n,1
 do l=1,n,1
! pause

    Ht(i,j)=Ht(i,j)+exp(cj*(energy2(i)-energy2(j))*time)*CONJG(S(k,i))*H0(k,l)*S(l,j)

 enddo

 enddo
 enddo
 enddo

 do i=1,n,1
 Ht(i,i)=Ht(i,i)-cj*eps

 enddo 

    A=Ht*(-1.0d0)
 CALL DLINCG (n,A,LDA,AINV,LDAINV)
 G=AINV


 do i=1,m,1
 rou1(i)=-1.0d0/pi*AIMAG(G(2*i-1,2*i-1))
 rou2(i)=-1.0d0/pi*AIMAG(G(2*i,2*i))

 rou(i)=rou1(i)+rou2(i)

 rout=rout+rou(i)
 enddo

     
   do i=1,n/2,1
   write(*,*) i,rou(i)
  !    pause

   WRITE(2,100) real(i),rou1(i),rou2(i),rou(i)

 enddo
! enddo !time


100   FORMAT(1X,1290F20.10) 

 !   enddo  ! eta loop

 DEALLOCATE(energy1)
 DEALLOCATE(energy2)
    DEALLOCATE(H)
    DEALLOCATE(EVEC)
 DEALLOCATE(state1)
 DEALLOCATE(state2)
    DEALLOCATE(EVAL)


! enddo !psaik loop
!enddo      !N loop
!enddo       !gama loop
  
   end

      SUBROUTINE ARRANGE(N,ARRAY,evec)
      IMPLICIT NONE
      INTEGER N,I,J,k
      REAL(8) ARRAY(N),AA
      COMPLEX(8) evec(n,n),temp(n)
      DO 12 I=2,N
      AA=ARRAY(I)
      do k=1,n
      temp(k)=evec(k,i)
      enddo
      DO 11 J=I-1,1,-1
      IF(ARRAY(J)<=AA)GOTO 15 
      ARRAY(J+1)=ARRAY(J)
      do k=1,n
      evec(k,j+1)=evec(k,j)
      enddo
   11 CONTINUE
      J=0
   15 ARRAY(J+1)=AA
      do k=1,n
      evec(k,j+1)=temp(k)
      enddo
   12 CONTINUE
      RETURN
      END