Added new version of code supporting vectorization

1 files changed, 0 insertions, 260 deletions

diff --git a/class-langevin-motion.f90 b/class-langevin-motion.f90
deleted file mode 100644
index b52396c..0000000
--- a/class-langevin-motion.f90
+++ /dev/null
@@ -1,260 +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
-
-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