! molecular dynamics in two dimensions
! Program to accompany "Computational Physics" by N. Giordano and H. Nakanishi
! Copyright Prentice Hall 1997, 2006
!
! Modified by H. Nakanishi for better initialization options and for
! better output viewing

program md
   option nolet
   library "sgfunc.trc"
   library "sglib.trc"
   dim x(0,0),y(0,0),energy(10000),vx(0),vy(0),mytime(10000),sqpos(10000),sqpair(10000)
   dim p(0),px(0),py(0),temperature(10000)
   dim xp(0),yp(0)
   call initialize(x,y,xp,yp,vx,vy,n_particles,len,dt,n_plot,out$)
   if out$ <> "" then
      open #3: name out$, organization text, create newold
      erase #3
   end if
   d2=len/400
   d1=d2*1.05
   t = 0    ! t = time
   open #1: screen 0,1,0,0.96
   open #2: screen 0,1,0.96,1
   window #1
   set window -1,len+1,-1,len+1
   set background color "white"
   clear
   set color "black"
   window #1
   call display_grid(len)
   call display(x,y,xp,yp,n_particles,len,d1,d2)
   window #2
   set background color "white"
   set color "black"
   print #2: "Time =",t
   i = 0
   j = 0
   n_p = 0
   r2 = 0
   x1 = 0
   y1 = 0
   x2 = x(2,2)-x(2,1)
   y2 = y(2,2)-y(2,1)
   do                      ! move forward one time step
      call update(x,y,vx,vy,n_particles,len,dt,x1,y1,r2,x2,y2,rp2)   ! 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
         window #1
         call display(x,y,xp,yp,n_particles,len,d1,d2)
         window #2
         clear
         print #2: "Time =",t;
         j = 0
         i = i + 1
         mytime(i) = t      ! calculate and record time, energy, and temperature
         sqpos(i) = r2
         sqpair(i) = rp2
         call calc_energy(x(,),y(,),vx(),vy(),n_particles,len,e,pot_e,temp)
         energy(i) = e
         temperature(i) = temp
         print #2: ", Temperature =",temp;
         if out$ <> "" then
            print #3, using "#.#####^^^^": t;
            print #3, using "<#":", ";
            print #3, using "#.#####^^^^": e;
            print #3, using "<#":", ";
            print #3, using "#.#####^^^^": temp;
            print #3, using "<#":", ";
            print #3, using "#.#####^^^^": r2;
            print #3, using "<#":", ";
            print #3, using "#.#####^^^^": rp2
         end if
         n_p = n_p + 1
      end if
      if key input then
            get key z
            c$ = chr$(z)
            if c$ = "q" then
                  exit do
            else if c$ = "c" then
                  window #1
                  clear
                  call display_grid(len)
            else if c$ = "p" then
                  get key z
            else if c$ = "+" then
                  for k=1 to n_particles
                      x(1,k)=x(2,k)-1.5*(x(2,k)-x(1,k))
                      y(1,k)=y(2,k)-1.5*(y(2,k)-y(1,k))
                  next k
            else if c$ = "-" then
                  for k=1 to n_particles
                      x(1,k)=x(2,k)-0.5*(x(2,k)-x(1,k))
                      y(1,k)=y(2,k)-0.5*(y(2,k)-y(1,k))
                  next k
            else if c$ = "1" then
                  for k=1 to n_particles
                      x(1,k)=x(2,k)-1.1*(x(2,k)-x(1,k))
                      y(1,k)=y(2,k)-1.1*(y(2,k)-y(1,k))
                  next k
            else if c$ = "2" then
                  for k=1 to n_particles
                      x(1,k)=x(2,k)-1.2*(x(2,k)-x(1,k))
                      y(1,k)=y(2,k)-1.2*(y(2,k)-y(1,k))
                  next k
            else if c$ = "3" then
                  for k=1 to n_particles
                      x(1,k)=x(2,k)-1.3*(x(2,k)-x(1,k))
                      y(1,k)=y(2,k)-1.3*(y(2,k)-y(1,k))
                  next k
            else if c$ = "4" then
                  for k=1 to n_particles
                      x(1,k)=x(2,k)-1.4*(x(2,k)-x(1,k))
                      y(1,k)=y(2,k)-1.4*(y(2,k)-y(1,k))
                  next k
            end if
      end if
   loop 
      ! now finished      prepare for plotting final results
   mat redim mytime(i), energy(i), temperature(i), sqpos(i), sqpair(i)
   close #1
   close #2
   if out$ <> "" then close #3
   set background color "white"
   call setcanvas("white")
   set color "black"
   call settitle("Energy time series")
   call sethlabel("Time (in L-J units)")
   call setvlabel("Total energy (in epsilon)")
   call datagraph(mytime,energy,0,1,"black")
   set cursor 1,1
   print "hit any key to proceed -> "
   get key z
   clear
   call settitle("Temperature time series")
   call sethlabel("Time (in L-J units)")
   call setvlabel("Temperature (in epsilon/k_B)")
   call datagraph(mytime,temperature,0,1,"black")
   set cursor 1,1
   print "hit any key to proceed -> "
   get key z
   clear
   call settitle("Tagged particle square displacement time series")
   call sethlabel("Time (in L-J units)")
   call setvlabel("r^2 (in sigma^2)")
   call datagraph(mytime,sqpos,0,1,"black")
   print "hit any key to proceed -> "
   get key z
   clear
   call settitle("Tagged pair square distance time series")
   call sethlabel("Time (in L-J units)")
   call setvlabel("dr^2 (in sigma^2)")
   call datagraph(mytime,sqpair,0,1,"black")
end

! initialize variables
sub initialize(x(,),y(,),xp(),yp(),vx(),vy(),n_particles,len,dt,n_plot,out$)
   input prompt "# particles => ": n_particles
   mat redim x(2,n_particles),y(2,n_particles), vx(n_particles), vy(n_particles)
   mat redim xp(n_particles),yp(n_particles)
   input prompt "side length (in units of sigma) => ": len
   input prompt "time step (in units given on p.235) => ": dt
   input prompt "plot/record after this many steps => ": n_plot
   input prompt "max initial speed => ": vmax
   input prompt "max initial relative displacement => ": dmax
   line input prompt "file name to save thermodynamic data => ": out$
   print "Type c (clear green tracks), p to pause, or q to plot time series and end." 
   print "Type + to increase speed 50% or - to decrease 50%, or";
   print "type N [N=1,...,4] to increase speed by 10N%."
   print "Hit any key to start..."
   get key z
   randomize
   sqn = sqr(n_particles)
   if (sqn = int(sqn)) then
      grid = len/sqn
   else 
      grid = len / int(sqn + 1)
   end if
   n = 0
   i = 0
   do while i < len        ! arrange particles on a roughly square lattice
      j = 0			! to keep them apart initially
      do while j < len
         n = n + 1
         if n <= n_particles then 
            x(2,n) = i + 0.5*grid + dmax * (rnd - 0.5) * grid * sqr(2)
            y(2,n) = j + 0.5*grid + dmax * (rnd - 0.5) * grid * sqr(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
end sub

! initialize velocities randomly
sub init_velocity(vx,vy,v0)
   vx = v0 * (rnd - 0.5) * sqr(2)
   vy = v0 * (rnd - 0.5) * sqr(2)
end sub

! 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,x1,y1,r2,x2,y2,rp2)
   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 i=1 then
         x1 = x1+x_new(1)-x(2,1)
         y1 = y1+y_new(1)-y(2,1)
         r2 = x1^2+y1^2
      else if i=2 then
         x2 = x2+x_new(2)-x(2,2)
         y2 = y2+y_new(2)-y(2,2)
         rp2 = (x2-x1)^2+(y2-y1)^2
      end if
      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

! 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 background color "white"
   set color "black"
   clear
   box lines 0,len,0,len
   set color "red"
end sub

! display particles
sub display(x(,),y(,),xp(),yp(),n_particles,len,d1,d2)
   for i = 1 to n_particles
      set color "green"
      box area xp(i)-d1,xp(i)+d1,yp(i)-d1,yp(i)+d1
      if i=1 then
         set color "blue"
      else if i=2 then 
         set color "black"
      else
         set color "red"
      end if
      box area x(2,i)-d2,x(2,i)+d2,y(2,i)-d2,y(2,i)+d2
      xp(i) = x(2,i)
      yp(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