Updated class code

13 files changed, 477 insertions, 277 deletions

diff --git a/class/Makefile b/class/Makefile
new file mode 100644
index 0000000..17979f9
--- /dev/null
+++ b/class/Makefile
@@ -0,0 +1,19 @@
+FC 		= gfortran
+FFLAGS 	= -O4 -fopenmp -fcheck=all
+NPROC 	= $(shell nproc)
+TARGET 	= langdyn.x
+
+all: clean $(TARGET) run
+
+$(TARGET): langdyn.f90
+	$(FC) $(FFLAGS) -o $(TARGET) langdyn.f90
+
+clean:
+	rm -f *.mod *.out *.gif $(TARGET)
+
+run: $(TARGET)
+	OMP_NUM_THREADS=$(NPROC) ./$(TARGET)
+
+plot: trajectories.out kinetic_energies.out
+	gnuplot movie_plot.gnu
+	gnuplot ke_plot.gnu
diff --git a/class/class-langevin-motion.f90 b/class/class-langevin-motion.f90
deleted file mode 100644
index fcdc77a..0000000
--- a/class/class-langevin-motion.f90
+++ /dev/null
@@ -1,192 +0,0 @@
-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
diff --git a/class/ke_compare.gnu b/class/ke_compare.gnu
new file mode 100644
index 0000000..e09d021
--- /dev/null
+++ b/class/ke_compare.gnu
@@ -0,0 +1,15 @@
+set term cairolatex pdf size 5in,3in
+set output "kinetic_energy.tex"
+#set term x11 persist
+set linetype 1 lc rgb '#E41A1C' pt 7 ps 0.3
+set linetype 2 lc rgb '#377EB8' pt 7 ps 0.3
+set linetype 3 lc rgb '#4DAF4A' pt 7 ps 0.3
+set linetype 4 lc rgb '#984EA3' pt 7 ps 0.3
+
+set xlabel "Time [s]"
+set ylabel "$\\sum$(Kinetic Energy) [J]"
+set offsets graph 0.1, graph 0.1, graph 0.1, graph 0.1
+set key outside top horizontal center
+
+plot 'kinetic_without_slowdown' using 1:2 with linespoints lt 1 title "No External Force", \
+    'kinetic_with_slowdown' using 1:($2-4.50993810539) with linespoints lt 2 title "External Force"
diff --git a/class/ke_plot.gnu b/class/ke_plot.gnu
new file mode 100644
index 0000000..c331f1c
--- /dev/null
+++ b/class/ke_plot.gnu
@@ -0,0 +1,14 @@
+#set term cairolatex pdf size 5in,3in
+#set output "kinetic_energy.tex"
+set term x11 persist
+set linetype 1 lc rgb '#E41A1C' pt 7 ps 0.3
+set linetype 2 lc rgb '#377EB8' pt 7 ps 0.3
+set linetype 3 lc rgb '#4DAF4A' pt 7 ps 0.3
+set linetype 4 lc rgb '#984EA3' pt 7 ps 0.3
+
+set xlabel "Time [s]"
+set ylabel "$\\sum$(Kinetic Energy) [J]"
+set offsets graph 0.1, graph 0.1, graph 0.1, graph 0.1
+set nokey
+
+plot 'kinetic_energies.out' using 1:2 with linespoints lt 1 notitle
diff --git a/class/langdyn.f90 b/class/langdyn.f90
new file mode 100644
index 0000000..b1880b2
--- /dev/null
+++ b/class/langdyn.f90
@@ -0,0 +1,282 @@
+module globals
+! Global variables
+use omp_lib                                    ! help the compiler find the OMP libraries
+implicit none
+integer :: n=1000
+double precision, parameter :: L=1.0d0
+double precision, parameter :: pi=2q0*asin(1q0) ! numerical constant
+end module globals
+
+module Langevin
+! Initialization and update rule for Langevin particles
+use globals
+implicit none
+logical, allocatable, dimension(:) :: is_tracked
+double precision :: dt,kT,g,m,sigma,eps,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.00001d0 ! time step size
+kT=0.1d0    ! energy
+g=0.1d0     ! drag coefficient
+m=1d0     ! mass of the particles, can be normalized to 1.
+sigma=1d-3              ! Potential parameters
+eps=1d0
+rc=sigma*2d0**(1d0/6d0) ! Effective particle size
+
+! Set auxiliary parameters
+pref1=g
+pref2=sqrt(24d0*kT*g/dt)
+
+end subroutine set_parameters
+subroutine initialize_particles
+integer :: i, j, k, nx, ix, iy
+double precision :: ran1(n),ran2(n),gr1(n),gr2(n),dist, lx, ly
+logical, allocatable, dimension(:,:) :: is_filled
+! Give particles initial position and velocity
+
+   call random_seed()
+   call random_number(ran1)                       ! uses the built-in PRNG, easy but not very accurate
+   call random_number(ran2)
+   !x=L*(ran1-0.5d0)
+   !x0=x
+   !y=L*(ran2-0.5d0)
+   !y0=y
+
+   !> Basic box stacking with no RNG
+   !nx = ceiling(sqrt(real(n)))
+   !dist = (0.9*L)/nx
+   !k = 1
+   !outer: do i = 1,nx
+   ! do j = 1,nx
+   !     x(k) = i*dist-L/2
+   !     y(k) = j*dist-L/2
+   !     k = k + 1
+   !     if (k.eq.n) exit outer
+   ! end do
+   !end do outer
+   
+   !> Box RNG on a finer mesh 
+   nx = ceiling(L/rc)
+   dist = 0.9*L/nx
+
+   allocate(is_filled(nx,nx))
+   k = 1
+
+   !> Assuming (nx)*(nx) grid, first scale random numbers to a max of nx
+   do while (k .lt. n)
+       !> Random numbers
+       call random_number(lx)
+       call random_number(ly)
+
+       !> Scale them to be grid points
+       ix = ceiling(lx*nx)
+       iy = ceiling(ly*nx)
+
+       !> Check if grid is occupied
+       if (is_filled(ix,iy)) cycle
+
+       !> If not, set the particle position and call it a day
+       x(k) = ix*rc - L/2
+       y(k) = iy*rc - L/2
+
+       !> And mark the spot as filled
+       is_filled(ix,iy) = .true.
+
+       !> Increment!!
+       k = k + 1
+   end do
+
+   ax=0d0
+   ay=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=gr1
+   vy=gr2
+
+end subroutine initialize_particles
+end module Langevin
+
+module domainDecomposition
+  use globals
+  use Langevin
+  implicit none
+  integer, parameter :: b=4
+  integer :: nbl(0:b*b-1,9)
+contains
+  include "neighbourList.f90"
+  include "putParticleInBox.f90"
+  include "sortParticles.f90"
+end module domainDecomposition
+
+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
+
+     if(abs(x(i)).gt.L/2 .or. abs(y(i)).gt.L/2) is_tracked(i)=.false.
+
+   end subroutine impose_BC
+end module BC
+
+program main
+  use globals
+  use domainDecomposition
+use Langevin
+use BC
+implicit none
+integer :: i,j,lim(0:b*b,2),s,ns,p1,p2
+double precision :: t,t_max,m1,m2,rx,ry,dij,F
+double precision :: wtime,begin,end,Tint,TsinceWrite
+double precision, allocatable, dimension(:) :: ran1,ran2,xscrap,yscrap,vxscrap,vyscrap
+
+! Open files
+open(11,file='trajectories.out')
+open(12,file='means.out')
+open(13,file='kinetic_energies.out')
+
+! Allocate arrays
+allocate(x(n),y(n),vx(n),vy(n),ax(n),ay(n),vhx(n),vhy(n),x0(n),y0(n),is_tracked(n),ran1(n),ran2(n),xscrap(n),yscrap(n),vxscrap(n),vyscrap(n))
+
+call buildNBL()
+
+is_tracked = .True.
+t=0d0
+t_max=2.0d0     ! integration time
+Tint = 0.001d0
+TsinceWrite=0d0
+call set_parameters
+call initialize_particles
+
+begin = omp_get_wtime()
+!call cpu_time(begin)
+
+! Conclusion: we need to re-order the loop like this:
+! a. update half-step velocities
+! b. update positions
+! c. compute accellerations/forces
+! d. update all velocities
+
+!$omp parallel
+do while(t.lt.t_max)
+   ! one thread: fetch pseudo-random numbers
+   ! one thread: update velocity, position, impose BC
+   !$omp sections
+   !$omp section  ! Write to disk
+       if(tSinceWrite.gt.Tint) then
+       !if(.true.) then
+         do i=1,n
+            write(11,*) xscrap(i),yscrap(i)
+         enddo
+         write(11,*) NEW_LINE('A')
+         write(12,*) t,sqrt(sum((xscrap-x)**2+(yscrap-y)**2, mask=is_tracked)/real(count(is_tracked),8))
+         write(13,*) t,sum(0.5*vx**2 + 0.5*vy**2, mask=is_tracked)
+         tSinceWrite=0d0
+      end if
+       xscrap=x
+       yscrap=y
+       vxscrap=vx
+       vyscrap=vy
+       !  write(12,*) t, sum(m*(vxscrap**2+vyscrap**2)/(2*n))
+
+   !$omp section
+      call random_number(ran1)
+      ran1=ran1-0.5d0
+      call random_number(ran2)
+      ran2=ran2-0.5d0
+
+   !$omp section
+      vhx=vx+0.5d0*ax*dt
+      vhy=vy+0.5d0*ay*dt
+      x=x+vhx*dt
+      y=y+vhy*dt
+      do j=1,n
+         call impose_BC(j)
+      end do
+      call order(x,y,vx,vy,x0,y0,vhx,vhy,is_tracked,lim) ! BUGFIX 11/03: half-step velocities must also be re-ordered.
+      ax=0d0                   ! Add forces here if any
+      ay=0d0                   ! Add forces here if any
+   !$omp end sections
+   !$omp flush
+   !$omp single
+      ax=ax-pref1*vhx+pref2*ran1
+      ay=ay-pref1*vhy+pref2*ran2
+   !$omp end single
+
+! Our first attempt at parallelization of the code: run the computation of the distances and interaction forces on multiple threads:
+!$omp do private(s,i,ns,p1,p2,rx,ry,dij,F)
+   do s=0,b*b-1
+      do i=1,9
+         if(nbl(s,i).eq.-1) exit
+         ns=nbl(s,i)  ! BUGFIX 11/03: the loop counter was used as sector index instead of the entries of nbl.
+         do p1=lim(s,1),lim(s,2)
+            do p2=lim(ns,1),lim(ns,2)
+               if(p1.eq.p2) cycle
+               rx=x(p2)-x(p1)
+               ry=y(p2)-y(p1)
+               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(p1)=ax(p1)+F*rx/(dij*m)
+                  ay(p1)=ay(p1)+F*ry/(dij*m)
+               end if
+            end do
+         end do
+      end do
+   end do
+   !$omp end do
+   !$omp single
+      vx=vhx+0.5d0*ax*dt
+      vy=vhy+0.5d0*ay*dt
+      t=t+dt
+      tSinceWrite=tSinceWrite+dt
+   !$omp end single
+end do
+!$omp end parallel
+
+end = omp_get_wtime()
+!call cpu_time(end)
+print *,'Wtime=',end-begin
+if(lim(b*b,2).gt.lim(b*b,1)) print '(a5,i7,a8,i8,a11)','Lost ',lim(b*b,2)-lim(b*b,1),' out of ',n,' particles.'
+! De-allocate arrays
+deallocate(x,y,vx,vy,ax,ay,x0,y0,is_tracked,ran1,ran2,xscrap,yscrap,vxscrap,vyscrap)
+! Close files
+close(11)
+close(12)
+close(13)
+
+end program main
diff --git a/class/makefile b/class/makefile
deleted file mode 100644
index 935d1c6..0000000
--- a/class/makefile
+++ /dev/null
@@ -1,19 +0,0 @@
-FC = gfortran
-FFLAGS = -Wall -Wextra
-
-all: test
-
-mod_globals.o: mod_globals.f90
-	$(FC) $(FFLAGS) -c mod_globals.f90
-
-mod_sectors.o: mod_sectors.f90
-	$(FC) $(FFLAGS) -c mod_sectors.f90
-
-test_sectors: mod_sectors_test.f90 mod_globals.o mod_sectors.o
-	$(FC) $(FFLAGS) -o test_sectors.x mod_sectors_test.f90 mod_globals.o mod_sectors.o
-	
-tests: test_sectors
-	./test_sectors.x
-
-clean:
-	rm -f *.o *.x *.mod
diff --git a/class/mean_plot.gnu b/class/mean_plot.gnu
new file mode 100644
index 0000000..e86137d
--- /dev/null
+++ b/class/mean_plot.gnu
@@ -0,0 +1,20 @@
+#set term cairolatex pdf size 5in,3in
+#set output "kinetic_energy.tex"
+set term x11 persist
+set linetype 1 lc rgb '#E41A1C' pt 7 ps 0.3
+set linetype 2 lc rgb '#377EB8' pt 7 ps 0.3
+set linetype 3 lc rgb '#4DAF4A' pt 7 ps 0.3
+set linetype 4 lc rgb '#984EA3' pt 7 ps 0.3
+
+set xlabel "Time [s]"
+set ylabel "RMSD [m]"
+set logscale xy
+set offsets graph 0.1, graph 0.1, graph 0.1, graph 0.1
+set nokey
+
+a=0
+cumsum(x)=(a=a+x,a)
+
+plot 'means.out' using 1:(cumsum($2)) with linespoints lt 1 notitle, \
+    'means.out' with linespoints notitle, \
+    x**2 with lines dt 2
diff --git a/class/mod_globals.f90 b/class/mod_globals.f90
deleted file mode 100644
index 7717c88..0000000
--- a/class/mod_globals.f90
+++ /dev/null
@@ -1,8 +0,0 @@
-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
diff --git a/class/mod_sectors_test.f90 b/class/mod_sectors_test.f90
deleted file mode 100644
index 1f490df..0000000
--- a/class/mod_sectors_test.f90
+++ /dev/null
@@ -1,38 +0,0 @@
-program sectors_test
-
-    use, intrinsic :: iso_fortran_env
-    use :: sectors
-    implicit none
-    
-    integer(int32), allocatable, dimension(:) :: neighbours_list
-
-    !> Testing get_neighbor_ids function
-    !> list = get_neighbor_ids(p,N)
-    !> input: p,n
-    !> Output: variable-length vector
-
-    print *, "Neighbours of 0:"
-    neighbours_list = get_neighbour_ids(0,4)
-    print *, neighbours_list
-
-    print *, "Neighbours of 1:"
-    neighbours_list = get_neighbour_ids(1,4)
-    print *, neighbours_list
-    
-    print *, "Neighbours of 5:"
-    neighbours_list = get_neighbour_ids(5,4)
-    print*, neighbours_list
-
-    print *, "Neighbours of 12:"
-    neighbours_list = get_neighbour_ids(12,4)
-    print*, neighbours_list
-
-    print *, "Neighbours of 13:"
-    neighbours_list = get_neighbour_ids(13,4)
-    print*, neighbours_list
-
-    print *, "Neighbours of 15:"
-    neighbours_list = get_neighbour_ids(15,4)
-    print*, neighbours_list
-
-end program sectors_test
diff --git a/class/movie_plot.gnu b/class/movie_plot.gnu
new file mode 100644
index 0000000..683753c
--- /dev/null
+++ b/class/movie_plot.gnu
@@ -0,0 +1,20 @@
+set term x11 noraise
+set xrange [-0.6:0.6]
+set yrange [-0.6:0.6]
+set size square
+set key outside top center
+set object rectangle from -0.5,-0.5 to 0.5,0.5 dt 3
+
+stats 'trajectories.out' nooutput
+
+do for [i=0:STATS_blocks-1] {
+    plot 'trajectories.out' index i with points pt 7 ps 0.5 title sprintf("Frame %i out of %i", i, STATS_blocks-1)
+    pause 0.0001
+}
+
+set term gif size 500,500 animate delay 0.001
+set output "wiggly.gif"
+
+do for [i=0:STATS_blocks-1:5] {
+    plot 'trajectories.out' index i with points pt 7 ps 0.5 title sprintf("Frame %i out of %i", i, STATS_blocks-1)
+}
diff --git a/class/mod_sectors.f90 b/class/neighbourList.f90
index d32a75a..8bf9c14 100644
--- a/class/mod_sectors.f90
+++ b/class/neighbourList.f90
@@ -1,29 +1,22 @@
-module sectors
-use globals
-implicit none
-contains
-
-subroutine ONETO2(id, box_size, x, y)
+subroutine ONETO2(id, 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)
+    y = id / b
+    x = mod(id, b)
 end subroutine ONETO2
 
-function TWOtoONE(x, y, box_size) result(id)
+function TWOtoONE(x, y) result(id)
     implicit none
-    integer, intent(in) :: x, y, box_size
+    integer, intent(in) :: x, y
     integer :: id
 
-    id = y * box_size + x
+    id = y * b + x
 end function TWOtoONE
 
-function get_neighbour_ids(p, N) result (neighbours)
+function get_neighbour_ids(p) 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
     integer :: x_cell, y_cell, x_test, y_test
@@ -33,7 +26,7 @@ function get_neighbour_ids(p, N) result (neighbours)
     k = 1
     max_list = -1
 
-    call ONETO2(p,N,x_cell,y_cell)
+    call ONETO2(p,x_cell,y_cell)
     max_list(k) = p
 
     !print *, max_list
@@ -48,18 +41,23 @@ function get_neighbour_ids(p, N) result (neighbours)
             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
+            if((x_test .ge. 0 .and. x_test .lt. b) .and. (y_test .ge. 0 .and. y_test .lt. b)) 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)
+                max_list(k) = TWOtoONE(x_test, y_test)
             end if
         end do
     end do
 
     !> Return a vector of all non-zero elements
-    neighbours = pack(max_list, max_list>= 0.0)
-
+!    neighbours = pack(max_list, max_list>= 0.0)
+    neighbours = max_list ! entries -1 must be skipped as they correspond to non-existent neighbours at the edge or corner
 end function get_neighbour_ids
     
-end module sectors
+subroutine buildNBL()
+  integer :: sector
+  do sector=0,b*b-1
+     nbl(sector,:) = get_neighbour_ids(sector)
+  end do
+end subroutine buildNBL
diff --git a/class/putParticleInBox.f90 b/class/putParticleInBox.f90
new file mode 100644
index 0000000..d98a116
--- /dev/null
+++ b/class/putParticleInBox.f90
@@ -0,0 +1,33 @@
+   ! intakes coordinates and outputs index
+   integer function indexfxn(px, py) result(idx)
+     double precision, intent(in) :: px, py
+     double precision :: boxSize,yCoord
+      integer :: ix, iy
+
+     !> Check if particle is in box
+     if(abs(px).GT.L/2 .OR. abs(py).GT.L/2) then
+         !> Particle is outside, don't track it
+         idx=b*b ! This is the "recycle bin" for untracked particles.
+     else
+      ! the sector size      
+      ! M = floor(L / rc)
+      yCoord = - py
+      ! coordinate + half length of boundary [-L/2,L/2] -> [0,L] 
+      ! divided by effective particle size to give coordinate
+      boxSize=L/real(b,8)
+      ix = floor((px + L/2.0d0)/boxSize)
+      iy = floor((yCoord + L/2.0d0)/boxSize)
+
+      ! for particles on the boundary
+      ix = min(max(ix, 0), b-1)
+      iy = min(max(iy, 0), b-1)
+      if (ix < 0 .or. ix >= b .or. iy < 0 .or. iy >= b) then
+         print *, "ERROR: particle outside sector grid"
+         print *, "px, py =", px, py
+         stop
+      end if
+
+      ! call rey's function
+      idx = TWOtoONE(ix,iy)
+     end if
+   end function indexfxn
diff --git a/class/sortParticles.f90 b/class/sortParticles.f90
new file mode 100644
index 0000000..1645b5d
--- /dev/null
+++ b/class/sortParticles.f90
@@ -0,0 +1,56 @@
+subroutine makeNR(nr,in)
+  integer :: particle,sector,nr(0:b*b),in(n)
+  nr = 0
+  do particle=1,n
+     sector=indexfxn(x(particle),y(particle))
+     if(sector.eq.b*b) is_tracked(particle)=.False.
+     in(particle)=sector
+     nr(sector)=nr(sector)+1
+  end do
+  return
+end subroutine makeNR
+
+subroutine order(x,y,vx,vy,x0,y0,vhx,vhy,is_tracked,lim)
+! Order the list of positions by sector and find starting and final index for each sector
+! In: x and y coordinates and velocities. Out: ordered lists x, y, vx and vy and array lim with one row 
+! for each sector, first column is start index, second is end index so that particles lim(k,1) through lim(k,2) reside in sector k.
+integer :: lim(0:b*b,2),in(n),nr(0:b*b),ct(0:b*b),k
+double precision :: x(n),y(n),d1(n),d2(n),d3(n),d4(n),d5(n),d6(n),d7(n),d8(n),vx(n),vy(n),x0(n),y0(n),vhx(n),vhy(n)
+logical :: d9(n), is_tracked(n)
+
+call makeNR(nr,in)
+
+! Set loop limits based on the number of particles in each sector
+lim(0,1)=1
+lim(0,2)=nr(0)
+
+do k=1,b*b
+   lim(k,1)=lim(k-1,2)+1
+   lim(k,2)=lim(k-1,2)+nr(k)
+end do
+
+! Re-order particle list
+d1=x
+d2=y
+d3=vx
+d4=vy
+d5=x0
+d6=y0
+d7=vhx
+d8=vhy
+d9=is_tracked
+ct=0
+do k=1,n
+   x(lim(in(k),1)+ct(in(k)))=d1(k)
+   y(lim(in(k),1)+ct(in(k)))=d2(k)
+   vx(lim(in(k),1)+ct(in(k)))=d3(k)
+   vy(lim(in(k),1)+ct(in(k)))=d4(k)
+   x0(lim(in(k),1)+ct(in(k)))=d5(k)
+   y0(lim(in(k),1)+ct(in(k)))=d6(k)
+   vhx(lim(in(k),1)+ct(in(k)))=d7(k)
+   vhy(lim(in(k),1)+ct(in(k)))=d8(k)
+   is_tracked(lim(in(k),1)+ct(in(k)))=d9(k)
+   ct(in(k))=ct(in(k))+1
+end do
+
+end subroutine order