diff options
Diffstat (limited to 'class-langevin-motion.f90')
| -rw-r--r-- | class-langevin-motion.f90 | 260 |
1 files changed, 260 insertions, 0 deletions
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 |