module UtilIsing
   implicit none
   integer,parameter::L=100             ! linear size of system
   integer,parameter::M=30000            ! no. of spin flips
   real,parameter::tempr=3.5            ! temperature
   integer*1,dimension(0:L-1)::s        ! spin variables s(j)
   real,dimension(-1:1)::W           ! exp(-beta*deltaE)
   integer,dimension(0:L-1)::rnn,lnn
end module


program ising
   use UtilIsing
  implicit none
   call initialize()
   !call snapshot()
   call energy()

end program

subroutine initialize()
   use UtilIsing
   implicit none
   integer::p,j
   call random_seed()
   do p=-1,1            ! tabulaton of botzman weights
      W(p)=exp(-4.0*p/tempr)
   enddo
   ! tabulation of rnn and lnn for right and left n.n. 
   lnn(0)=L-1;  rnn(0)=1;
   lnn(L-1)=L-2;rnn(L-1)=0;        
   do j=1,L-2
      lnn(j)=j-1
      rnn(j)=j+1
   enddo
end subroutine

subroutine metropolis()
   use UtilIsing
   implicit none
   integer::j,p,delE
   real::r
   call random_number(r); j=int(r*L)         ! 0 <= j < L
   delE=2*s(j)*(s(lnn(j))+s(rnn(j)))
   p=delE/4
   call random_number(r)
   if(r<W(p)) s(j)=-s(j)
end subroutine

subroutine snapshot()
   use UtilIsing
   implicit none
   integer::E0, step,j,p,delE,counter
   real::r
   s=1; counter=0
   open(15,file="snapshot.dat")
   do step=1,M
      call random_number(r); j=int(r*L)         ! 0 <= j < L
      delE=2*s(j)*(s(lnn(j))+s(rnn(j)))
      p=delE/4
      call random_number(r)
      if(r<W(p)) then
         counter=counter+1
         s(j)=-s(j)
         do j=0,L-1
            if(s(j)==1) write(15,*) counter,j
         enddo
      endif
   enddo
   close(15)
end subroutine


subroutine energy()
   use UtilIsing
   implicit none
   integer::step,j,p,E,delE,counter
   real::r
   open(16,file="energy.dat")
   s=1;    ! Ferromagnetic start
   ! High temperature start below:
   !do j=0,L-1; call random_number(r);if (r<0.5) s(j)=-1;enddo
   E=0; do j=0,L-1
      E=E-s(j)*s(rnn(j))
   enddo
   counter=0
   do step=1,M
      call random_number(r); j=int(r*L)
      delE=2*s(j)*(s(lnn(j))+s(rnn(j)))
      p=delE/4
      call random_number(r)
      if(r<W(p)) then
         counter=counter+1
         s(j)=-s(j)
         E=E+delE
         write(16,*) counter,real(E)/real(L)
      endif
   enddo
   close(16)
end subroutine
