program n_particles_LJ_PBC_rand_TP
      parameter(n=8)
      real x(n),y(n),x_old(n),y_old(n),vx(n),vy(n)
      real force_x(n),force_y(n),force_x_old(n),force_y_old(n)
      integer pgopen

      do i=1,n
        call random_number(r)
        vx(i)=2*(r-0.5)
        write(*,*)'The initial vx of particle',i,'is',vx(i)
        call random_number(r)
        vy(i)=2*(r-0.5)
        write(*,*)'The initial vy of particle',i,'is',vy(i)
c       write(*,*)'What is initial vx of particle',i,'of',n,'?'
c       read(*,*) vx(i)
c       write(*,*)'What is initial vy of particle',i,'of',n,'?'
c       read(*,*) vy(i)
      end do

      x_range = 10.0
      y_range = 10.0

C In this code, n must be >= 2.

C Diagonal line-up type inital position
      do i=1,n
       x(i) = 8.0*(i-1)/(n-1)
       y(i) = 8.0*(i-1)/(n-1)
       write(*,*)'position of particle',i,'is (',x(i),y(i),')'
      end do

C Cheker board line-up type initial position
      do i=1,n
       x(i) = (x_range/3)*mod(i,3)
       y(i) = (y_range/3)*(i/3)
      end do


      dis_limit = 1.0
c     goto 200
      call random_coord(x(1),y(1),x_range,y_range)
      i_count=1
      do while(i_count<n)
c       write(*,*)'i_count=',i_count
 111      call random_coord(x_temp,y_temp,x_range,y_range)
c         dis_min = dis_limit
          if(i_count>=2)then
            call PBC_distance(x(i_count),y(i_count),x(i_count-1),
     &      y(i_count-1),x2mx1,y2my1,dis,x_range,y_range)
            dis_min = dis
          end if
          do j=1,i_count
            call PBC_distance(x(j),y(j),x_temp,y_temp,x2mx1,y2my1,dis,
     &      x_range,y_range)
            if(i_count==1) dis_min = dis
            write(*,*)'x and y and dis=',x_temp,y_temp,dis
            if(dis<dis_min)dis_min=dis
          end do
          write(*,*)'dis_min=',dis_min,'i_count=',i_count
c         if(dis_min<=dis_limit) goto 111
          if(dis_min<=dis_limit)then
           write(*,*)'Close-contact distance',dis_min,'too short'
           write(*,*)'accept(1) or retry(2)?'
           read(*,*) ikey
           if(ikey==2)then
             dis_min = dis_limit
             goto 111
           end if
          end if 
          i_count=i_count+1
          x(i_count)=x_temp
          y(i_count)=y_temp
          write(*,*)'The i_count is',i_count 
      end do
 200  continue


      if( pgopen('/xwin') <= 0 ) stop
      call pgpap(4.0,1.0)
      call pgenv(-0.15,x_range+0.15,-0.15,y_range+0.15,1,-1)
      call pgsfs(2)

      dt = 0.0001
      dt2 = dt*dt

      write(*,*) 'Before master loop'

      Tcum=0.0
      Pcum=0.0

      do 1000 i_step =1, 1000000

c find forces and work(virial) from pairs
        do i=1,n
          force_x(i)=0.0
          force_y(i)=0.0
        end do
        virial=0.0
        U=0.0
        do i=1,n-1
          do j=i+1,n
            call lj_force(x(i),y(i),x(j),y(j),force_x_tmp,force_y_tmp,
     &      x_range,y_range,x2mx1,y2my1,pot)
            force_x(i)=force_x(i)+force_x_tmp
            force_y(i)=force_y(i)+force_y_tmp
            force_x(j)=force_x(j)-force_x_tmp
            force_y(j)=force_y(j)-force_y_tmp
            virial = virial + x2mx1*force_x_tmp + y2my1*force_y_tmp
            U=U+pot
          end do
        end do
c 

       do i=1,n

c       call wall_reflect(x(i),y(i),vx(i),vy(i),x_range,y_range)
c       call PBC(x(i),y(i),x_range,y_range)
        call strong_PBC(x(i),y(i),x_range,y_range)

        x(i) = x(i) + vx(i) * dt + 0.5 * force_x(i) * dt2
        y(i) = y(i) + vy(i) * dt + 0.5 * force_y(i) * dt2
        vx(i) = vx(i) + (force_x(i)+force_x_old(i))/2 * dt
        vy(i) = vy(i) + (force_y(i)+force_y_old(i))/2 * dt

        call pgbbuf
        call pgsci(0)
        call pgcirc(x_old(i),y_old(i),0.1)
        call pgsci(1)
        call pgcirc(x(i),y(i),0.1)
        call pgebuf

        x_old(i) = x(i)
        y_old(i) = y(i)
        force_x_old(i) = force_x(i)
        force_y_old(i) = force_y(i)
       enddo

c Temperature and Pressure
       ek_cum = 0.0
       do i=1,n
       ek = ek + 0.5*(vx(i)*vx(i) + vy(i)*vy(i))
       end do
c      Ek=(d/2)nkT => T=Ek/(nk) (d=2, unit 1/k)
c      ek=ek/n
       area=x_range*y_range
c      P=(n*k*T + 0.5*virial)/area 
       ek_cum = ek_cum + ek
       V_cum = V_cum + Virial
       T = ek_cum/i_step
       P = T + 0.5*V_cum/i_step
       P = P/area
       if ( mod(i_step,1000) == 0 ) write(*,*)'T=',T,'P=',P

 1000 continue
      call pgclos
      end


      subroutine lj_force(x1,y1,x2,y2,force_x,force_y,x_range,y_range,
     &x2mx1,y2my1,pot)
       sigma = 1.0
       epsilon = 1.0
c      r = sqrt( (x2-x1)**2 + (y2-y1)**2 )
       call PBC_distance(x1,y1,x2,y2,x2mx1,y2my1,dis,x_range,y_range)
       r=dis
       s_o_r = sigma /r
       force = (-1)*24*epsilon*(2*s_o_r**12-s_o_r**6)/r
       force_x = force*(x2mx1)/r
       force_y = force*(y2my1)/r
       pot=4*epsilon*(s_o_r**12-s_o_r**6)
      end subroutine

      subroutine PBC_distance(x1,y1,x2,y2,x2mx1,y2my1,dis,x_range,
     &y_range)
      dis_min=sqrt( (0.5*x_range)**2 + (0.5*y_range)**2 )
       do i=-1,1
         do j=-1,1
           x2mx1_temp=i*x_range+x2-x1
           y2my1_temp=j*y_range+y2-y1
           r=sqrt( (x2mx1_temp)**2 + (y2my1_temp)**2 )
           if(r < dis_min)then
             dis_min = r
             x2mx1=x2mx1_temp
             y2my1=y2my1_temp
             dis=r
           endif
         end do
       end do
      end subroutine

      subroutine wall_reflect(x,y,vx,vy,x_range,y_range)
        if(x>x_range) vx = (-1)*vx
        if(x<0.0) vx = (-1)*vx
        if(y>y_range) vy = (-1)*vy
        if(y<0.0) vy = (-1)*vy
      end subroutine


      subroutine PBC(x,y,x_range,y_range)
        if(x>x_range) x = x - x_range
        if(x<0.0)     x = x + x_range
        if(y>y_range) y = y - y_range
        if(y<0.0)     y = y + y_range
      end subroutine

      subroutine random_coord(x1,y1,x_range,y_range)
        call random_number(temp)
c       x1=x_range*(temp-0.5)
        x1=x_range*(temp)
        call random_number(temp)
c       y1=y_range*(temp-0.5)
        y1=y_range*(temp)
      end subroutine

      subroutine strong_PBC(x,y,x_range,y_range)
 201    if(x>x_range) x = x - x_range
        if(x>x_range) goto 201
 202    if(x<0.0)     x = x + x_range
        if(x<0.0)     goto 202
 203    if(y>y_range) y = y - y_range
        if(y>y_range) goto 203
 204    if(y<0.0)     y = y + y_range
        if(y<0.0)     goto 204
      end subroutine