!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      program main
      implicit none
      integer(8),parameter :: nclasses=64
      integer(8),parameter :: n1x=50
      integer(8),parameter :: n2x=50
      real(8)   ,parameter :: hint=-1.d0  ! spin-spin interaction
      real(8)   ,parameter :: kbt=1.d0
      integer(8),parameter :: nstep=10**6 ! number of mc steps
      integer(8)           :: nmem(nclasses,n2x)
      integer(8)           :: nmem1
      integer(8)           :: state(n1x,n2x)  !occupation
      integer(8)           :: nnarr(n1x,n2x)  !#(nearest neighbors)
      integer(8)           :: class(4,n1x,n2x)
      integer(8)           :: i,istep,icount
      integer(8)           :: nn,nn1,nn2,isvar,idir
      integer(8)           :: i1mm,i1m,i1,i1p,i1pp
      integer(8)           :: i2mm,i2m,i2,i2p,i2pp
      integer(8)           :: j1mm,j1m,j1,j1p,j1pp
      integer(8)           :: j2mm,j2m,j2,j2p,j2pp
      integer(8)           :: iclass,thisclass
      integer(8)           :: thisdir
      real(8)              :: x(nclasses)
      real(8)              :: w(nclasses)
      real(8)              :: p1(nclasses)
      real(8)              :: b
      real(8)              :: deltae,e1,e2
      real(8)              :: svar
      real(8)              :: ran
!     ************************************************************************************
!     ====================================================================================
!     ==  prepare initial state                                                         ==
!     ====================================================================================
      do i1=1,n1x
        do i2=1,n2x
          call random_number(ran)
          state(i1,i2)=nint(ran)
        enddo
      enddo
!
!     ====================================================================================
!     ==  set up neighorlist                                                            ==
!     ====================================================================================
      call neighborlist(n1x,n2x,state,nnarr)
!
!     ====================================================================================
!     ==  set up process table                                                          ==
!     == a process is characterized by the number of neighbors at the initial and the   ==
!     == final state                                                                    ==
!     ====================================================================================
!     ==  class= 1:   nn1=0, nn2=1                                                      ==
!     ==  class= 2:   nn1=1, nn2=1                                                      ==
!     ==  class= 3:   nn1=2, nn2=1                                                      ==
!     ==  class= 4:   nn1=3, nn2=1                                                      ==
!     ==  class= 5:   nn1=0, nn2=2                                                      ==
!     ==  class= 6:   nn1=1, nn2=2                                                      ==
!     ==    ...                                                                         ==
!     ==  class=15:   nn1=2, nn2=4                                                      ==
!     ==  class=16:   nn1=3, nn2=4                                                      ==
!     ====================================================================================
      i=0
      do nn2=1,4
        do nn1=0,3
          i=i+1
          e1=hint*real(nn1)
          e2=hint*real(nn2-1)
          deltae=e2-e1
          p1(i)=exp(-deltae/kbt)       ! probability for a spin flip
        enddo
      enddo
      svar=1.d0/sum(p1)
      p1(:)=p1(:)*svar
!
!     == report probabilities ===========================================================
      print*,'====== probability per class ======='
      do i1=1,4
        write(*,fmt='(4f20.10)')p1(1+4*(i1-1):4*i1)
      enddo
!
!     ====================================================================================
!     ==  determine process classes                                                     ==
!     ==  class(1,i1,i2) process type for jump from (i1,i2) to the right (i2->i2+1)     ==
!     ==  class(2,i1,i2) process type for jump from (i1,i2) to the left  (i2->i2-1)     ==
!     ==  class(3,i1,i2) process type for jump from (i1,i2) up           (i1->i1+1)     ==
!     ==  class(4,i1,i2) process type for jump from (i1,i2) down         (i1->i1-1)     ==
!     ==  count the number of processes in each class                                   ==
!     ====================================================================================
      call processlist(n1x,n2x,state,class,nmem)
!
!     ====================================================================================
!     ==  monte carlo loop                                                              ==
!     ====================================================================================
      do istep=1,nstep
        if(modulo(istep,nstep/100_8).eq.0)print*,'percent finished ',100*istep/nstep
!
!       ==================================================================================
!       == select a class                                                              ==
!       ==================================================================================
        do i=1,nclasses
          w(i)=real(sum(nmem(i,:)))
        enddo
        w(:)=w(:)/sum(w(:))
        x(:)=p1(:)*w(:)
        x(:)=x(:)/sum(x(:))
        call random_number(ran)
        svar=0.d0
        do i=1,nclasses
          svar=svar+x(i)
          if(svar.gt.ran) then
            thisclass=i
            exit
          end if
        enddo
!
!       ==================================================================================
!       == select a random site from this class                                         ==
!       ==================================================================================
        nmem1=sum(nmem(thisclass,:))
        call random_number(ran)
        isvar=1+int(ran*nmem1)
        isvar=min(isvar,nmem1)   ! atom number in this class
!
        icount=0
        do i2=1,n2x
          if(icount+nmem(thisclass,i2).lt.isvar) then
            icount=icount+nmem(thisclass,i2)
          else
            do i1=1,n1x
              do idir=1,4
                if(class(idir,i1,i2).eq.thisclass) then
                  icount=icount+1
                  if(icount.eq.isvar) then
                    thisdir=idir
                    goto 1000
                  end if
                end if
              enddo
            enddo
          end if
        enddo
1000    continue
!
!       ==================================================================================
!       == move atom
!       ==================================================================================
!       == neighborhood of initial state        
        i1pp=1+modulo(i1+1,n1x)
        i1p =1+modulo(i1,n1x)
        i1m =1+modulo(i1-2,n1x)
        i1mm=1+modulo(i1-3,n1x)
        i2pp=1+modulo(i2+1,n2x)
        i2p =1+modulo(i2,n2x)
        i2m =1+modulo(i2-2,n2x)
        i2mm=1+modulo(i2-3,n2x)
!
!       == final state
        j1=i1
        j2=i2
        if(thisdir.eq.1) then
          j2=i2p
        else if(thisdir.eq.2) then
          j2=i2m
        else if(thisdir.eq.3) then
          j1=i1p
        else if(thisdir.eq.4) then
          j1=i1m
        end if
!
!       == neighborhood of final state==================================================
        j1pp=1+modulo(j1+1,n1x)
        j1p =1+modulo(j1,n1x)
        j1m =1+modulo(j1-2,n1x)
        j1mm=1+modulo(j1-3,n1x)
        j2pp=1+modulo(j2+1,n2x)
        j2p =1+modulo(j2,n2x)
        j2m =1+modulo(j2-2,n2x)
        j2mm=1+modulo(j2-3,n2x)
!
!       ==  consistency check ============================================================
        if(state(i1,i2).eq.0) then
          print*,'class ',class(thisdir,i1,i2)
          stop 'error initial state unoccupied'
        end if
        if(state(j1,j2).eq.1) then
          print*,'class ',class(thisdir,i1,i2)
          print*,'x1 ',i1,i2
          print*,'x2 ',j1,j2
         stop 'error final state occupied'
        end if

!       == perform jump ==================================================================
        state(i1,i2)=0
        state(j1,j2)=1
!
!       == update neighborlist ===========================================================
        nnarr(i1m,i2)=nnarr(i1m,i2)-1
        nnarr(i1p,i2)=nnarr(i1p,i2)-1
        nnarr(i1,i2p)=nnarr(i1,i2p)-1
        nnarr(i1,i2m)=nnarr(i1,i2m)-1
        nnarr(j1m,j2)=nnarr(j1m,j2)+1
        nnarr(j1p,j2)=nnarr(j1p,j2)+1
        nnarr(j1,j2p)=nnarr(j1,j2p)+1
        nnarr(j1,j2m)=nnarr(j1,j2m)+1
!
!       ==================================================================================
!       == update process list                                                          ==        
!       ==================================================================================
        do idir=1,4
          call updateprocesslist(idir,i1,i2,n1x,n2x,state,nnarr,class,nmem)
          call updateprocesslist(idir,j1,j2,n1x,n2x,state,nnarr,class,nmem)
        enddo
!
!=========================================================================
!       == two steps straight right,left,up,down
        call updateprocesslist(1_8,i1,i2mm,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(2_8,i1,i2pp,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,i1mm,i2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,i1pp,i2,n1x,n2x,state,nnarr,class,nmem)

        call updateprocesslist(2_8,i1,i2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(1_8,i1,i2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,i1m,i2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,i1p,i2,n1x,n2x,state,nnarr,class,nmem)

!       ==  left lower diagonal
        call updateprocesslist(1_8,i1m,i2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,i1m,i2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(2_8,i1m,i2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,i1,i2m,n1x,n2x,state,nnarr,class,nmem)

!       == right lower diagonal
        call updateprocesslist(2_8,i1m,i2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,i1m,i2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(1_8,i1m,i2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,i1,i2p,n1x,n2x,state,nnarr,class,nmem)

!       == left upper diagonal
        call updateprocesslist(1_8,i1p,i2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,i1p,i2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(2_8,i1p,i2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,i1,i2m,n1x,n2x,state,nnarr,class,nmem)
!
!       == right upper diagonal
        call updateprocesslist(2_8,i1p,i2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,i1p,i2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(1_8,i1p,i2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,i1,i2p,n1x,n2x,state,nnarr,class,nmem)
!
!=========================================================================
!       == two steps straight right,left,up,down
        call updateprocesslist(1_8,j1,j2mm,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(2_8,j1,j2pp,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,j1mm,j2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,j1pp,j2,n1x,n2x,state,nnarr,class,nmem)

        call updateprocesslist(2_8,j1,j2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(1_8,j1,j2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,j1m,j2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,j1p,j2,n1x,n2x,state,nnarr,class,nmem)

!       ==  left lower diagonal
        call updateprocesslist(1_8,j1m,j2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,j1m,j2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(2_8,j1m,j2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,j1,j2m,n1x,n2x,state,nnarr,class,nmem)

!       == right lower diagonal
        call updateprocesslist(2_8,j1m,j2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,j1m,j2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(1_8,j1m,j2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,j1,j2p,n1x,n2x,state,nnarr,class,nmem)

!       == left upper diagonal
        call updateprocesslist(1_8,j1p,j2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,j1p,j2m,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(2_8,j1p,j2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,j1,j2m,n1x,n2x,state,nnarr,class,nmem)
!
!       == right upper diagonal
        call updateprocesslist(2_8,j1p,j2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(4_8,j1p,j2p,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(1_8,j1p,j2,n1x,n2x,state,nnarr,class,nmem)
        call updateprocesslist(3_8,j1,j2p,n1x,n2x,state,nnarr,class,nmem)

      enddo
!
!     ====================================================================================
!     ==  analyze result                                                                ==
!     ====================================================================================
      open(11,file='diffusion.out')
      call plotstate(11,n1x,n2x,state)
      close(11) 
      stop
      end
!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      subroutine neighborlist(n1x,n2x,state,nnarr)
      implicit none
      integer(8),intent(in) :: n1x
      integer(8),intent(in) :: n2x
      integer(8),intent(in) :: state(n1x,n2x)
      integer(8),intent(out):: nnarr(n1x,n2x)
      integer(8)            :: i1,i2
      integer(8)            :: i1p,i1m,i2p,i2m
!     ************************************************************************************
      do i1=1,n1x
        do i2=1,n2x
          i1p=1+modulo(i1,n1x)
          i1m=1+modulo(i1-2,n1x)
          i2p=1+modulo(i2,n2x)
          i2m=1+modulo(i2-2,n2x)
          nnarr(i1,i2)=state(i1p,i2)+state(i1m,i2)+state(i1,i2p)+state(i1,i2m)
        enddo
      enddo          
      return
      end
!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      subroutine updateprocesslist(idir,i1,i2,n1x,n2x,state,nnarr,class,nmem)
!     ====================================================================================
!     ==  update process classes                                                        ==
!     **  assumes that occupations "state" and neighborlist "nnarr" are updated         **
!     ==  class(1,i1,i2) process type for jump from (i1,i2) to the right (i2->i2+1)     ==
!     ==  class(2,i1,i2) process type for jump from (i1,i2) to the left  (i2->i2-1)     ==
!     ==  class(3,i1,i2) process type for jump from (i1,i2) up           (i1->i1+1)     ==
!     ==  class(4,i1,i2) process type for jump from (i1,i2) down         (i1->i1-1)     ==
!     ==  count the number of processes in each class                                   ==
!     ====================================================================================
      implicit none
      integer(8),intent(in) :: idir
      integer(8),intent(in) :: i1
      integer(8),intent(in) :: i2
      integer(8),intent(in) :: n1x
      integer(8),intent(in) :: n2x
      integer(8),intent(in) :: state(n1x,n2x)
      integer(8),intent(in):: nnarr(n1x,n2x)
      integer(8),intent(inout):: class(4,n1x,n2x)
      integer(8),intent(inout):: nmem(64,n2x)
      integer(8)            :: j1,j2
      integer(8)            :: nn1,nn2
!     ************************************************************************************
!
!     == disconnect all processes ========================================================
      if(class(idir,i1,i2).ne.0) then
        nmem(class(idir,i1,i2),i2)=nmem(class(idir,i1,i2),i2)-1
        class(idir,i1,i2)=0
      end if
!
!     == no hops from this site if not occupied
      if(state(i1,i2).eq.0) return
!
!     == select final position
      j1=i1
      j2=i2
      if(idir.eq.1) then
        j2=1+modulo(i2,n2x)
      else if(idir.eq.2) then
        j2=1+modulo(i2-2,n2x)
      else if(idir.eq.3) then
        j1=1+modulo(i1,n1x)
      else
        j1=1+modulo(i1-2,n1x)
      end if

!     == include hop
      if(state(j1,j2).eq.0) then
        nn1=nnarr(i1,i2)
        nn2=nnarr(j1,j2)
        class(idir,i1,i2)=1+nn1+4*(nn2-1)
        nmem(class(idir,i1,i2),i2)=nmem(class(idir,i1,i2),i2)+1
      end if
      return
      end
!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      subroutine processlist(n1x,n2x,state,class,nmem)
!     ====================================================================================
!     ==  determine process classes                                                     ==
!     ==  class(1,i1,i2) process type for jump from (i1,i2) to the right (i2->i2+1)     ==
!     ==  class(2,i1,i2) process type for jump from (i1,i2) to the left  (i2->i2-1)     ==
!     ==  class(3,i1,i2) process type for jump from (i1,i2) up           (i1->i1+1)     ==
!     ==  class(4,i1,i2) process type for jump from (i1,i2) down         (i1->i1-1)     ==
!     ==  count the number of processes in each class                                   ==
!     ====================================================================================
      implicit none
      integer(8),intent(in) :: n1x
      integer(8),intent(in) :: n2x
      integer(8),intent(in) :: state(n1x,n2x)
      integer(8),intent(out):: class(4,n1x,n2x)
      integer(8),intent(out):: nmem(64,n2x)
      integer(8)            :: i1,i2,idir
      integer(8)            :: i1p,i1m,i2p,i2m
      integer(8)            :: nn1,nn2
      integer(8)            :: nnarr(n1x,n2x)
!     ************************************************************************************
      call neighborlist(n1x,n2x,state,nnarr)
      nmem(:,:)=0
      class(:,:,:)=0
      do i1=1,n1x
        do i2=1,n2x
          if(state(i1,i2).eq.0) then
            class(:,i1,i2)=0
            cycle
          end if
          i1p=1+modulo(i1,n1x)
          i1m=1+modulo(i1-2,n1x)
          i2p=1+modulo(i2,n2x)
          i2m=1+modulo(i2-2,n2x)
          nn1=nnarr(i1,i2)
!         == hop to the right side
          idir=1
          if(state(i1,i2p).eq.1) then
            class(idir,i1,i2)=0  ! no jumps from this site possible
          else
            nn2=nnarr(i1,i2p)
            class(idir,i1,i2)=1+nn1+4*(nn2-1)
            nmem(class(idir,i1,i2),i2)=nmem(class(idir,i1,i2),i2)+1
          end if
!         == hop to the left
          idir=2
          if(state(i1,i2m).eq.1) then
            class(idir,i1,i2)=0
          else
            nn2=nnarr(i1,i2m)
            class(idir,i1,i2)=1+nn1+4*(nn2-1)
            nmem(class(idir,i1,i2),i2)=nmem(class(idir,i1,i2),i2)+1
          end if
!         == hop up
          idir=3
          if(state(i1p,i2).eq.1) then
            class(idir,i1,i2)=0
          else
            nn2=nnarr(i1p,i2)
            class(idir,i1,i2)=1+nn1+4*(nn2-1)
            nmem(class(idir,i1,i2),i2)=nmem(class(idir,i1,i2),i2)+1
          end if
!         == hop down
          idir=4
          if(state(i1m,i2).eq.1) then
            class(idir,i1,i2)=0
          else
            nn2=nnarr(i1m,i2)
            class(idir,i1,i2)=1+nn1+4*(nn2-1)
            nmem(class(idir,i1,i2),i2)=nmem(class(idir,i1,i2),i2)+1
          end if
        enddo
      enddo
      return
      end
!
!     ...1.........2.........3.........4.........5.........6.........6.........7.........8
      subroutine plotstate(nfil,n1x,n2x,state)
!     ************************************************************************************
      implicit none
      integer  ,intent(in) :: nfil
      integer(8),intent(in) :: n1x
      integer(8),intent(in) :: n2x
      integer(8),intent(in) :: state(n1x,n2x)
      integer(8)            :: i1,i2
      character(n2x)        :: string
!     ************************************************************************************
      do i2=1,n2x
        string(i2:i2)='-'
      enddo
      write(nfil,*)'|'//string//'|'
      do i1=1,n1x
        string=''
        do i2=1,n2x
          if(state(i1,i2).eq.1)string(i2:i2)='o'
        enddo
        write(nfil,*)'|'//string//'|'
      enddo
      do i2=1,n2x
        string(i2:i2)='-'
      enddo
      write(nfil,*)'|'//string//'|'
      return
      end
