! molecular dynamics in two dimensions
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program md
   option nolet
   library "sgfunc*","sglib*"
   dim x(0,0), y(0,0), energy(10000), vx(0), vy(0), mytime(10000)
   dim p(0),px(0),py(0),temperature(10000)
   call initialize(x,y,vx,vy,n_particles,len,dt,n_plot)
   t = 0    ! t = time
   i = 0
   j = 0
   n_p = 0
   do                      ! move forward one time step
      call update(x,y,vx,vy,n_particles,len,dt)   ! update system
      t = t + dt
      j = j + 1             ! use to keep track of how often to plot on screen
      if j >= n_plot then   ! and record values for later plotting
         call display(x,y,n_particles)
         j = 0
         i = i + 1
         mytime(i) = t      ! calculate and record time, energy, and temperature
         call calc_energy(x(,),y(,),vx(),vy(),n_particles,len,e,pot_e,temp)
         energy(i) = e
         temperature(i) = temp
         n_p = n_p + 1
      end if
      if key input then
            get key z
            c$ = chr$(z)
            if c$ = "q" then exit do
            if c$ = "c" then clear
            if c$ = "p" then 
                  call display_grid(len)
                  call display(x,y,n_particles)
            end if
      end if
   loop 
      ! now finished      prepare for plotting final results
   mat redim mytime(i), energy(i), temperature(i)
   call datagraph(mytime,energy,4,0,"red")
   set cursor 1,1
   print "hit any key to proceed -> "
   get key z
   call datagraph(mytime,temperature,4,0,"red")
end
! move forward one time step
! x(i,n),y(i,n) = position of particle i
! n = 1,2,3 = oldest, current, and new positions
! vx,vy = velocity components   n_particles = number of particles
! len = size of box (use periodic boundary conditions)   dt = time step
sub update(x(,),y(,),vx(),vy(),n_particles,len,dt)
   dim x_new(100),y_new(100)
   for i = 1 to n_particles
      call force(i,x(,),y(,),n_particles,len,fx,fy)  ! compute the forces
      x_new(i) = 2 * x(2,i) - x(1,i) + fx * dt^2     ! use Verlet method
      y_new(i) = 2 * y(2,i) - y(1,i) + fy * dt^2
      vx(i) = (x_new(i) - x(1,i)) / (2 * dt)         ! keep track of velocities
      vy(i) = (y_new(i) - y(1,i)) / (2 * dt)
      if x_new(i) < 0 then    ! periodic boundary conditions
         x_new(i) = x_new(i) + len
         x(2,i) = x(2,i) + len
      else if x_new(i) > len then 
         x_new(i) = x_new(i) - len 
         x(2,i) = x(2,i) - len 
      end if
      if y_new(i) < 0 then 
         y_new(i) = y_new(i) + len
         y(2,i) = y(2,i) + len
      else if y_new(i) > len then 
         y_new(i) = y_new(i) - len 
         y(2,i) = y(2,i) - len 
      end if
   next i
   for i = 1 to n_particles                     ! update current and old values
      x(1,i) = x(2,i)                           
      x(2,i) = x_new(i)
      y(1,i) = y(2,i)
      y(2,i) = y_new(i)
   next i
end sub
! initialize variables
sub initialize(x(,),y(,),vx(),vy(),n_particles,len,dt,n_plot)
   n_particles = 20
   mat redim x(2,n_particles),y(2,n_particles), vx(n_particles), vy(n_particles)
   len = 10
   dt = 0.02
   n_plot = 3    ! plot and record after every third time step
   vmax = 1
   grid = len / int(sqr(n_particles) + 1)
   n = 0
   i = 1
   do while i < len        ! arrange particles on a roughly square lattice
      j = 1                ! to keep them apart initially
      do while j < len
         n = n + 1
         if n <= n_particles then 
            x(2,n) = i + (rnd - 0.5) * grid / 2
            y(2,n) = j + (rnd - 0.5) * grid / 2
            call init_velocity(vx(n),vy(n),vmax)   ! give them all some initial
            x(1,n) = x(2,n) - vx(n) * dt           ! velocity
            y(1,n) = y(2,n) - vy(n) * dt
         end if
         j = j + grid
      loop
      i = i + grid
   loop
   set background color "white"
   set color "black"
   clear
   set window -1,len+1,-1,len+1
   call display_grid(len)
   call display(x,y,n_particles)
end sub
! calculate current energy = potential + kinetic
! also compute temperature via equipartition
sub calc_energy(x(,),y(,),vx(),vy(),n_particles,len,e,pot_e,temp)
   pot_e = 0   ! potential energy
   k_e = 0     ! kinetic energy
   for i = 1 to n_particles
      for j = i+1 to n_particles
         call find_r(i,j,1,x,y,len,r,dx,dy)  ! find spacing of two atoms
         if r < 30 then                      ! using nearest separation rule
            call pair_force(r,f,u)
            pot_e = pot_e + u
         end if
      next j
      k_e = k_e + (vx(i)^2 + vy(i)^2) / 2    ! kinetic energy
   next i
   e = k_e + pot_e
   temp = k_e / n_particles                  ! equipartition in two dimensions
end sub
! display box edges
sub display_grid(len)
   set color "black"
   clear
   box lines 0,len,0,len
   set color "red"
end sub
! display particles
sub display(x(,),y(,),n_particles)
   for i = 1 to n_particles
      plot x(2,i),y(2,i)
   next i 
end sub
! compute forces on all of the particles
sub force(n,x(,),y(,),n_particles,len,fx,fy)
   fx = 0
   fy = 0
   for i = 1 to n_particles
      if i <> n then
         call find_r(i,n,2,x,y,len,r,dx,dy)
         if r < 30 then
            call pair_force(r,f,u)
            fx = fx + f * dx / r
            fy = fy + f * dy / r
         end if
      end if
   next i
end sub
! Lennard-Jones force
sub pair_force(r,f,u)
   u = 4 * (1/r^12 - 1/r^6)
   f = 24 * (2/r^13 - 1/r^7)
end sub
! find spacing taking periodic boundary conditions into account
sub find_r(i,n,flag,x(,),y(,),len,r,dx,dy)
   dx = x(flag,n) - x(flag,i)
   dy = y(flag,n) - y(flag,i)
   if abs(dx) > len / 2 then dx = dx - sgn(dx) * len
   if abs(dy) > len / 2 then dy = dy - sgn(dy) * len
   r = sqr(dx^2 + dy^2)
end sub
! initialize velocities randomly
sub init_velocity(vx,vy,v0)
   vx = v0 * (rnd - 0.5)
   vy = v0 * (rnd - 0.5)
end sub