      PROGRAM md
      IMPLICIT NONE
      REAL*8 E, area, ke, kecum, pe, pecum, t, vcum,
     +       virial, Lx, Ly, R2cum, ax, ay, dr, dt, dt2, gcum,
     +       nbin, vx, vy, x, xsave, y, ysave
      INTEGER IC1_, IC2_, IR_, N, ncum
      CHARACTER*80 flag
      LOGICAL Ltmp1_
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL initial(t,ke,kecum,pecum,vcum,area)
      CALL set_up_windows(IC1_,IC2_)
      CALL accel(pe,virial)
      E = ke + pe                           ! total energy
      ncum = 0                              ! number of times data accumulated
      flag = ''
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL show_positions(flag,IC2_)
         CALL Verlet(t,ke,pe,virial)
         CALL show_output(t,ke,pe,virial,kecum,vcum,ncum,area,IC1_)
         Ltmp1_ = .NOT.(flag .EQ. 'stop')
      END DO
      CALL save_config
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(t,ke,kecum,pecum,vcum,area)
      IMPLICIT NONE
      REAL*8 t, ke, kecum, pecum, vcum, area, Lx, Ly,
     +       R2cum, ax, ay, dr, dt, dt2, gcum, nbin, vx, vy, x, xsave,
     +       y, ysave, DATA_(64)
      INTEGER IO1_, IR_, N, i
      CHARACTER*80 file, heading, response, start
      LOGICAL Ltmp1_
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      INTEGER DP_
      DATA DATA_/
     +    1.09D0,0.98D0,-0.33D0,0.78D0,3.12D0,5.25D0,0.12D0,-1.19D0,
     +    0.08D0,2.38D0,-0.08D0,-0.10D0,0.54D0,4.08D0,-1.94D0,-0.56D0,
     +    2.52D0,4.39D0,0.75D0,0.34D0,3.03D0,2.94D0,1.70D0,-1.08D0,
     +    4.25D0,3.01D0,0.84D0,0.47D0,0.89D0,3.11D0,-1.04D0,0.06D0,
     +    2.76D0,0.31D0,1.64D0,1.36D0,3.14D0,1.91D0,0.38D0,-1.24D0,
     +    0.23D0,5.71D0,-1.58D0,0.55D0,1.91D0,2.46D0,-1.55D0,-0.16D0,
     +    4.77D0,0.96D0,-0.23D0,-0.83D0,5.10D0,4.63D0,-0.31D0,0.65D0,
     +    4.97D0,5.88D0,1.18D0,1.48D0,3.90D0,0.20D0,0.46D0,-0.51D0/
      DATA DP_/1/, IO1_/11/
      dt = 0.01D0
      dt2 = dt*dt
      response = ''
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         WRITE(*,'(A,$)') 'read data statements (d) or file (f)? '
         READ(*,'(A)') start
         IF((start .EQ. 'd') .OR. (start .EQ. 'D')) THEN
            response = 'ok'
            N = 16
            Lx = 6
            Ly = 6
            DO i = 1, N
               x(i) = DATA_(DP_)
               DP_ = DP_ + 1
               y(i) = DATA_(DP_)
               DP_ = DP_ + 1
               vx(i) = DATA_(DP_)
               DP_ = DP_ + 1
               vy(i) = DATA_(DP_)
               DP_ = DP_ + 1
            END DO
         ELSE IF((start .EQ. 'f') .OR. (start .EQ. 'F')) THEN
            response = 'ok'
            WRITE(*,'(A,$)') 'file name = '
            READ(*,'(A)') file
            OPEN(IO1_, file = file)
            READ(IO1_,*) N
            READ(IO1_,*) Lx
            READ(IO1_,*) Ly
            READ(IO1_,'(A)') heading
            DO i = 1, N
               READ(IO1_,*) x(i),y(i)
            END DO
            READ(IO1_,'(A)') heading
            DO i = 1, N
               READ(IO1_,*) vx(i),vy(i)
            END DO
            CLOSE(IO1_)
         ELSE
            WRITE(*,*)
            WRITE(*,*) 'd or f are the only acceptable responses.'
         END IF
         Ltmp1_ = response .NE. 'ok'
      END DO
      CALL GWclear(IR_, -1)
      ke = 0                                ! kinetic energy
      DO i = 1, N
         ke = ke + vx(i)*vx(i) + vy(i)*vy(i)
      END DO
      ke = 0.5D0*ke
      area = Lx*Ly
      t = 0                                 ! time
*     initialize sums
      kecum = 0
      pecum = 0
      vcum = 0
      END
C
      SUBROUTINE set_up_windows(IC1_,IC2_)
      IMPLICIT NONE
      REAL*8 Lx, Ly, R2cum, ax, ay, dr, dt, dt2, gcum,
     +       nbin, vx, vy, x, xsave, y, ysave
      REAL xwin, ywin, GWaspect
      INTEGER IR_, N, IC1_, IC2_
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      CALL GWvport(IR_, 0.0, 0.90, GWaspect(1), 1.0) ! numerical output
      CALL GWindow(IR_, 0.0, 3.0, 80.0, 0.0)
      CALL GWsavevp(IR_, IC1_)
      CALL headings
      CALL GWvport(IR_, 0.02*GWaspect(1), 0.02, GWaspect(1), 0.90)
                                                     ! particle trajectories
      CALL compute_aspect_ratio(REAL(Lx),xwin,ywin)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetpen(IR_, 0, -1, -1, 4)
      CALL GWrect(IR_, 0.0, 0.0, REAL(Lx), REAL(Ly))
      CALL GWanchor(IR_, 1)
      END
C
      SUBROUTINE headings
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      INTEGER IR_
      WRITE(Stmp1_, '(A)') 'time steps '
      CALL GWputtxt(IR_, 1.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') 'time '
      CALL GWputtxt(IR_, 12.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') 'energy '
      CALL GWputtxt(IR_, 24.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') '<T> '
      CALL GWputtxt(IR_, 36.0, 1.0, Stmp1_)
      WRITE(Stmp1_, '(A)') '<P>'
      CALL GWputtxt(IR_, 48.0, 1.0, Stmp1_)
      END
C
      SUBROUTINE Verlet(t,ke,pe,virial)
      IMPLICIT NONE
      REAL*8 xnew, ynew, t, ke, pe, virial, Lx, Ly, R2cum, ax, ay,
     +       dr, dt, dt2, gcum, nbin, vx, vy, x, xsave, y, ysave, pbc
      INTEGER N, i
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      DO i = 1, N
         xnew = x(i) + vx(i)*dt + 0.5D0*ax(i)*dt2
         ynew = y(i) + vy(i)*dt + 0.5D0*ay(i)*dt2
*        partially update velocity using old acceleration
         vx(i) = vx(i) + 0.5D0*ax(i)*dt
         vy(i) = vy(i) + 0.5D0*ay(i)*dt
         x(i) = pbc(xnew,Lx)
         y(i) = pbc(ynew,Ly)
      END DO
      CALL accel(pe,virial)                 ! new acceleration
      ke = 0
      DO i = 1, N
*        complete the update of the velocity using new acceleration
         vx(i) = vx(i) + 0.5D0*ax(i)*dt
         vy(i) = vy(i) + 0.5D0*ay(i)*dt
         ke = ke + vx(i)*vx(i) + vy(i)*vy(i)
      END DO
      ke = 0.5D0*ke
      t = t + dt
      END
C
      FUNCTION pbc(pos,L)
      IMPLICIT NONE
      REAL*8 pbc, pos, L
      IF(pos .LT. 0) THEN
         pbc = pos + L
      ELSE IF(pos .GT. L) THEN
         pbc = pos - L
      ELSE
         pbc = pos
      END IF
      END
C
      SUBROUTINE accel(pe,virial)
      IMPLICIT NONE
      REAL*8 dx, dy, fxij, fyij, pot, pe, virial, Lx, Ly, R2cum, ax,
     +       ay, dr, dt, dt2, gcum, nbin, vx, vy, x, xsave, y, ysave,
     +       separation
      INTEGER N, i, j
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      DO i = 1, N
         ax(i) = 0
         ay(i) = 0
      END DO
      pe = 0
      virial = 0
      DO i = 1, N - 1
         DO j = i + 1, N
            dx = separation(x(i) - x(j),Lx)
            dy = separation(y(i) - y(j),Ly)
*           acceleration = force because mass = 1 in reduced units
            CALL force(dx,dy,fxij,fyij,pot)
            ax(i) = ax(i) + fxij
            ay(i) = ay(i) + fyij
            ax(j) = ax(j) - fxij            ! Newton's third law
            ay(j) = ay(j) - fyij
            pe = pe + pot
            virial = virial + dx*fxij + dy*fyij
         END DO
      END DO
      END
C
      SUBROUTINE force(dx,dy,fx,fy,pot)
      IMPLICIT NONE
      REAL*8 f_over_r, r2, rm2, rm6, dx, dy, fx, fy, pot
      r2 = dx*dx + dy*dy
      rm2 = 1/r2
      rm6 = rm2*rm2*rm2
      f_over_r = 24*rm6*(2*rm6 - 1)*rm2
      fx = f_over_r*dx
      fy = f_over_r*dy
      pot = 4*(rm6*rm6 - rm6)
      END
C
      FUNCTION separation(ds,L)
      IMPLICIT NONE
      REAL*8 separation, ds, L
      IF(ds .GT. 0.5D0*L) THEN
         separation = ds - L
      ELSE IF(ds .LT. -0.5D0*L) THEN
         separation = ds + L
      ELSE
         separation = ds
      END IF
      END
C
      SUBROUTINE show_positions(flag,IC2_)
      IMPLICIT NONE
      REAL*8 Lx, Ly, R2cum, ax, ay, dr, dt, dt2,
     +       gcum, nbin, vx, vy, x, xsave, y, ysave
      INTEGER IR_, N, i, k, IC2_
      LOGICAL TBkeyinput
      CHARACTER*80 flag
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      CALL GWselvp(IR_, IC2_)
      IF(TBkeyinput()) THEN
         CALL TBgetkey(k)
         IF(k .EQ. ICHAR('r')) THEN
            CALL GWerase(IR_, 1, 0)
            CALL GWrefresh(IR_)
            flag = ''
         ELSE IF(k .EQ. ICHAR('s')) THEN
            flag = 'stop'
         ELSE IF(k .EQ. ICHAR('n')) THEN
            flag = 'no_show'
         END IF
      END IF
      IF(flag .NE. 'no_show') THEN
         CALL GWerase(IR_, 1, 0)
         DO i = 1, N
            CALL GWsetpxl(IR_, REAL(x(i)), REAL(y(i)), 13)
         END DO
      END IF
      END
C
      SUBROUTINE show_output(t,ke,pe,virial,kecum,vcum,ncum,area,IC1_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 E, mean_ke, p, t, ke, pe, virial, kecum, vcum, area, Lx,
     +       Ly, R2cum, ax, ay, dr, dt, dt2, gcum, nbin, vx, vy, x,
     +       xsave, y, ysave
      INTEGER IR_, N, ncum, IC1_
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, 2)
      CALL GWerase(IR_, 2, -1)
      ncum = ncum + 1
      WRITE(Stmp1_, '(I6)') ncum
      CALL GWputtxt(IR_, 1.0, 2.0, Stmp1_)
      WRITE(Stmp1_, '(F6.2)') t             ! time
      CALL GWputtxt(IR_, 12.0, 2.0, Stmp1_)
      E = ke + pe                           ! total energy
      WRITE(Stmp1_, '(F8.3)') E
      CALL GWputtxt(IR_, 24.0, 2.0, Stmp1_)
      kecum = kecum + ke
      vcum = vcum + virial
      mean_ke = kecum/ncum                  ! still need to divide by N
      p = mean_ke + (0.5D0*vcum)/ncum       ! mean pressure * area
      p = p/area
      WRITE(Stmp1_, '(F5.2)') mean_ke/N     ! mean kinetic temperature
      CALL GWputtxt(IR_, 36.0, 2.0, Stmp1_)
      WRITE(Stmp1_, '(F5.2,$)') p
      CALL GWputtxt(IR_, 48.0, 2.0, Stmp1_)
      CALL GWsetogn(IR_, 0, N + 1)
      END
C
      SUBROUTINE save_config
      IMPLICIT NONE
      REAL*8 Lx, Ly, R2cum, ax, ay, dr, dt, dt2, gcum, nbin, vx, vy,
     +       x, xsave, y, ysave
      INTEGER IR_, IO1_, N, i
      CHARACTER*80 file
      COMMON /GLBLD/Lx, Ly, R2cum(100), ax(36), ay(36), dr, dt, dt2,
     +             gcum(0:1000), nbin, vx(36), vy(36), x(36),
     +             xsave(100), y(36), ysave(100)
      COMMON /GLBLI/N
      DATA IO1_/11/
      CALL GWinput(IR_, "name of saved configuration = ", file)
      IF(IR_ .EQ. 0) RETURN
      OPEN(IO1_, file = file)
      WRITE(IO1_,*) N
      WRITE(IO1_,*) Lx
      WRITE(IO1_,*) Ly
      WRITE(IO1_,'(2A12)') 'x(i)', 'y(i)'
*     comma added between outputs on the same line so that file
*     can be read by True BASIC
      DO i = 1, N
         WRITE(IO1_, '(3X,F9.4,'', '',F9.5)') x(i), y(i)
      END DO
      WRITE(IO1_, '(3X,2A9)') 'vx(i)', 'vy(i)'
      DO i = 1, N
         WRITE(IO1_, '(3X,F9.4,'', '',F9.5)') vx(i), vy(i)
      END DO
      CLOSE(IO1_)
      END
