module langevin
! Initialization and update rule for Langevin particles
use globals
implicit none
double precision :: dt,kT,g,m,eps,sigma,rc                   ! time step size and physical parameters
double precision :: pref1,pref2                 ! auxiliary parameters
double precision, allocatable, dimension(:) :: x,y,vx,vy,ax,ay,vhx,vhy,x0,y0 ! particle positions, accellerations, velocities, half-step velocities, initial positions
contains
subroutine set_parameters

! Set time step and physical parameters
dt=0.0005d0 ! time step size
kT=10d0    ! energy
g=1d0     ! draf coefficient
m=1d0     ! mass of the particles, can be normalized to 1.
eps=1d0
sigma=1d-3
rc=sigma*2d0**(1d0/6d0)
!print *, rc

! Set auxiliary parameters
pref1=g
pref2=sqrt(24d0*kT*g/dt)

end subroutine set_parameters
subroutine initialize_particles
integer :: i
double precision :: ran1,ran2,gr1,gr2
! Give particles initial position and velocity
do i=1,n
   call random_number(ran1)                       ! uses the built-in PRNG, easy but not very accurate
   call random_number(ran2)
   x(i)=L*(ran1-0.5d0)
   x0(i)=x(i)
   y(i)=L*(ran2-0.5d0)
   y0(i)=y(i)
   ax(i)=0d0
   ay(i)=0d0
   call random_number(ran1)
   call random_number(ran2)
   gr1=sqrt(kT/(m))*sqrt(-2*log(ran1))*cos(2*pi*ran2) ! Box-Mueller transform
   gr2=sqrt(kT/(m))*sqrt(-2*log(ran1))*sin(2*pi*ran2)
   vx(i)=gr1
   vy(i)=gr2
end do

end subroutine initialize_particles
end module langevin

module BC
  ! Subroutines related to the boundary conditions
  use globals
  use langevin
  implicit none
contains
  subroutine impose_BC(i)
    integer :: i
    !> Recall we are inside an LxL box centered at the origin

    !> Top boundary
    if(y(i) .GT. L/2) then
        y(i) = L - y(i)
        vhy(i) = -vhy(i)
    end if

    !> Left boundary
    if(x(i) .LT. -L/2) then
        x(i) = -L - x(i)
        vhx(i) = -vhx(i)
    end if

    !> Bottom boundary
    if(y(i) .LT. -L/2) then
        y(i) = -L - y(i)
        vhy(i) = -vhy(i)
    end if

    !> Right boundary
    if(x(i) .GT. L/2) then
        x(i) = L - x(i)
        vhx(i) = -vhx(i)
    end if

    !> Final check
    if(abs(x(i)).GT.L/2 .OR. abs(y(i)).GT.L/2) then
        !> Particle is still outside, don't track it
        is_tracked(i) = .FALSE.
    endif

  end subroutine impose_BC
end module BC

program main
use globals
use langevin
use BC
implicit none
integer :: i,j
double precision :: t,t_max,m1,m2,dij,rx,ry,F
double precision :: wtime,begin,end
double precision, dimension(n) :: ran1, ran2

! Open files
!open(12,file='means')
open(23,file='particle_A')
open(24,file='particle_B')
open(25,file='interactions')

! Allocate arrays
allocate(x(n),y(n),vx(n),vy(n),ax(n),ay(n),vhx(n),vhy(n),x0(n),y0(n))

t=0d0
t_max=10d0     ! integration time

call set_parameters
call initialize_particles

call cpu_time(begin)

!> 1. Update half-velocity
!> 2. Update position
!> 3. Compute accelerations and forces
!> 4. Update all velocities

do while(t.lt.t_max)
    vhx = vx+0.5d0*ax*dt
    vhy = vy+0.5d0*ay*dt

    x = x + vhx*dt
    y = y + vhy*dt

    do i = 1,n
        call impose_BC(i)
    end do

    call random_number(ran1)
    ran1 = ran1-0.5d0
    call random_number(ran2)
    ran2 = ran2-0.50

    ax = 0d0
    ay = 0d0

    !ax = ax - pref1*vhx + pref2*ran1
    !ay = ay - pref1*vhy + pref2*ran2
    ax = ax - pref1*vhx + 0.d0*ran1
    ay = ay - pref1*vhy + 0.d0*ran2

    do i=1,n
      do j=1,n
        if(j.ne.i) then
          if(.not.is_tracked(j)) cycle
          rx=x(j)-x(i)
          ry=y(j)-y(i)
          dij = sqrt(rx**2 + ry**2)
          if(dij.lt.rc) then
              F = 4d0*eps*( -12d0*sigma**12/dij**13 + 6d0*sigma**6/dij**7)
              ax(i)=ax(i)+F*rx/(dij*m)
              ay(i)=ay(i)+F*ry/(dij*m)
          endif
        end if
      end do
    end do
      
    vx=vhx+0.5d0*ax*dt
    vy=vhy+0.5d0*ay*dt

    t = t + dt

    do i=1,n
        print '(e24.12,e24.12)', x(i), y(i)
    end do

    print *
    print *

end do

call cpu_time(end)
!print *,'Wtime=',end-begin
!print *, count(is_tracked), "particles conserved out of ", size(is_tracked)

! De-allocate arrays
deallocate(x,y,vx,vy,ax,ay,x0,y0)
! Close files
!close(11)
!close(12)
close(23)
close(24)
close(25)

end program main