myprojects-hg-reborn

Fork

program excitations

implicit none

real(kind=8), allocatable :: L1(:,:),L0(:,:),Marr(:,:),ini(:),diag(:),nondiag(:)
real(kind=8) :: epsi,mass,sigma,norm,lambda
integer:: Nsites,i, unitout,j,k,Ntot
character*12:: name2
! parameters for the lapack routine
character*1 :: jobvl,jobvr
integer :: info,lwork
real(kind=8), allocatable :: wr(:),wi(:),vr(:,:),vl(:,:),work(:),temp(:)
external :: DGEEV,saveresultx


Nsites=501

mass=dble(9.0d0)
print*,"The mass is ",mass


epsi=dble(1.0d0)
print*,"epsi is ",epsi
lambda=dble(6.99713d0)

print*,"Lambda is ", lambda

sigma=dble(1.0d0)
print*,"sigma is ",sigma
unitout=18



allocate(L1(Nsites,Nsites),L0(Nsites,Nsites),Marr(2*Nsites,2*Nsites),ini(Nsites))
allocate(diag(Nsites),nondiag(Nsites-1))

L0=dble(0.0d0)    !I initialize the (sparse) building blocks of Marr
L1=dble(0.0d0)
Marr=dble(0.0d0)  ! and Marr itself

!Now I read ini (normalized to one). Phi is the ground state I determine using another code.

open(18,file="ground-state",status='unknown')
 
                    do i=1,Nsites
                       read(18,*) ini(i)
                    end do
              close(18)




!as a test, I work out the initial normalization of the state and I tehn renormalize it to one.

norm=dble(0.0d0)

do i=1,Nsites

norm=norm+ini(i)*ini(i)

enddo

print*,"norm is", norm

!Now I make sure that the ground state is normalized to one

ini=ini/dsqrt(norm)

 do i=1,Nsites


     diag(i)=dble(2.0d0)*epsi-((dabs(ini(i)))**(dble(2.0d0)*sigma))*mass + lambda
  enddo



do i=1,Nsites-1

nondiag(i)=-epsi
enddo              

! Now I can build L0 and L1

!First L0

do i=1,Nsites

L0(i,i)=diag(i)

enddo


do i=1,Nsites-1

L0(i,i+1)=nondiag(i)
L0(i+1,i)=nondiag(i)

enddo
print*,"L0(11,12) and L0(12,11) and L0(12,14) are ",L0(11,12),"and ",L0(12,11), "and ",L0(12,14)


!Now I deal with L1

do i=1,Nsites


     diag(i)=dble(2.0d0)*epsi-dble(3.0d0)*((dabs(ini(i)))**(dble(2.0d0)*sigma))*mass + lambda
  enddo


do i=1,Nsites

L1(i,i)=diag(i)

enddo


do i=1,Nsites-1

L1(i,i+1)=nondiag(i)
L1(i+1,i)=nondiag(i)

enddo


! Now I can build up my matrix Marr

do i=1,Nsites
do j=1,Nsites

Marr(i,j+Nsites)=L0(i,j)

enddo
enddo


! and the other block is


! Now I can build up my matrix Marr

do i=1,Nsites
do j=1,Nsites

Marr(i+Nsites,j)=L1(i,j)

enddo
enddo


print*,"Marr(1,1), Marr(2,2), Marr(3,3) are ",Marr(1,1),Marr(2,2),Marr(3,3)
print*,"Marr(151,151), Marr(152,152), Marr(52,152) are ",Marr(151,151),Marr(152,152),Marr(52,152)
print*,"Marr(1,102), Marr(2,103), Marr(51,152) are ",Marr(1,102),Marr(2,103),Marr(51,152)

! Now I can set up the diagonalization of this non-symmetric matrix


!character*1 :: jobvl,jobvr
!integer :: lda,ldvl,ldvr,info,lwork
!real(kind=8), allocatable :: wr(:),wi(:),vr(:,:),vl(:,:),work(:)

jobvl="N"
jobvr="V"
Ntot=2*Nsites
allocate(wr(Ntot),wi(Ntot),vl(Ntot,Ntot),vr(Ntot,Ntot))
lwork=4*Ntot
allocate(work(lwork),temp(Ntot))
! SUBROUTINE DGEEV( JOBVL, JOBVR, N, A,	LDA, WR, WI, VL, LDVL, VR, LDVR,
!		    WORK, LWORK, INFO )

call DGEEV( JOBVL, JOBVR, Ntot, Marr,Ntot, WR, WI, VL,Ntot, VR, Ntot,  work, lwork, info )



print*,"Info is ",info



call saveresultx ("real--eigval",unitout,wr,Ntot)
call saveresultx ("imag--eigval",unitout,wi,Ntot)

do i=1,Ntot

do k=1,Ntot

temp(k)=vr(k,i)

enddo

write (name2, '("eigenvect",I3.3)')i
          call saveresultx (name2,unitout,temp,Ntot)

enddo




print*,"So far so good"

end program excitations




subroutine saveresultx (name,unitout,data,NN)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      implicit none


      integer :: K, unitout, NN
      real(kind=8) :: Data(NN)
      character*12 :: name


!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      open (unitout,file=name)
      print *,'Result saved in file:',name


      do K=1,NN
               if (dabs(Data(K)) .LE. 1.0D-99) Data(K)=0.0D0
               write (unitout,100) Data(K)
      end do


      close (unitout)  
      return


100   format (1E15.7E2)  ! last E2 determines two digits in Exponent
                         !  to avoid confusion in gnuplot
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      end subroutine saveresultx   
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC