program LennardJonesLiquid
       implicit none
       integer L,N,MCS,i,j,k,d, m, sim_time, check_time
       parameter (L=60, MCS=230000, sim_time = 30000, check_time = 100)
       real r_c, x_step, y_step,r_step, ran1, P, dr, ro, PI, T
       parameter (r_c = 2.5, x_step = 0.01,y_step = 0.01,r_step = 0.01,
     &                  dr = 0.1 ,PI = 3.141592653, ro = 0.82, T = 0.7)
       ! 0.8 - liquid & 0.93 - solid
       real xnew, ynew, dx, dy, r_ij2, Rnew2, dU
       real dx_new, dy_new
       real, allocatable:: x(:), y(:)
       d =-1
       N = int( ro*L**2 + 0.0)
       write(*,*) 'N = ', N
       allocate (x(N), y(N))
       call Liquid(L,N,x,y,ro)
       call showmatrix(N,x,y)
        do k=1, MCS
          do i=1, N
             dU = 0
             xnew = x(i) + (ran1(d)-0.5)*x_step
             ynew = y(i) + (ran1(d)-0.5)*y_step
             if(xnew.gt.L) xnew = xnew - L
             if(ynew.gt.L) ynew = ynew - L
             if(xnew.lt.0.) xnew = xnew + L
             if(ynew.lt.0.) ynew = ynew + L
              do j=1, N
               if(j.ne.i) then
                  dx = abs( x(i) - x(j) )
                  dy = abs( y(i) - y(j) )
                  if(dx.gt.L/2.) dx = L - dx
                  if(dy.gt.L/2.) dy = L - dy
                  r_ij2 = dx**2 + dy**2

                  dx_new = abs( xnew-x(j) )
                  dy_new = abs( ynew-y(j) )
                  if(dx_new.gt.L/2.) dx_new = L - dx_new
                  if(dy_new.gt.L/2.) dy_new = L - dy_new
                  Rnew2 = dx_new**2 + dy_new**2
                  if(Rnew2.le.1.0) goto 7
                  !r^2<=1 then r<=1, sqrt() not needed

                  if(r_ij2.le.r_c.and.Rnew2.le.r_c) then
                     dU = dU + 4.0*(1/Rnew2**6 - 1/r_ij2**6)
                  endif
               endif
              enddo
              if(ran1(d).le.min(1.0,exp(-dU/T))) then
                  x(i) = xnew
                  y(i) = ynew
              endif
7           continue
          enddo
          write(*,*) k
        enddo
        call showmatrix(N,x,y)
        write(*,*) 'Done'
        pause

       end program

       subroutine Liquid(L,N,x,y,ro)
        integer L, N, M, k
        real x(N), y(N), d, ro
        x(1) = 0.5
        y(1) = 0.5
        M = (L-1)*(L-1)
        k = int(N/2.)
        if(N.ge.M) then
           d = 1.0
        else
           d = 1.0 + 1./(ro*L)
        endif
        write (*,*) 'd = ',d
        do i=2, k
          y(i) = y(i-1)
          x(i) = x(i-1) + d
          if(x(i).gt.L-0.5) then
             y(i) = y(i) + d
             x(i) = 0.5
          endif
        enddo
        y(k) = L - 0.5
        x(k) = L - 0.5
        do i=k+1, N
          y(i) = y(i-1)
          x(i) = x(i-1) - d
          if(x(i).lt.0.5) then
             y(i) = y(i) - d
             x(i) = L - 0.5
          endif
        enddo
       end

       subroutine showmatrix(n,x,y)
       integer n
       real x(n), y(n)
       open(99,file='matrix-0.1,L=60,T=2.txt')
       do i=1, n
          write(99,*) x(i), y(i)
       enddo
       close(99)
       end


       FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM

      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END