PROGRAM ideal
*     demon algorithm for the one-dimensional, ideal classical gas 
      IMPLICIT NONE
      REAL*8 E, Ecum, Ed, Edcum, dvmax, v(100)
      INTEGER N, accept, imcs, mcs
      CALL initial(N,v,mcs,E,Ed,Ecum,Edcum,accept,dvmax)
      DO imcs = 1, mcs
         CALL change(N,v,E,Ed,Ecum,Edcum,accept,dvmax)
      END DO
      CALL averages(N,Ecum,Edcum,mcs,accept)
      END
C
      SUBROUTINE initial(N,v,mcs,E,Ed,Ecum,Edcum,accept,dvmax)
      IMPLICIT NONE
      REAL*8 vinitial, v(100), E, Ed, Ecum, Edcum, dvmax
      REAL dmy, rnd
      INTEGER i, N, mcs, accept
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'number of particles = '
      READ(*,*) N
      WRITE(*,'(A,$)') 'initial energy of system = '
      READ(*,*) E
      Ed = 0                                ! initial demon energy
      WRITE(*,'(A,$)') 'number of MC steps per particle = '
      READ(*,*) mcs
      WRITE(*,'(A,$)') 'maximum change in velocity = '
      READ(*,*) dvmax
*     divide energy equally among particles
      vinitial = sqrt(2*E/N)                ! mass unity
*     all particles have same initial velocities
      DO i = 1, N
         v(i) = vinitial
      END DO
*     initialize sums
      Ecum = 0
      Edcum = 0
      accept = 0
      END
C
      SUBROUTINE change(N,v,E,Ed,Ecum,Edcum,accept,dvmax)
      IMPLICIT NONE
      REAL*8 de, dv, vtrial, v(100), E, Ed, Ecum, Edcum, dvmax
      REAL rnd
      INTEGER i, ipart, N, accept
      DO i = 1, N
         dv = (2*rnd(0) - 1)*dvmax          ! trial change in velocity
         ipart = int(rnd(0)*N + 1)          ! select random particle
         vtrial = v(ipart) + dv             ! trial velocity
*        trial energy change
         de = 0.5D0*(vtrial*vtrial - v(ipart)*v(ipart))
         IF(de .LE. Ed) THEN
            v(ipart) = vtrial
            accept = accept + 1
            Ed = Ed - de
            E = E + de
         END IF
      END DO
*     accumulate data after each Monte Carlo step per particle
      Ecum = Ecum + E
      Edcum = Edcum + Ed
      END
C
      SUBROUTINE averages(N,Ecum,Edcum,mcs,accept)
      IMPLICIT NONE
      REAL*8 Edave, Esave, accept_prob, norm, Ecum, Edcum
      INTEGER N, mcs, accept
      norm = 1.0D0/mcs
      Edave = Edcum*norm                    ! mean demon energy
      norm = norm/N
      accept_prob = accept*norm             ! acceptance probability
*     system averages per particle
      Esave = Ecum*norm                     ! mean energy per system particle
      WRITE(*,*) 'mean demon energy =', Edave
      WRITE(*,*) 'mean system energy per particle =', Esave
      WRITE(*,*) 'acceptance probability =', accept_prob
      END