From 03cd33443240ae62be82fb8f98cf52df1a441182 Mon Sep 17 00:00:00 2001
From: Connor Moore <connor@hhmoore.ca>
Date: Sat, 21 Feb 2026 23:36:50 -0500
Subject: Added neighbour list function to class code.

---
 class-langevin-motion.f90 | 260 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 260 insertions(+)
 create mode 100644 class-langevin-motion.f90

diff --git a/class-langevin-motion.f90 b/class-langevin-motion.f90
new file mode 100644
index 0000000..b52396c
--- /dev/null
+++ b/class-langevin-motion.f90
@@ -0,0 +1,260 @@
+module globals
+! Global variables
+implicit none
+integer, parameter :: n=100
+double precision, parameter :: pi=2q0*asin(1q0) ! numerical constant
+double precision, parameter :: L=1.0
+logical, dimension(n) :: is_tracked = .TRUE.
+end module globals
+
+module sectors
+use globals
+implicit none
+contains
+
+subroutine ONETO2(id, box_size, x, y)
+    implicit none
+    integer, intent(in)  :: id
+    integer, intent(in)  :: box_size
+    integer, intent(out) :: x, y
+
+    y = id / box_size
+    x = mod(id, box_size)
+end subroutine ONETO2
+
+function TWOtoONE(x, y, box_size) result(id)
+    implicit none
+    integer, intent(in) :: x, y, box_size
+    integer :: id
+
+    id = y * box_size + x
+end function TWOtoONE
+
+function get_neighbour_ids(p, N) result (neighbours)
+    integer, intent(in) :: p !> The sequential position
+    integer, intent(in) :: N !> The grid size (i.e., NxN sectors)
+    integer, allocatable, dimension(:) :: neighbours !> The list of neighbours
+    integer :: i,j,k=0
+    integer :: x_cell, y_cell, x_test, y_test
+    integer :: max_list(8)=-1
+
+    call ONETO2(p,N,x_cell,y_cell)
+
+    !> Start by just getting all of the neighbours
+    do i = -1, 1
+        do j = -1,1
+            if(i == 0 .and. j == 0) cycle !> Skip the cell itself
+            
+            !> The "test" coordinate
+            x_test = x_cell + i
+            y_test = y_cell + j
+
+            !> If the coordinates are real (within the expected range), add them to the list:
+            if((x_test .ge. 0 .and. x_test .lt. N) .and. (y_test .ge. 0 .and. y_test .lt. N)) then
+                !> Increment the number of correct coordinates
+                k = k + 1
+                !> Flip a zero to the new coordinate in the max list
+                max_list(k) = TWOtoONE(x_test, y_test, N)
+            end if
+        end do
+    end do
+
+    !> Return a vector of all non-zero elements
+    neighbours = pack(max_list, max_list.ge.0)
+
+end function get_neighbour_ids
+    
+end module sectors
+
+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
-- 
cgit v1.2.3