New report and code for final submission.HEADmaster

9 files changed, 0 insertions, 522 deletions

diff --git a/class/Makefile b/class/Makefile
deleted file mode 100644
index 17979f9..0000000
--- a/class/Makefile
+++ /dev/null
@@ -1,19 +0,0 @@
-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/ke_compare.gnu b/class/ke_compare.gnu
deleted file mode 100644
index e09d021..0000000
--- a/class/ke_compare.gnu
+++ /dev/null
@@ -1,15 +0,0 @@
-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
deleted file mode 100644
index c331f1c..0000000
--- a/class/ke_plot.gnu
+++ /dev/null
@@ -1,14 +0,0 @@
-#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
deleted file mode 100644
index b1880b2..0000000
--- a/class/langdyn.f90
+++ /dev/null
@@ -1,282 +0,0 @@
-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/mean_plot.gnu b/class/mean_plot.gnu
deleted file mode 100644
index e86137d..0000000
--- a/class/mean_plot.gnu
+++ /dev/null
@@ -1,20 +0,0 @@
-#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/movie_plot.gnu b/class/movie_plot.gnu
deleted file mode 100644
index 683753c..0000000
--- a/class/movie_plot.gnu
+++ /dev/null
@@ -1,20 +0,0 @@
-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/neighbourList.f90 b/class/neighbourList.f90
deleted file mode 100644
index 8bf9c14..0000000
--- a/class/neighbourList.f90
+++ /dev/null
@@ -1,63 +0,0 @@
-subroutine ONETO2(id, x, y)
-    implicit none
-    integer, intent(in)  :: id
-    integer, intent(out) :: x, y
-
-    y = id / b
-    x = mod(id, b)
-end subroutine ONETO2
-
-function TWOtoONE(x, y) result(id)
-    implicit none
-    integer, intent(in) :: x, y
-    integer :: id
-
-    id = y * b + x
-end function TWOtoONE
-
-function get_neighbour_ids(p) result (neighbours)
-    integer, intent(in) :: p !> The sequential position
-    integer, allocatable, dimension(:) :: neighbours !> The list of neighbours
-    integer :: i,j,k
-    integer :: x_cell, y_cell, x_test, y_test
-    integer :: max_list(9)
-
-    !> Don't save the values between calls
-    k = 1
-    max_list = -1
-
-    call ONETO2(p,x_cell,y_cell)
-    max_list(k) = p
-
-    !print *, max_list
-
-    !> 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. 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)
-            end if
-        end do
-    end do
-
-    !> Return a vector of all non-zero elements
-!    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
-    
-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
deleted file mode 100644
index d98a116..0000000
--- a/class/putParticleInBox.f90
+++ /dev/null
@@ -1,33 +0,0 @@
-   ! 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
deleted file mode 100644
index 1645b5d..0000000
--- a/class/sortParticles.f90
+++ /dev/null
@@ -1,56 +0,0 @@
-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