Updated class code

1 files changed, 0 insertions, 192 deletions

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