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 3D 
c *****   model of the Lennard-Jones 6-12 fluid.                            
c *****  Here epsilon and sigma are both assumed to be 1
c *****  This code is not optimized for performance, rather to illustrate
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=0.1, beta=1./T)
c **** Maximum displacement
       parameter(dxmax = 0.3  )
c **** Box size
       parameter(bx=15, by=15)
c **** Cut-off radius
       parameter(rcut=3, r2cut=rcut*rcut, rcutdx = rcut+dxmax )
c **** Number of particles
       parameter(Npart=100)
c **** Random number generator's seed
       parameter(iseed=1234567)

c **** Set up arrays 
       dimension x(2,Npart) 

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, 

       Nlat =  int(sqrt(real(Npart)))
       dlatx = bx/Nlat 
       dlaty = by/Nlat 

       ipart = 0
       do ix = 1,Nlat
       do iy = 1,Nlat
        ipart=ipart+1
        x(1,ipart) = ix * dlatx
        x(2,ipart) = iy * dlaty
       enddo
       enddo

c **** Initial Energy 

       Energy = 0.0d0
       do i = 1,Npart-1
       do j = i+1,Npart
        dx = x(1,i)-x(1,j) 
        dy = x(2,i)-x(2,j) 
c **** Here are the Periodic Boundary Conditions
        dx=dx - bx * anint(dx/bx)
        dy=dy - by * anint(dy/by)
        r2ij=dx*dx+dy*dy
c ****  Cut off. Go to the next j if too far away
        if(r2ij.gt.r2cut) goto 99
        r2ij=1.0d0/(r2ij*r2ij*r2ij)
        Energy = Energy -4.* ( r2ij - r2ij*r2ij )
99      continue
       enddo
       enddo

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
c *** Loop through the particles
       do i = 1,Npart
c *** It's random displacement
       dispx = (rand() - 0.5d0) * dxmax 
       dispy = (rand() - 0.5d0) * dxmax 
       xnew = x(1,i) + dispx
       ynew = x(2,i) + dispy

c **** Calculate the energy change 
c **** w = Eijnew - Eijold
       w = 0.0d0
       do  j = 1,Npart  
c **** No interactions with itself       
        if(i.eq.j) goto 90
c **** Old interactions contributions
        dxold = x(1,i)-x(1,j) 
        dyold = x(2,i)-x(2,j) 
c **** Here are the Periodic Boundary Conditions
        dxold = dxold - bx * anint(dxold/bx)
        dyold = dyold - by * anint(dyold/by)
        r2ijold=dxold*dxold+dyold*dyold
       
c **** Save the trouble if really too far away
c **** pick the next j
        if(dsqrt(r2ijold).gt.rcutdx) goto 90

c **** New interactions contributions
        dxnew = xnew  -x(1,j) 
        dynew = ynew  -x(2,j) 
c **** Here are the Periodic Boundary Conditions
        dxnew = dxnew - bx * anint(dxnew/bx)
        dynew = dynew - by * anint(dynew/by)
        r2ijnew=dxnew*dxnew+dynew*dynew
c **** w = w + Enew
c ****  Cut off. No energy if too far away 
        if(r2ijnew.lt.r2cut) then
         r2ijnew=1.0d0/(r2ijnew*r2ijnew*r2ijnew)
          w = w  - 4.* ( r2ijnew - r2ijnew*r2ijnew )
        endif 
c **** w = w - Eold
c ****  Cut off. No energy if too far away 
        if(r2ijold.lt.r2cut) then
         r2ijold=1.0d0/(r2ijold*r2ijold*r2ijold)
          w = w  + 4.* ( r2ijold - r2ijold*r2ijold )
        endif 
c **** end loop over the j's
90      continue
        enddo

c **** do the Metroplis acceptance criterium
       if( ( w.lt.0.d0 ).or.( exp(-beta*w).gt.rand() ) ) then
c **** accept the move, upgrade x_i, total energy, and energy
c **** Bring the particle back into the box if it ran away
              xnew = xnew + bx * int(xnew/bx)
              ynew = ynew + by * int(ynew/by)
              x(1,i)=xnew
              x(2,i)=ynew
              Energy = Energy + w
          endif
c **** end loop over the i's
       enddo

c **** increment averages
         avE = avE + Energy 
         avE2= avE2+ Energy * Energy

c **** print out the running average
         if(mod(icycle,Nprint).eq.1) then
             Erunav = avE / icycle 
             E2runav = avE2 / icycle 
             E2runav = sqrt((E2runav-Erunav*Erunav))
             write(6,2) icycle, Erunav, E2runav
2      format(1x,i8,"cycles Average Energy = ",e12.5,
     >       " Mean.Sqd.Dev", e12.5)
         endif
c **** End loop over the MC cycles
       enddo

c *** That's all folks

        end
