program jastrow
      implicit none 
      INTEGER(4) nh,sztot
      INTEGER(4) ngen,nscra,nbra,ncorr
      INTEGER(4) i,j,jn,iout,jout,indvic
      INTEGER(4) nacc
      REAL(8) jperp,ener,alpha
      REAL(8) ratio,zeta,rata,rnd

      INTEGER(4), dimension(8) :: iseed
      INTEGER(4), dimension(:,:), allocatable :: ivic
      INTEGER(4), dimension(:), allocatable :: iconf

      REAL(8), dimension(:), allocatable :: tabpip
      REAL(8), dimension(:), allocatable :: szsz
      REAL(8), dimension(:,:), allocatable :: vpot,vjas

      namelist /lattice/ nh
      namelist /parameters/ jperp,sztot
      namelist /wavefunction/ alpha
      namelist /montecarlo/ iseed,ngen,nbra,nscra,ncorr

! reading part
      nh=0
      jperp=1.d0
      sztot=0
      alpha=1.d0

      ngen=0
      nbra=0
      nscra=0
      ncorr=0

      read(5,lattice)
      read(5,parameters)
      read(5,wavefunction)

      if(nh==0) then
       write(6,*) 'nh must be specified'
       stop
      endif

      write(6,*) ' Number of sites                    :',nh
      write(6,*) ' Total Sz                           :',sztot
      write(6,*) ' Anisotropy of the super-exchange   :',jperp
      write(6,*) ' Variational parameter alpha        :',alpha

      read(5,montecarlo)

      if(ngen==0) then
       write(6,*) 'ngen must be specified'
       stop
      endif
      if(nbra==0) then
       write(6,*) 'nbra must be specified'
       stop
      endif
      if(nscra==0) then
       write(6,*) 'nscra must be specified'
       stop
      endif
      write(6,*)
      write(6,*) ' Number of measures                 :',ngen
      write(6,*) ' Number of MC steps between measures:',nbra
      write(6,*) ' Number of measures between upscra  :',nscra
      if(ncorr==0) then
       write(6,*) ' No correlation functions computed'
      else
       write(6,*) ' Correlation functions computed'
      endif

      call random_seed(put=iseed)

      open(unit=11,file='fort.11',form='formatted',status='unknown')
      open(unit=12,file='fort.12',form='unformatted',status='unknown')
      if(ncorr/=0) then
       open(unit=13,file='fort.13',form='unformatted',status='unknown')
      endif

      rewind(11)
      rewind(12)
      if(ncorr/=0) rewind(13)

      ALLOCATE(ivic(nh,2))
      ALLOCATE(vpot(nh,nh))
      ALLOCATE(vjas(nh,nh))
      ALLOCATE(szsz(nh))
      ALLOCATE(tabpip(nh))
      ALLOCATE(iconf(nh))

! table of nearest neighbors
      call neighbors(nh,ivic)

! pseudo-potential of the Jastrow
      call pseudo(nh,alpha,vpot,vjas)

! random initialization of spins
      call init(nh,sztot,iconf)  

! main VMC loop
      nacc=0

      do i=1,ngen

       if(mod(i,nscra)==1) then
        call upscratch(nh,iconf,vjas,tabpip)
       endif

       do j=1,nbra

! nearest-neighbor spin flip
        call random_number(rnd)
        iout=rnd*nh+1
        call random_number(rnd)
        indvic=rnd*2+1
        jout=ivic(iout,indvic)

        call ratiovar(iout,jout,nh,iconf,tabpip,vjas,ratio)

        call random_number(rnd)
        zeta=1.d0-rnd
        if(ratio**2>zeta) then
         nacc=nacc+1
         call upjastrow(nh,iout,jout,iconf,tabpip,vjas)
        endif

       enddo

       call localenergy(nh,iconf,ivic,tabpip,vjas,jperp,ener)

       write(11,*) i,ener/nh
       write(12) i,ener/nh

       if(ncorr/=0) then
        call spinspinz(nh,iconf,szsz)
        write(13) i,(szsz(j),j=1,nh)
       endif

      enddo

      rata=dble(nacc)/(ngen*nbra)
      write(6,*) 
      write(6,*) 'accept. rate off diagonal moves =',rata

      close(11)
      close(12)
      if(ncorr/=0) close(13)

      DEALLOCATE(ivic)
      DEALLOCATE(vpot)
      DEALLOCATE(vjas)
      DEALLOCATE(szsz)
      DEALLOCATE(tabpip)
      DEALLOCATE(iconf)

      stop 
      end

!======================================================================

      subroutine neighbors(nh,ivic)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,il,ir
      INTEGER(4) ivic(nh,2)

      do i=1,nh
       ir=0
       if(i==nh) ir=nh
       il=0
       if(i==1) il=nh
       ivic(i,1)=i+1-ir        ! right
       ivic(i,2)=i-1+il        ! left
      enddo

      return
      end

      subroutine pseudo(nh,alpha,vpot,vjas)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,j
      REAL(8) pi,dist
      REAL(8) alpha
      REAL(8) vpot(nh,nh),vjas(nh,nh)

      pi=dacos(-1.d0)

      do i=1,nh
       do j=1,nh
        dist=nh/pi*dsin(pi*abs(i-j)/nh)
        vpot(i,j)=0.d0
        if(i/=j) vpot(i,j)=dlog(dist**2)
       enddo
      enddo

      do i=1,nh
       do j=1,nh
        vjas(i,j)=dexp(alpha*vpot(i,j))
       enddo
      enddo

      return
      end

      subroutine init(nh,sztot,iconf)
      implicit none
      INTEGER(4) nh,sztot,nup
      INTEGER(4) i,jr,ib
      INTEGER(4) iconf(nh)
      REAL(8) rnd

      nup=nh/2+sztot 

      do i=1,nh
       iconf(i)=-1
      enddo

      ib=0
      do while(ib<nup) 
       call random_number(rnd)
       jr=rnd*nh+1
       if(iconf(jr)==-1) then
        iconf(jr)=1
        ib=ib+1
       endif
      enddo

      return
      end

      subroutine upscratch(nh,iconf,vjas,tabpip)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,j
      INTEGER(4) iconf(nh)
      REAL(8) tabpip(nh),vjas(nh,nh)

! tabpip(i)=exp[\sum_k v_{i,k} (S^z_k+1/2)]
! J=exp[1/2 \sum_{i,k} v_{i,k} (S^z_i+1/2)(S^z_k+1/2)]

      do i=1,nh
       tabpip(i)=1.d0
      enddo

      do j=1,nh
       if(iconf(j)==1) then
        do i=1,nh
         tabpip(i)=tabpip(i)*vjas(i,j)            
        enddo
       endif
      enddo

      return
      end
 
      subroutine ratiovar(i,j,nh,iconf,tabpip,vjas,ratio)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,j
      REAL(8) ratio
      INTEGER(4) iconf(nh)
      REAL(8) tabpip(nh),vjas(nh,nh)

      if(iconf(i)==-1.and.iconf(j)==1) then
       ratio=tabpip(i)/tabpip(j)*vjas(i,i)/vjas(i,j)
      elseif(iconf(i)==1.and.iconf(j)==-1) then
       ratio=tabpip(j)/tabpip(i)*vjas(i,i)/vjas(i,j)
      else
       ratio=0.d0
      endif

      return
      end

      subroutine localenergy(nh,iconf,ivic,tabpip,vjas,jperp,ener)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,jn
      REAL(8) jperp,ener
      INTEGER(4) iconf(nh),ivic(nh,2)
      REAL(8) tabpip(nh),vjas(nh,nh)

      ener=0.d0

      do i=1,nh
       jn=ivic(i,1)
       ener=ener+0.25*jperp*iconf(i)*iconf(jn)
      enddo

! the minus sign is due to the Marshall sign rule
      do i=1,nh
       jn=ivic(i,1)
       if(iconf(i)==-1.and.iconf(jn)==1) then
        ener=ener-0.5d0*tabpip(i)/tabpip(jn)*vjas(i,i)/vjas(i,jn)
       elseif(iconf(i)==1.and.iconf(jn)==-1) then
        ener=ener-0.5d0*tabpip(jn)/tabpip(i)*vjas(i,i)/vjas(i,jn)
       endif
      enddo

      return
      end

      subroutine spinspinz(nh,iconf,szsz)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,j,k
      INTEGER(4) iconf(nh)
      REAL(8) szsz(nh)

      do i=1,nh
       szsz(i)=0.d0
       do j=1,nh
        k=mod(j+i-1,nh)+1
        if(iconf(j)==iconf(k)) then
         szsz(i)=szsz(i)+0.25d0
        else
         szsz(i)=szsz(i)-0.25d0
        endif
       enddo
       szsz(i)=szsz(i)/nh
      enddo

      return
      end

      subroutine upjastrow(nh,iout,jout,iconf,tabpip,vjas)
      implicit none
      INTEGER(4) nh
      INTEGER(4) i,jn,iout,jout
      INTEGER(4) iconf(nh)
      REAL(8) tabpip(nh),vjas(nh,nh)
 
      if(iconf(iout)==1.and.iconf(jout)==-1) then
       do i=1,nh
        tabpip(i)=tabpip(i)/vjas(i,iout)
        tabpip(i)=tabpip(i)*vjas(i,jout)
       enddo
       iconf(iout)=-1
       iconf(jout)=1
      elseif(iconf(iout)==-1.and.iconf(jout)==1) then
       do i=1,nh
        tabpip(i)=tabpip(i)*vjas(i,iout)
        tabpip(i)=tabpip(i)/vjas(i,jout)
       enddo
       iconf(iout)=1
       iconf(jout)=-1
      else
       write(6,*) 'Wrong updating'
       stop
      endif

      return
      end