您所在位置：主题：求助！这个程序是干嘛的？计算物理方向
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