program n_particles_LJ_PBC
      parameter(n=3)
      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
        write(*,*)'What is initial vx of particle',i,'of',n,'?'
        read(*,*) vx(i)
        write(*,*)'What is initial vy of particle',i,'of',n,'?'
        read(*,*) vy(i)
      end do

      x_range = 10.0
      y_range = 10.0

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

      do i=1,n
       x(i) = 5.0*(i-1)/(n-1)
       y(i) = 5.0*(i-1)/(n-1)
      end do

      if( pgopen('/xwin') <= 0 ) stop
      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

      do i_step =1, 10000000

c find forces from pairs
        do i=1,n
          force_x(i)=0.0
          force_y(i)=0.0
        end do
        do i=1,n-1
          do j=i+1,n
c           call lj_force(x(i),y(i),x(j),y(j),force_x_tmp,force_y_tmp)
            call PBC_lj_force(x(i),y(i),x(j),y(j),force_x_tmp,
     &      force_y_tmp,x_range,y_range)
            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
          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)
        call 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
c call force again at new position (to do)
c       call PBC_lj_force(x(i)
c       vx(i) = vx(i) + (force_x(i)+force_x_old(i))/2 * dt
c       vy(i) = vy(i) + (force_y(i)+force_y_old(i))/2 * dt
        vx(i) = vx(i) + force_x(i) * dt
        vy(i) = vy(i) + force_y(i) * dt
c       vx(i) = vx(i) + force_x_old(i) * dt
c       vy(i) = vy(i) + force_y_old(i) * 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
      enddo
      call pgclos
      end


      subroutine PBC_distance(x1,y1,x2,y2,dis,x2mx1,y2my1,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 PBC_lj_force(x1,y1,x2,y2,force_x,force_y,x_range,
     &y_range)
       sigma = 1.0
       epsilon = 1.0
c      r = sqrt( (x2-x1)**2 + (y2-y1)**2 )
       call PBC_distance(x1,y1,x2,y2,dis,x2mx1,y2my1,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
c      force_x = force*(x2-x1)/r
c      force_y = force*(y2-y1)/r
       force_x = force* x2mx1 /r
       force_y = force* y2my1 /r
      end subroutine

      subroutine lj_force(x1,y1,x2,y2,force_x,force_y)
       sigma = 1.0
       epsilon = 1.0
       r = sqrt( (x2-x1)**2 + (y2-y1)**2 )
       s_o_r = sigma / r
       force = (-1)*24*epsilon*(2*s_o_r**12-s_o_r**6)/r
       force_x = force*(x2-x1)/r
       force_y = force*(y2-y1)/r
      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