!--------------!
 module system
!--------------!
 save

 integer :: nn
 integer :: nb
 integer :: zz

 integer, allocatable :: site(:,:)
 integer, allocatable :: hamx(:)
 integer, allocatable :: nxop(:)
 real(8), allocatable :: hamz(:)
 real(8), allocatable :: jbnd(:)

 complex(8), allocatable :: ff(:)
 complex(8), allocatable :: f0(:)
 complex(8), allocatable :: f1(:)
 complex(8), allocatable :: f2(:)

 real(8), allocatable :: sdat(:,:)

 end module system
!-----------------!

!--------------------!
 program timevolution
!--------------------!
 use system; implicit none

 integer :: i,j,t,nt,wf,a0,bins,confs
 real(8) :: pn,hh,jj,dh1,dh2,dj1,dj2,tt,dt

 open(10,file='read.in',status='old')
 read(10,*)nn,pn
 read(10,*)tt,dt,wf
 read(10,*)bins,confs
 close(10)

! nn = number of spins
! pn = probability of negative couplings
! tt = integration time
! dt = integration time step
! wf = output results every wf time step
! bins  = number of bins
! confs = number of disorder configs per bin
 
 call initran(1)

 call lattice()
 call makehamiltonian()

 nt=int(tt/dt+0.1d0)    ! total number of RK integration steps
 allocate(sdat(2,nt/wf)); sdat=0.d0

 do j=1,bins
    do i=1,confs
       call couplings(pn,a0)
       call initstate()   
       do t=0,nt-1
          call ramp(t,nt,hh,jj,dh1,dh2,dj1,dj2)
          call timestep(dt,hh,jj,dh1,dh2,dj1,dj2)               
          if (mod(t+1,wf)==0) call savedata((t+1)/wf,a0)
       enddo
    enddo
    call writedata(nt/wf,confs)
 enddo

 end program timevolution
!------------------------!

!-------------------------!
 subroutine savedata(t,a0)
!-------------------------!
 use system; implicit none

 integer :: a,t,a0
 real(8) :: ez,p0

 sdat(1,t)=sdat(1,t)+2.d0*abs(ff(a0))**2
 sdat(2,t)=sdat(2,t)+sum(hamz(:)*abs(ff(:))**2)-hamz(a0)

 end subroutine savedata
!-----------------------!

!------------------------------!
 subroutine writedata(nt,confs)
!------------------------------!
 use system; implicit none

 integer :: i,nt,confs

 sdat=sdat/dble(confs)
 open(10,file='res.dat',position='append')
 do i=1,nt
    write(10,'(i6,2f18.12)')i,sdat(:,i)
 enddo 
 close(10)
 sdat=0.d0

 end subroutine writedata
!------------------------!

!-------------------------------------------!
 subroutine ramp(t,nt,hh,jj,dh1,dh2,dj1,dj2)
!-------------------------------------------!
 use system; implicit none

 integer :: t,nt
 real(8) :: ss,ds,hh,jj,dh1,dh2,dj1,dj2

 ss=dble(t)/dble(nt)
 jj=ss
 hh=1.d0-ss
 ds=1.d0/dble(nt)    
 dj1=0.5d0*ds
 dh1=-dj1
 dj2=ds
 dh2=-dj2

 end subroutine ramp
!-------------------!

!---------------------------------------------!
 subroutine timestep(dt,hh,jj,dh1,dh2,dj1,dj2)
!------------------------------------------------------------------------------!
! Uses the RK method to integrate the Schrodinger equation by one time step dt !
!------------------------------------------------------------------------------!
 use system; implicit none

 integer :: i,j
 real(8) :: hh,jj,dh1,dh2,dj1,dj2,dt

 f2=ff
 f0=ff
 call hoperation(hh,jj)
 f1=-dt*(0,1)*f1
 f2=f2+f1/6.d0
 f0=ff+0.5d0*f1
 call hoperation(hh+dh1,jj+dj1)
 f1=-dt*(0,1)*f1
 f2=f2+f1/3.d0
 f0=ff+0.5d0*f1
 call hoperation(hh+dh1,jj+dj1)
 f1=-dt*(0,1)*f1
 f2=f2+f1/3.d0
 f0=ff+f1
 call hoperation(hh+dh2,jj+dj2)
 f1=-dt*(0,1)*f1
 f2=f2+f1/6.d0
 ff=f2/sqrt(dot_product(f2,f2))

 ! Explicit symmetrization of the state
 do i=0,2**nn-1
    j=ieor(i,zz)
    if (i<j) ff(j)=ff(i)
 enddo

 end subroutine timestep
!-----------------------!

!----------------------------!
 subroutine hoperation(hh,jj)
!----------------------------!
 use system; implicit none

 integer :: i,j,a,b,c
 real(8) :: hh,jj

 f1(:)=jj*hamz(:)*f0(:)
 j=0
 do a=0,2**nn-1
    f1(hamx(j+1:j+nxop(a)))=f1(hamx(j+1:j+nxop(a)))-hh*f0(a)
    j=j+nxop(a)
 enddo

! Take care of spin-reflection related states not included above, 
! using spin-reflection symmetry

 do b=0,2**nn-1           
    c=ieor(b,zz)          
    if (b<c) f1(c)=f1(b)  
 enddo                    

 end subroutine hoperation
!-------------------------!

!----------------------------!
 Subroutine makehamiltonian()
!-------------------------------------------------------!
! allocates the storage arrays for the hamiltoniand and !
! makes the list of the non-zero off-diagonal elements  !
! (their locations in the matrix)                       !
!-------------------------------------------------------!
 use system; implicit none

 integer :: i,j,a,b,c,hz

 allocate(nxop(0:2**nn-1))    ! number of off-diagonal elements for each basis state
 allocate(hamx(nn*2**(nn-1))) ! location of off-diagonal H elements compact-stored here
 allocate(hamz(0:2**nn-1))    ! diagonal H elements will be stored here

! zz is the number used to spin-reflect a state with ieor().
 zz=0
 do i=0,nn-1
    zz=ibset(zz,i)
 enddo

 j=0
 do a=0,2**nn-1
    nxop(a)=0
    do i=0,nn-1
       if (btest(a,i)) then
          b=ibclr(a,i)
       else
          b=ibset(a,i)
       endif

      ! Spin-reflect the state b, giving c.
      ! Only include Hamiltonian term if b<c.
       c=ieor(b,zz)
       if (b<c) then
          j=j+1
          hamx(j)=b
          nxop(a)=nxop(a)+1
       endif
    enddo

    ! Sort the Hamiltonian terms, can speed up the execution
    if (nxop(a)>1) call sort(nxop(a),hamx(j-nxop(a)+1:j))

 enddo

 end subroutine makehamiltonian
!------------------------------!

!--------------------!                                                                                                        
 subroutine sort(n,a)
!----------------------------!
! Based on Numerical Recipes.
!----------------------------!
 implicit none

 integer :: i,j,l,n,ir,rra,a(n)

 l=n/2+1
 ir=n
 do
    if (l>1) then
       l=l-1
          rra=a(l)
       else
       rra=a(ir)
       a(ir)=a(1)
       ir=ir-1
       if (ir==1) then
          a(1)=rra
          return
       endif
    endif
    i=l
    j=l+l
    do
       if (j<=ir) then
          if (j<ir) then
             if (a(j)<a(j+1)) j=j+1
          endif
          if (rra<a(j)) then
             a(i)=a(j)
             i=j
             j=j+j
          else
             j=ir+1
          endif
       else
          exit
       endif
    enddo
    a(i)=rra
 enddo

 end subroutine sort
!-------------------!                                                                                                            

!----------------------!
 subroutine initstate()
!----------------------!
 use system; implicit none

 if (.not.allocated(ff)) then
    allocate(ff(0:2**nn-1))
    allocate(f0(0:2**nn-1))
    allocate(f1(0:2**nn-1)) 
    allocate(f2(0:2**nn-1)) 
 endif

 ff(:)=1.d0
 ff(:)=ff(:)/sqrt(dot_product(ff(:),ff(:)))
 
 end subroutine initstate
!------------------------!

!--------------------!
 subroutine lattice()
!------------------------------------!
! generates the lattice connectivity !
! (sites connected by interactions)  !
!------------------------------------!
 use system; implicit none

 integer :: i,j,b

 nb=(nn*(nn-1))/2
 allocate(jbnd(nb))
 allocate(site(2,nb))
 b=0
 do i=0,nn-1
    do j=i+1,nn-1
       b=b+1
       site(1,b)=i
       site(2,b)=j
    enddo
 enddo

 end subroutine lattice
!----------------------!

!---------------------------!
 subroutine couplings(pn,a0)
!---------------------------!
 use system; implicit none

 integer :: i,a,b,a0
 real(8) :: pn,hz

 real(8), external :: ran

 do i=1,nb
    jbnd(i)=ran()
    if (ran()<pn) jbnd(i)=-jbnd(i)
 enddo

 do a=0,2**nn-1
    hz=0.d0
    do i=1,nb
       if (btest(a,site(1,i)).eqv.btest(a,site(2,i))) then
          hz=hz-jbnd(i)
       else
          hz=hz+jbnd(i)
       endif
    enddo
    hamz(a)=hz/(dble(nn-1)/2.d0)
 enddo

 hz=0.d0
 do a=0,2**nn-1
    b=ieor(a,zz)
    if (a<b.and.hamz(a)<hz) then
       hz=hamz(a)
       a0=a
    endif
 enddo

 end subroutine couplings
!------------------------!

!--------------------!
 subroutine cleanup()
!--------------------!
 use system; implicit none

 deallocate(site)
 deallocate(jbnd)
 deallocate(hamz)
 deallocate(hamx)
 deallocate(nxop)
 deallocate(ff)
 deallocate(f0)
 deallocate(f1)
 deallocate(f2)

 end subroutine cleanup
!----------------------!

!----------------------!
 real(8) function ran()
!----------------------------------------------!
! 64-bit congruental generator                 !
! iran64=oran64*2862933555777941757+1013904243 !
!----------------------------------------------!
 implicit none

 real(8)    :: dmu64
 integer(8) :: ran64,mul64,add64
 common/bran64/dmu64,ran64,mul64,add64

 ran64=ran64*mul64+add64
 ran=0.5d0+dmu64*dble(ran64)

 end function ran
!----------------!

!---------------------!
 subroutine initran(w)
!---------------------!
 implicit none

 integer(8) :: irmax
 integer(4) :: w,nb,b

 real(8)    :: dmu64
 integer(8) :: ran64,mul64,add64
 common/bran64/dmu64,ran64,mul64,add64
      
 irmax=2_8**31
 irmax=2*(irmax**2-1)+1
 mul64=2862933555777941757_8
 add64=1013904243
 dmu64=0.5d0/dble(irmax)

 open(10,file='seed.in',status='old')
 read(10,*)ran64
 close(10)
 if (w.ne.0) then
    open(10,file='seed.in',status='unknown')
    write(10,*)abs((ran64*mul64)/5+5265361)
    close(10)
 endif

 end subroutine initran
!----------------------!