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, ii
      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_/2/
      DATA IC2_/1/
      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.
c     DO WHILE(Ltmp1_)
      do ii = 1, 2000
         CALL show_positions(flag,IC2_)
         CALL Verlet(t,ke,pe,virial)
         CALL show_output(t,ke,pe,virial,kecum,vcum,ncum,area,IC1_)
c        call check_stop (ncum,flag)
         Ltmp1_ = .NOT.(flag .EQ. 'stop')
      END DO
      CALL save_config
c     CALL GWquit(IR_)
      call pgclos
      END

      subroutine check_stop (nucm, flag)
      integer ncum
      character*80 flag
      if ( mod(ncum,1000) .eq. 0 ) then
        open (unit=21, file='md_stop.txt', status='old',err= 201)
        read (*,*) flag
        close (21)
      endif
 201  continue
      end

C mhl done
      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

c mhl added
      do i=1,n
        xsave(i) = x(i)
        ysave(i) = y(i)
      end do

      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
c     initialize sums
      kecum = 0
      pecum = 0
      vcum = 0
      END

C mhl done
      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_, pgopen
      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

      if( pgopen('/xwin') .le. 0) stop 'Open device fail.'

      call pgpap (10.0, 0.5)
      call pgiden
  
      call pgstbg(0)

      call pgsubp (2,1)

      call pgenv (0.0, real(Lx), 0.0, real(Ly), 0, -1)
      call pgsch (1.5)
      call pglab ('', '', 'Program for MD siulation.')
      call pgpanl (2,1)
      call headings

      END


C mhl done
      SUBROUTINE headings
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      INTEGER IR_

c     call pgsch (2.0)
      call pgsch (1.5)

      WRITE(Stmp1_, '(A)') 'time steps '
      call pgptxt(2.5, 4.5, 0.0, 1.0, 'time steps = ')
      WRITE(Stmp1_, '(A)') 'time '
      call pgptxt(2.5, 3.5, 0.0, 1.0, 'time = ')
      WRITE(Stmp1_, '(A)') 'energy '
      call pgptxt(2.5, 2.5, 0.0, 1.0, 'energy = ')
      WRITE(Stmp1_, '(A)') '<T> '
      call pgptxt(2.5, 1.5, 0.0, 1.0, '<T> = ')
      WRITE(Stmp1_, '(A)') '<P>'
      call pgptxt(2.5, 0.5, 0.0, 1.0, '<P> = ')
      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
c        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
c        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)
c           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
      real radius
      parameter (radius = 0.10)
      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 pgpanl (1,1)
        call pgsci (1)

c       call pgbbuf
        do i = 1, n
          call pgsci (0)
          call pgcirc (real(xsave(i)), real(ysave(i)), radius)
          call pgsci (1)
          call pgbox ('bc',0.0, 0,'bc',0.0 ,0)
          call pgsci (mod(i,4)+1)
          call pgcirc ( real(x(i)), real(y(i)), radius )
        end do
c       call pgebuf

        do i = 1, n
          xsave(i) = x(i)
          ysave(i) = y(i)
        end do

      END

C mhl done
      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 pgpanl (2,1)
c     call pgsch (2.0)
      call pgsch (1.5)

      ncum = ncum + 1

c     call pgbbuf

      WRITE(Stmp1_, '(I10)') ncum
c     call pgtext (1.5,4.5,'time steps = '//Stmp1_)
      call pgptxt (2.5,4.5, 0.0, 0.0, Stmp1_)
      WRITE(Stmp1_, '(F6.2)') t             ! time
c     call pgtext (1.5,3.5,'time = '//Stmp1_)
      call pgptxt (2.5,3.5, 0.0, 0.0, Stmp1_)
      E = ke + pe                           ! total energy
      WRITE(Stmp1_, '(F8.3)') E
c     call pgtext (1.5,2.5,'<E> ='//Stmp1_)
      call pgptxt (2.5,2.5, 0.0, 0.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
c     call pgtext (1.5,1.5,'<T> ='//Stmp1_)
      call pgptxt (2.5,1.5, 0.0, 0.0, Stmp1_)
      WRITE(Stmp1_, '(F5.2,$)') p
c     call pgtext (1.5,0.5,'Pressure ='//Stmp1_)
      call pgptxt (2.5,0.5, 0.0, 0.0, Stmp1_)
      
c     call pgebuf

      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 filename
      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/
      filename = 'md_save.dat'
      OPEN(IO1_, file = filename)
      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