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

 integer :: lx
 integer :: ly
 integer :: nn
 integer :: nb
 integer :: nd
 integer :: mm

 integer, allocatable :: oper(:)
 integer, allocatable :: lbnd(:)
 integer, allocatable :: rbnd(:)
 integer, allocatable :: lspn(:)
 integer, allocatable :: rspn(:)
 integer, allocatable :: frst(:)
 integer, allocatable :: last(:)
 integer, allocatable :: vrtx(:)
 integer, allocatable :: bsite(:,:)
 integer, allocatable :: dsite(:,:)
 integer, allocatable :: xy(:,:)
 integer, allocatable :: xyi(:,:)
 real(8), allocatable :: hamp(:,:) 

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

!------------!
 module mdata
!------------!
 save

 integer :: nmsr=0.d0
 integer, allocatable :: plbnd(:)
 integer, allocatable :: prbnd(:)
 integer, allocatable :: loopnmbr(:)
 real(8), allocatable :: scor(:)

 end module mdata
!----------------!

!================!
 program probasic
!=====================!
! Anders Sandvik, 2012
!---------------------!
use system; use mdata; implicit none

 integer :: i,j,init,bins,mstps,istps

 open(10,file='read.in',status='old')
 read(10,*)lx,ly,mm
 read(10,*)init,bins,mstps,istps
 close(10)

 call initran(1)
 call makelattice()
 call initsystem(init)

 if (istps>0) then
    do i=1,istps
       call montecarlosweep()
     enddo
    call writeconf()
 endif

 allocate(scor(0:lx/2)); scor=0.d0

 do j=1,bins    
    do i=1,mstps
       call montecarlosweep()
       call measure() 
    enddo
    call writebindata()
    call writeconf()
 enddo
   
 call deallocateall()

 end program probasic
!====================!

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

 call diagonalupdate()
 call makevertexlist()
 call loopupdate()
 call stateupdate()

 end subroutine montecarlosweep
!------------------------------!

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

 integer :: i,b
 real(8) :: ran

 external :: ran

 lspn(:)=rspn(:)
 do i=0,mm-1
    if (mod(oper(i),2)==0) then       
       10 b=int(ran()*dble(nb))+1
       if (lspn(bsite(1,b))/=lspn(bsite(2,b))) then
          oper(i)=2*b
       else
          goto 10
       endif
    else
       b=oper(i)/2
       lspn(bsite(1,b))=-lspn(bsite(1,b))
       lspn(bsite(2,b))=-lspn(bsite(2,b))     
    endif
 enddo

 end subroutine diagonalupdate
!-----------------------------!

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

 integer :: i,b,p,s1,s2,v0,v1,v2,v3,v4

 last(:)=-1
 frst(:)=-1
 do p=0,mm-1
    v0=4*p
    b=oper(p)/2
    s1=bsite(1,b)
    s2=bsite(2,b)
    v1=last(s1)
    v2=last(s2)
    if (v1/=-1) then
       vrtx(v1)=v0
       vrtx(v0)=v1
    else
       frst(s1)=v0
    endif
    if (v2/=-1) then
       vrtx(v2)=v0+1
       vrtx(v0+1)=v2
    else
       frst(s2)=v0+1
    endif
    last(s1)=v0+2
    last(s2)=v0+3       
 enddo

 v0=4*mm
 do i=1,nn
    v1=v0+1; v2=v0+2; v3=v0+3
    if (frst(i)>=0) then
       vrtx(frst(i))=v0
       vrtx(v0)=frst(i)
       vrtx(last(i))=v2
       vrtx(v2)=last(i)
    else
       vrtx(v0)=v2
       vrtx(v2)=v0
    endif
    vrtx(v1)=4*mm+4*rbnd(i)-3
    vrtx(v3)=4*mm+4*lbnd(i)-1
    v0=v0+4
 enddo 

 end subroutine makevertexlist
!-----------------------------!

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

 integer :: b,p,v0,v1,v2
 real(8) :: ran
 external :: ran

 do v0=0,4*mm+4*nn-2,2
    if (vrtx(v0)>=0) then
       v1=v0
       if (ran()<0.5d0) then
          call visitloop()
       else
          call fliploop()
       endif
    endif
 enddo

 contains

   !------------------------!
    subroutine visitloop()
   !------------------------!

    do
       vrtx(v1)=-1
       v2=ieor(v1,1)
       v1=vrtx(v2)
       vrtx(v2)=-1
       if (v1==v0) exit
    enddo

    end subroutine visitloop
   !------------------------! 

   !---------------------!
    subroutine fliploop()
   !---------------------!

    do
       p=v1/4
       if (p<mm) then
          oper(p)=ieor(oper(p),1)
       else
          p=p-mm+1
          if (mod(v1,4)<2) then
             rspn(p)=-rspn(p)
          else
             lspn(p)=-lspn(p)
          endif
       endif
       vrtx(v1)=-1
       v2=ieor(v1,1)
       v1=vrtx(v2)
       vrtx(v2)=-1
       if (v1==v0) exit
    enddo

    end subroutine fliploop
   !-----------------------!

 end subroutine loopupdate
!-------------------------!

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

 integer :: i,d,x1,x2,y1,y2,s1,s2,s3,s4
 real(8) :: ran,w1,w2

 external :: ran

 do i=1,nn
    d=int(ran()*nd)+1
    s1=dsite(1,d)
    s2=dsite(2,d)
    s3=lbnd(s1)
    s4=lbnd(s2)
    if (lspn(s1)==lspn(s4)) cycle
    x1=abs(xy(1,s3)-xy(1,s1)); x1=min(x1,lx-x1)
    y1=abs(xy(2,s3)-xy(2,s1)); y1=min(y1,ly-y1)
    x2=abs(xy(1,s4)-xy(1,s2)); x2=min(x2,lx-x2)
    y2=abs(xy(2,s4)-xy(2,s2)); y2=min(y2,ly-y2)
    w1=hamp(x1,y1)*hamp(x2,y2)
    x1=abs(xy(1,s4)-xy(1,s1)); x1=min(x1,lx-x1)
    y1=abs(xy(2,s4)-xy(2,s1)); y1=min(y1,ly-y1)
    x2=abs(xy(1,s3)-xy(1,s2)); x2=min(x2,lx-x2)
    y2=abs(xy(2,s3)-xy(2,s2)); y2=min(y2,ly-y2)
    w2=hamp(x1,y1)*hamp(x2,y2)
    if (ran()<w2/w1) then
       lbnd(s1)=s4
       lbnd(s4)=s1
       lbnd(s2)=s3
       lbnd(s3)=s2
    endif
 enddo      
 do i=1,nn
    d=int(ran()*nd)+1
    s1=dsite(1,d)
    s2=dsite(2,d)
    s3=rbnd(s1)
    s4=rbnd(s2)
    if (rspn(s1)==rspn(s4)) cycle
    x1=abs(xy(1,s3)-xy(1,s1)); x1=min(x1,lx-x1)
    y1=abs(xy(2,s3)-xy(2,s1)); y1=min(y1,ly-y1)
    x2=abs(xy(1,s4)-xy(1,s2)); x2=min(x2,lx-x2)
    y2=abs(xy(2,s4)-xy(2,s2)); y2=min(y2,ly-y2)
    w1=hamp(x1,y1)*hamp(x2,y2)
    x1=abs(xy(1,s4)-xy(1,s1)); x1=min(x1,lx-x1)
    y1=abs(xy(2,s4)-xy(2,s1)); y1=min(y1,ly-y1)
    x2=abs(xy(1,s3)-xy(1,s2)); x2=min(x2,lx-x2)
    y2=abs(xy(2,s3)-xy(2,s2)); y2=min(y2,ly-y2)
    w2=hamp(x1,y1)*hamp(x2,y2)
    if (ran()<w2/w1) then
       rbnd(s1)=s4
       rbnd(s4)=s1
       rbnd(s2)=s3
       rbnd(s3)=s2
    endif
 enddo      

 end subroutine stateupdate
!--------------------------!

!---------------------------!
 subroutine propagatebonds()
!---------------------------!
 use system; use mdata; implicit none

 integer :: i,b,s1,s2,s3,s4

 prbnd(:)=rbnd(:)
 do i=0,mm/2-1
    b=oper(i)/2
    s1=bsite(1,b)
    s2=bsite(2,b)
    s3=prbnd(s1)      
    s4=prbnd(s2)
    prbnd(s1)=s2
    prbnd(s2)=s1
    prbnd(s3)=s4
    prbnd(s4)=s3
 enddo
 plbnd(:)=lbnd(:)
 do i=mm-1,mm/2,-1
    b=oper(i)/2
    s1=bsite(1,b)
    s2=bsite(2,b)
    s3=plbnd(s1)
    s4=plbnd(s2)
    plbnd(s1)=s2
    plbnd(s2)=s1
    plbnd(s3)=s4
    plbnd(s4)=s3
 enddo

 end subroutine propagatebonds
!-----------------------------!

!----------------------!
 subroutine makeloops()
!----------------------!
 use system; use mdata; implicit none

 integer :: i,s,s1,nl

 nl=0
 loopnmbr(:)=0
 do s1=1,nn
    if (loopnmbr(s1)==0) then    
       nl=nl+1
       s=s1
       do 
          loopnmbr(s)=nl
          s=prbnd(s)
          loopnmbr(s)=nl
          s=plbnd(s)
          if (s==s1) exit
       enddo
    endif
 enddo   

 end subroutine makeloops
!------------------------!

!--------------------!
 subroutine measure()
!--------------------!
 use system; use mdata; implicit none

 integer :: i,j,r,x,y,lj

 call propagatebonds()
 call makeloops()

 do j=1,nn
    x=xy(1,j)
    y=xy(2,j)
    lj=loopnmbr(j)
    do r=0,lx/2
       i=xyi(mod(x+r,lx),y)      
       if (lj==loopnmbr(i)) scor(r)=scor(r)+0.75d0
    enddo
 enddo

 nmsr=nmsr+1

 end subroutine measure
!----------------------!

!-------------------------!
 subroutine writebindata()
!-------------------------!
 use system; use mdata; implicit none

 integer :: r

 scor=scor/(dble(nn)*dble(nmsr))
 open(10,file='cor.dat',status='unknown',position='append')
 do r=0,lx/2
    write(10,'(i6,f18.14)')r,scor(r)
 enddo
 close(10)

 nmsr=0
 scor=0.d0

 end subroutine writebindata
!---------------------------!

!----------------------!
 subroutine writeconf()
!----------------------!
 use system; implicit none
 
 integer :: i

 open(10,file='conf',status='replace')
 do i=1,nn
    write(10,'(2i3,2i7)')lspn(i),rspn(i),lbnd(i),rbnd(i)
 enddo
 do i=0,mm-1
    write(10,*)oper(i) 
 enddo
 close(10)

 end subroutine writeconf
!------------------------!

!---------------------!
 subroutine readconf()
!---------------------!
 use system; implicit none
 
 integer :: i

 open(10,file='conf',status='old')
 do i=1,nn
    read(10,*)lspn(i),rspn(i),lbnd(i),rbnd(i)
 enddo
 do i=0,mm-1
    read(10,*)oper(i)
 enddo
 close(10)

 end subroutine readconf
!-----------------------!

!---------------------------!
 subroutine initsystem(init)
!---------------------------!
 use system; use mdata; implicit none
 
 integer :: i,j,a,x,y,init
 real(8) :: ran

 external :: ran

 mm=mm*nn/2
 allocate(oper(0:mm))
 allocate(vrtx(0:4*mm+4*nn))
 allocate(lspn(nn))
 allocate(rspn(nn))
 allocate(lbnd(nn))
 allocate(rbnd(nn))
 allocate(plbnd(nn))
 allocate(prbnd(nn))
 allocate(frst(nn))
 allocate(last(nn))
 allocate(loopnmbr(nn))

 if (init==0) then
    do i=1,nn
       lspn(i)=(-1)**(xy(1,i)+xy(2,i))   
    enddo
    rspn=lspn
    do y=0,ly-1
    do x=0,lx-1,2
       i=xyi(x,y)
       j=xyi(x+1,y)
       lbnd(i)=j
       lbnd(j)=i
    enddo
    enddo
    rbnd=lbnd
    do i=0,mm-1
       oper(i)=2*(int(ran()*nb)+1)
    enddo
 else
    call readconf()
 endif


 allocate(hamp(0:lx/2,0:ly/2))
 if (ly==1) then
    a=2
 elseif (ly==2) then
    a=4
 else
    a=3
 endif
 do y=0,ly/2
 do x=0,lx/2
    hamp(x,y)=1.d0/sqrt(dble(x**2+y**2))**a
 enddo
 enddo

 end subroutine initsystem
!-------------------------!

!------------------------!
 subroutine makelattice()
!----------------------------------------------------------------------------------------!
! Constructs the lattice, in the form of the list of sites connected by nearest-neighbor
! bonds 'bsite' (defining the hamiltonian) and next-nearest-neighbor sites 'dsite' (uses.
! in the sampling of the valence-bond trial state.
! There are three types of lattices and boundary conditions, depending on ly:
! ly=1  : 1D chain (periodic only in x-direction)
! ly=2  : 2-leg ladder (periodic only in x-direction)
! ly>2  : generic 2d system (periodic in both directions)
! Note that the system is frustrated for ly>2 and odd and the program would not give
! correct results in that case (unless one changes to open boundaries in the y-direction)
!----------------------------------------------------------------------------------------!
 use system; implicit none

 integer :: s,x1,x2,y1,y2

 nn=lx*ly

 allocate(xy(2,nn))
 allocate(xyi(0:lx-1,0:ly-1))

 do s=1,nn
    x1=mod(s-1,lx)
    y1=(s-1)/lx
    xy(1,s)=x1
    xy(2,s)=y1
    xyi(x1,y1)=s
 enddo

 if (ly==1) then        ! 1d chain
    nb=lx
    nd=lx
 elseif (ly==2) then    ! 2-leg ladder
    nb=nn+lx
    nd=nn
 else                   ! general 2d periodic lattice
    nb=2*nn
    nd=2*nn
 endif

 allocate(bsite(2,nb))
 allocate(dsite(2,nd))

 if (ly==1) then        ! 1D chain (periodic only in x-direction)
    do x1=0,lx-1
       s=xyi(x1,0)
       x2=mod(x1+1,lx)
       bsite(1,s)=s
       bsite(2,s)=xyi(x2,0)
       x2=mod(x1+2,lx)
       dsite(1,s)=s
       dsite(2,s)=xyi(x2,0)
    enddo
 elseif (ly==2) then    ! 2-leg ladder chain (periodic only in x-direction)
    do y1=0,ly-1
    do x1=0,lx-1
       s=xyi(x1,y1)
       x2=mod(x1+1,lx)
       y2=y1
       bsite(1,s)=s
       bsite(2,s)=xyi(x2,y2)
       if (y1==0) then
          x2=x1
          bsite(1,s+nn)=s
          bsite(2,s+nn)=xyi(x2,1)       
          x2=mod(x1+1,lx)
          dsite(1,s)=s
          dsite(2,s)=xyi(x2,1)
          x2=mod(x1-1+lx,lx)
          dsite(1,s+lx)=s
          dsite(2,s+lx)=xyi(x2,1)
       endif
    enddo
    enddo
 else                   ! general fully periodic 2D system
    do y1=0,ly-1
    do x1=0,lx-1
       s=xyi(x1,y1)
       x2=mod(x1+1,lx)
       y2=y1
       bsite(1,s)=s
       bsite(2,s)=xyi(x2,y2)
       x2=x1
       y2=mod(y1+1,ly)
       bsite(1,s+nn)=s
       bsite(2,s+nn)=xyi(x2,y2)       
       x2=mod(x1+1,lx)
       y2=mod(y1+1,ly)
       dsite(1,s)=s
       dsite(2,s)=xyi(x2,y2)
       x2=mod(x1-1+lx,lx)
       dsite(1,s+nn)=s
       dsite(2,s+nn)=xyi(x2,y2)
    enddo
    enddo
 endif

 end subroutine makelattice
!--------------------------!

!--------------------------!
 subroutine deallocateall()
!--------------------------! 
 use system; use mdata; implicit none

 deallocate(oper)
 deallocate(lspn)
 deallocate(rspn)
 deallocate(lbnd)
 deallocate(rbnd)
 deallocate(plbnd)
 deallocate(prbnd)
 deallocate(loopnmbr)
 deallocate(frst)
 deallocate(last)
 deallocate(vrtx)
 deallocate(hamp)
 deallocate(bsite)
 deallocate(dsite)
 deallocate(xyi)
 deallocate(xy)
 deallocate(scor)

 end subroutine deallocateall
!----------------------------!

!----------------------!
 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
!----------------------!