c ***************************************************************************
c Copyright(C) M.Kotelyanskii, November, 1997
c This code is just an example, it can be freely distributed and modified
c as soon as this notice remains in the file.
c It has been tested on the IBM Risc6000 under AIX. Other systems may require
c changes, particularly the functions calls for the random numbers generation.
c Even though it works properly to the best of the author's knowledge,
c he does not assume any responsibility for any damage caused by it
c ***************************************************************************
c ***************************************************************************
c *****  This is example of the Metropolis Monte Carlo method for the 2D 
c *****  Ising model with coupling constant J=-1 (ferromagnetic)       
c *****  This code is not optimized for performance. It just illustrates
c *****  the algorithm.
c ***************************************************************************

       implicit real*8(a-h,o-z)

c **** Specify parameters:
c **** Number of trials
       parameter(Ncycle=10000, Nprint=Ncycle/10)
c **** Temperature
       parameter(T=3.01, beta=1./T)
c **** Lattice size
       parameter(Nx=4, Ny=4)
c **** Random number generator's seed
       parameter(iseed=1234567)

c **** Set up Latiice
       dimension s(Nx,Ny)

c **** set counters
c **** Running Averages    
c **** Energy
       avE  = 0.0d0
       avE2 = 0.0d0
c **** Magnetization 
       avM  = 0.0d0
       avM2 = 0.0d0

c **** Set up the initial lattice configuration, 

       do ix = 1,Nx
       do iy = 1,Ny
        s(ix,iy) = 1.0d0
       enddo
       enddo

c **** Initial Energy and magnetization

       Energy = 0.0d0
       Am     = 0.0d0
       do ix = 1,Nx
       do iy = 1,Ny
        sn=0.0d0
c **** Here to be changed for the Periodic Boundary Conditions
        if(ix.gt.1 ) sn=sn+s(ix-1,iy)
        if(ix.lt.Nx) sn=sn+s(ix+1,iy)
        if(iy.gt.1 ) sn=sn+s(ix,iy-1)
        if(iy.lt.Ny) sn=sn+s(ix,iy+1)
        Energy=Energy - sn * s(ix,iy) 
        Am = Am - s(ix,iy)
       enddo
       enddo

c **** We've counted each pair twice, correct for it
       Energy=Energy/2.0

       write(6,1) beta, Energy,Am
1     format(" beta= ",e12.5," Initial Energy= ",e12.5," Mag-n= ",e12.5)

c *** set seed for the random numbers generator (notice, that the name for srand
c *** may be different on different architectures )
       call srand(iseed)

c *** Start MC cycle (assuming rand() returns random real*1 from 0 to 1)

       do icycle = 1, Ncycle
       do npicks = 1,Nx*Ny
c *** Pick the node to flip at random
         ix = int(rand()*Nx)+1
         iy = int(rand()*Ny)+1
c *** Calculate the energy change if flipped
         sn=0.0d0
         if(ix.gt.1 ) sn=sn+s(ix-1,iy)
         if(ix.lt.Nx) sn=sn+s(ix+1,iy)
         if(iy.gt.1 ) sn=sn+s(ix,iy-1)
         if(iy.lt.Ny) sn=sn+s(ix,iy+1)
         w =  2.0 * sn * s(ix,iy) 
c **** do the Metroplis acceptance criterium
          if(( w.lt.0.d0 ).or.( exp(-beta*w).gt.rand() )) then
c **** accept the move, upgrade s_i, total magnetization, and energy
              s(ix,iy)=-s(ix,iy) 
              Energy = Energy + w
              Am = Am - 2.0 * s(ix,iy)
          endif
c **** end loop over the lattice
       enddo


c **** increment averages
         avM = avM + Am
         avM2= avM2+ Am * Am
         avE = avE + Energy 
         avE2= avE2+ Energy * Energy

c **** print out the running average
         if(mod(icycle,Nprint).eq.1) then
             Erunav = avE / icycle /Nx/Ny
             E2runav = avE2 / icycle /Nx/Ny/Nx/Ny
             E2runav = sqrt((E2runav-Erunav*Erunav))
             runavM = avM / icycle /Nx/Ny
             runavM2 = avM2 / icycle /Nx/Ny/Nx/Ny
             runavM2 = sqrt((runavM2-runavM*runavM))
        write(6,2) icycle, Energy,Am, Erunav, E2runav, runavM, runavM2
2      format(1x,i8,"cycle Current Total Energy ",e12.5," Mag-n ",e12.5,
     >"  Average Energy = ",e12.5, " Mean.Sqd.Dev ", e12.5,
     >" Average Mag-n = ",e12.5," Mean.Sqd.Dev ", e12.5 )
           endif

c **** End MC cycle
       enddo


c *** That's all folks

        end
