!--------------! module system !--------------! implicit none integer :: nn integer :: nrep integer :: smax integer, allocatable :: repr(:) integer, allocatable :: type(:) integer, allocatable :: peri(:) integer, allocatable :: mtrf(:) integer, allocatable :: ntrf(:) real(8), allocatable :: mat(:,:) real(8), allocatable :: vec(:,:) real(8), allocatable :: enr(:) real(8), allocatable :: spn(:) real(8), allocatable :: ffun(:,:,:) real(8), allocatable :: gfun(:,:) end module system !-----------------! !===================! program hchain_mkpz !===================! use system; implicit none integer :: k,p,z,p1,p2,rm open(10,file='read.in',status='old') read(10,*)nn,rm close(10) smax=2**nn-1 allocate(repr(rm)) allocate(type(rm)) allocate(peri(rm)) allocate(mtrf(rm)) allocate(ntrf(rm)) do k=0,nn/2 call fgfunctions(k) if (k==0.or.k==nn/2) then p1=-1; p2=+1 else p1=+1; p2=+1 endif do p=p2,p1,-2 do z=+1,-1,-2 call makebasis(k,p,z) if (nrep/=0) then allocate(mat(nrep,nrep)) allocate(vec(nrep,nrep)) allocate(enr(nrep)) allocate(spn(nrep)) call hamiltonian(p,z) call diagonalize(nrep,mat,vec,enr) call spinsquared(p,z) call transform(nrep,mat,vec,spn) spn(:)=0.5d0*abs(sqrt(1.d0+4.d0*spn(:))-1.d0) call writedata(k,p,z) deallocate(mat) deallocate(vec) deallocate(enr) deallocate(spn) endif enddo enddo enddo deallocate(repr) deallocate(type) deallocate(peri) deallocate(mtrf) deallocate(ntrf) deallocate(ffun) deallocate(gfun) end program hchain_mkpz !=======================! !---------------------------! subroutine writedata(k,p,z) !---------------------------! use system; implicit none integer :: i,k,p,z open(10,file='eig.dat',position='append') write(10,10)k,p,z,nrep 10 format(3i5,i10) do i=1,nrep write(10,20)i,enr(i),spn(i) 20 format(i5,' ',2f16.10) enddo close(10) open(10,file='low.dat',position='append') write(10,30)k,p,z,enr(1),spn(1),nrep 30 format(3i5,2f16.10,i10) close(10) end subroutine writedata !------------------------! !---------------------------! subroutine hamiltonian(p,z) !---------------------------! use system; implicit none integer :: i,j,a,b,p,z,l,q,g,ia,ib,sa,sb,na,nb real(8) :: matelement,ez mat(:,:)=0.d0 do ia=1,nrep sa=repr(ia) if (ia>1.and.sa==repr(ia-1)) then cycle elseif (ia0) then if (ib>1.and.repr(ib)==repr(ib-1)) then ib=ib-1 nb=2 elseif (ib1.and.sa==repr(ia-1)) then cycle elseif (ia0) then if (ib>1.and.repr(ib)==repr(ib-1)) then ib=ib-1 nb=2 elseif (ib1.d-8) pass=.false. elseif (tp==nn.and.tz/=nn) then ca=3; n=tz if (gfun(z,n)<1.d-8) pass=.false. elseif (tp==nn.and.tz==nn) then ca=4; m=tpz if (gfun(s*p*z,m)<1.d-8) pass=.false. if (s==-1.and.gfun(-s*p*z,m)>1.d-8) pass=.false. elseif (tp/=nn.and.tz/=nn) then ca=5; m=tp; n=tz if (gfun(z,n)*gfun(s*p,m)<1.d-8) pass=.false. if (s==-1.and.gfun(-s*p,m)>1.d-8) pass=.false. endif if (pass) then nrep=nrep+1 repr(nrep)=a type(nrep)=2*ca+(s+1)/2 peri(nrep)=ra mtrf(nrep)=m ntrf(nrep)=n endif pass=.true. enddo endif enddo end subroutine makebasis !------------------------! !---------------------------------------------! subroutine checkstate(a0,k,ra,tp,tz,tpz,pass) !---------------------------------------------! use system; implicit none integer :: i,k,t,m,a0,at,az,ra,tp,tz,tpz logical :: pass pass=.false. m=0 do i=0,nn-1 if (btest(a0,i)) m=m+1 enddo if (m/=nn/2) return ra=nn tz=nn at=a0 do t=1,nn-1 if (btest(at,0)) then at=at/2 at=ibset(at,nn-1) else at=at/2 endif az=smax-at if (atrepr(a)) then amin=a+1 else return endif if (amin>amax) then a=-1 return endif enddo end subroutine findstate !------------------------! !-------------------------! subroutine fgfunctions(k) !-------------------------! use system; implicit none real(8), parameter :: pi=3.14159265358979d0 integer :: i,k real(8) :: kk kk=dfloat(2*k)*pi/dfloat(nn) if (.not.allocated(ffun)) then allocate(ffun(-1:1,-1:1,-nn:nn)) allocate(gfun(-1:1,-nn:nn)) endif do i=-nn,nn ffun(-1,-1,i)=+cos(kk*i) ffun(+1,+1,i)=+cos(kk*i) ffun(+1,-1,i)=+sin(kk*i) ffun(-1,+1,i)=-sin(kk*i) gfun(-1,i)=1.d0-cos(kk*i) gfun(+1,i)=1.d0+cos(kk*i) enddo end subroutine fgfunctions !--------------------------! !-------------------------------------! subroutine diagonalize(n,mat,vec,eig) !-------------------------------------! implicit none integer :: n,info real(8) :: mat(n,n),vec(n,n),eig(n),work(n*(3+n/2)) vec=mat call dsyev('V','U',n,vec,n,eig,work,n*(3+n/2),info) end subroutine diagonalize !--------------------------! !-----------------------------------! subroutine transform(n,mat,vec,dia) !-----------------------------------! implicit none integer :: i,n real(8) :: mat(n,n),vec(n,n),dia(n) mat=matmul(mat,vec) mat=matmul(transpose(vec),mat) do i=1,n dia(i)=mat(i,i) enddo end subroutine transform !------------------------!