PROGRAM md
PUBLIC x(36),y(36),vx(36),vy(36),ax(36),ay(36)
PUBLIC N,Lx,Ly,dt,dt2
PUBLIC gcum(0:1000),nbin,dr
PUBLIC xsave(100),ysave(100),R2cum(100)
LIBRARY "csgraphics"
CALL initial(t,ke,kecum,pecum,vcum,area)
CALL set_up_windows(#1,#2)
CALL accel(pe,virial)
LET E = ke + pe                   ! total energy
LET ncum = 0                      ! number of times data accumulated
LET flag$ = ""
DO
   CALL show_positions(flag$,#2)
   CALL Verlet(t,ke,pe,virial)
   CALL show_output(t,ke,pe,virial,kecum,vcum,ncum,area,#1)
LOOP until flag$ = "stop"
CALL save_config
END

SUB initial(t,ke,kecum,pecum,vcum,area)
    DECLARE PUBLIC x(),y(),vx(),vy()
    DECLARE PUBLIC N,Lx,Ly,dt,dt2
    DECLARE PUBLIC gcum(),nbin,dr
    DECLARE DEF pbc
    LET dt = 0.01
    LET dt2 = dt*dt
    LET response$ = ""
    DO
       INPUT prompt "read data statements (d) or file (f)? ": start$
       IF (start$ = "d") or (start$ = "D") then
          LET response$ = "ok"
          LET N = 16
          LET Lx = 6
          LET Ly = 6
          DATA 1.09,0.98,-0.33,0.78,3.12,5.25,0.12,-1.19
          DATA 0.08,2.38,-0.08,-0.10,0.54,4.08,-1.94,-0.56
          DATA 2.52,4.39,0.75,0.34,3.03,2.94,1.70,-1.08
          DATA 4.25,3.01,0.84,0.47,0.89,3.11,-1.04,0.06
          DATA 2.76,0.31,1.64,1.36,3.14,1.91,0.38,-1.24
          DATA 0.23,5.71,-1.58,0.55,1.91,2.46,-1.55,-0.16
          DATA 4.77,0.96,-0.23,-0.83,5.10,4.63,-0.31,0.65
          DATA 4.97,5.88,1.18,1.48,3.90,0.20,0.46,-0.51
          FOR i = 1 to N
              READ x(i),y(i),vx(i),vy(i)
          NEXT i
       ELSE IF (start$ = "f") or (start$ = "F") then
          LET response$ = "ok"
          INPUT prompt "file name = ": file$
          OPEN #1: name file$, access input
          INPUT #1: N
          INPUT #1: Lx
          INPUT #1: Ly
          INPUT #1: heading$
          FOR i = 1 to N
              INPUT #1: x(i),y(i)
          NEXT i
          INPUT #1: heading$
          FOR i = 1 to N
              INPUT #1: vx(i),vy(i)
          NEXT i
          CLOSE #1
       ELSE
          PRINT
          PRINT "d or f are the only acceptable responses."
       END IF
    LOOP until response$ = "ok"
    CLEAR
    LET ke = 0                    ! kinetic energy
    FOR i = 1 to N
        LET ke = ke + vx(i)*vx(i) + vy(i)*vy(i)
    NEXT i
    LET ke = 0.5*ke
    LET area = Lx*Ly
    LET t = 0                     ! time
    ! initialize sums
    LET kecum = 0
    LET pecum = 0
    LET vcum = 0
END SUB

SUB set_up_windows(#1,#2)
    DECLARE PUBLIC Lx,Ly
    OPEN #1: screen 0,1,0.90,1.0  ! numerical output
    OPEN #2: screen 0.02,1,0.02,0.90   ! particle trajectories
    CALL compute_aspect_ratio(Lx,xwin,ywin)
    SET WINDOW 0,xwin,0,ywin
    BOX LINES 0,Lx,0,Ly
    CALL headings(#1)
END SUB

SUB headings(#1)
    WINDOW #1
    SET CURSOR 1,1
    PRINT using "##########": "time steps";
    PRINT,
    PRINT using "###.##": "time";
    PRINT,
    PRINT using "---#.###": "energy";
    PRINT,
    PRINT using "##.##": "<T>";
    PRINT,
    PRINT using "##.##": "<P>"
END SUB

SUB Verlet(t,ke,pe,virial)
    DECLARE PUBLIC x(),y(),vx(),vy(),ax(),ay()
    DECLARE PUBLIC N,Lx,Ly,dt,dt2
    DECLARE DEF pbc
    FOR i = 1 to N
        LET xnew = x(i) + vx(i)*dt + 0.5*ax(i)*dt2
        LET ynew = y(i) + vy(i)*dt + 0.5*ay(i)*dt2
        ! partially update velocity using old acceleration
        LET vx(i) = vx(i) + 0.5*ax(i)*dt
        LET vy(i) = vy(i) + 0.5*ay(i)*dt
        LET x(i) = pbc(xnew,Lx)
        LET y(i) = pbc(ynew,Ly)
    NEXT i
    CALL accel(pe,virial)         ! new acceleration
    LET ke = 0
    FOR i = 1 to N
        ! complete the update of the velocity using new acceleration
        LET vx(i) = vx(i) + 0.5*ax(i)*dt
        LET vy(i) = vy(i) + 0.5*ay(i)*dt
        LET ke = ke + vx(i)*vx(i) + vy(i)*vy(i)
    NEXT i
    LET ke = 0.5*ke
    LET t = t + dt
END SUB

FUNCTION pbc(pos,L)
    IF pos < 0 then
       LET pbc = pos + L
    ELSE IF pos > L then
       LET pbc = pos - L
    ELSE
       LET pbc = pos
    END IF
END DEF

SUB accel(pe,virial)
    DECLARE PUBLIC x(),y(),ax(),ay()
    DECLARE PUBLIC N,Lx,Ly
    DECLARE DEF separation
    FOR i = 1 to N
        LET ax(i) = 0
        LET ay(i) = 0
    NEXT i
    LET pe = 0
    LET virial = 0
    FOR i = 1 to N - 1           ! compute total force on particle i
        FOR j = i + 1 to N        ! due to particles j > i
            LET dx = separation(x(i) - x(j),Lx)
            LET dy = separation(y(i) - y(j),Ly)
            ! acceleration = force because mass = 1 in reduced units
            CALL force(dx,dy,fxij,fyij,pot)
            LET ax(i) = ax(i) + fxij
            LET ay(i) = ay(i) + fyij
            LET ax(j) = ax(j) - fxij   ! Newton's third law
            LET ay(j) = ay(j) - fyij
            LET pe = pe + pot
            LET virial = virial + dx*fxij + dy*fyij
        NEXT j
    NEXT i
END SUB

SUB force(dx,dy,fx,fy,pot)
    LET r2 = dx*dx + dy*dy
    LET rm2 = 1/r2
    LET rm6 = rm2*rm2*rm2
    LET f_over_r = 24*rm6*(2*rm6 - 1)*rm2
    LET fx = f_over_r*dx
    LET fy = f_over_r*dy
    LET pot = 4*(rm6*rm6 - rm6)
END SUB

FUNCTION separation(ds,L)
    IF ds > 0.5*L then
       LET separation = ds - L
    ELSE IF ds < -0.5*L then
       LET separation = ds + L
    ELSE
       LET separation = ds
    END IF
END DEF

SUB show_positions(flag$,#2)
    DECLARE PUBLIC x(),y(),N,Lx,Ly
    IF key input then
       GET KEY k
       IF k = ord("r") then
          WINDOW #2
          CLEAR
          SET COLOR "black"
          BOX LINES 0,Lx,0,Ly
          LET flag$ = ""
       ELSEIF k = ord("s") then
          LET flag$ = "stop"
       ELSEIF k = ord("n") then
          LET flag$ = "no_show"
       END IF
    END IF
    IF flag$ <> "no_show" then
       WINDOW #2
       SET COLOR "red"
       FOR i = 1 to N
           PLOT x(i),y(i)
       NEXT i
    END IF
END SUB

SUB show_output(t,ke,pe,virial,kecum,vcum,ncum,area,#1)
    WINDOW #1
    DECLARE PUBLIC N,Lx,Ly
    SET CURSOR 2,1
    SET COLOR "black/white"
    LET ncum = ncum + 1
    PRINT using "######": ncum;
    PRINT,
    PRINT using "###.##": t;      ! time
    PRINT,
    LET E = ke + pe               ! total energy
    PRINT using "---#.###": E;
    PRINT,
    LET kecum = kecum + ke
    LET vcum = vcum + virial
    LET mean_ke = kecum/ncum      ! still need to divide by N
    LET p = mean_ke + (0.5*vcum)/ncum  ! mean pressure * area
    LET p = p/area
    PRINT using "##.##": mean_ke/N;    ! mean kinetic temperature
    PRINT,
    PRINT using "##.##": p;
END SUB

SUB save_config
    DECLARE PUBLIC x(),y(),vx(),vy()
    DECLARE PUBLIC N,Lx,Ly
    INPUT prompt "file name of saved configuration = ": file$
    OPEN #1: name file$,access output,create new
    PRINT #1: N
    PRINT #1: Lx
    PRINT #1: Ly
    PRINT #1: "x(i)","y(i)"
    ! comma added between outputs on the same line so that file
    ! can be read by True BASIC
    FOR i = 1 to N
        PRINT #1, using "----.####, ----.#####": x(i),y(i)
    NEXT i
    PRINT #1: "vx(i)","vy(i)"
    FOR i = 1 to N
        PRINT #1, using "----.####, ----.#####": vx(i),vy(i)
    NEXT i
    CLOSE #1
END SUB