      PROGRAM boltzmann
      IMPLICIT NONE
      REAL*8 E, P(-100:100), accept, accum(3), beta, delta, delv, v
      INTEGER imcs, mcs, nequil, nbin
*     Metropolis algorithm for a particle in one dimension
      CALL initial(v,E,beta,mcs,nequil,delta,nbin,delv)
      DO imcs = 1, nequil
         CALL metropolis(v,E,beta,delta,accept)
      END DO
      CALL initialize_sums(P,accum,accept,nbin)
      DO imcs = 1, mcs
         CALL metropolis(v,E,beta,delta,accept)
*        accumulate data after each trial change
         CALL mydata(P,accum,v,E,delv)
      END DO
      CALL output(P,accum,mcs,accept,nbin,delv)
      END
C
      SUBROUTINE initial(v0,E0,beta,mcs,nequil,delta,nbin,delv)
      IMPLICIT NONE
      REAL*8 T, vmax, v0, E0, beta, delta, delv
      REAL dmy, rnd
      INTEGER mcs, nequil, nbin
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'temperature = '
      READ(*,*) T
      beta = 1/T
      WRITE(*,'(A,$)') 'number of Monte Carlo steps = '
      READ(*,*) mcs
      nequil = int(0.1D0*mcs)
      WRITE(*,'(A,$)') 'initial speed = '
      READ(*,*) v0
      E0 = 0.5D0*v0*v0                      ! initial kinetic energy
      WRITE(*,'(A,$)') 'maximum change in velocity = '
      READ(*,*) delta
      vmax = 10*sqrt(T)
      nbin = 20                             ! number of bins
      delv = vmax/nbin                      ! velocity interval
      END
C
      SUBROUTINE initialize_sums(P,accum,accept,nbin)
      IMPLICIT NONE
      REAL*8 P(-100:100), accum(3), accept
      INTEGER i, ibin, nbin
      DO ibin = -nbin, nbin
         P(ibin) = 0
      END DO
      DO i = 1, 3
         accum(i) = 0
      END DO
      accept = 0
      END
C
      SUBROUTINE metropolis(v,E,beta,delta,accept)
      IMPLICIT NONE
      REAL*8 dE, dv, vtrial, v, E, beta, delta, accept
      REAL rnd
      dv = (2*rnd(0) - 1)*delta             ! trial velocity change
      vtrial = v + dv                       ! trial velocity
      dE = 0.5D0*(vtrial*vtrial - v*v)      ! trial energy change 
      IF(dE .GT. 0) THEN
         IF(exp(-beta*dE) .LT. rnd(0)) THEN
            RETURN
         END IF
      END IF
      v = vtrial
      accept = accept + 1
      E = E + dE
      END
C
      SUBROUTINE mydata(P,accum,v,E,delv)
      IMPLICIT NONE
      REAL*8 P(-100:100), accum(3), v, E, delv
      INTEGER ibin
      accum(1) = accum(1) + E
      accum(2) = accum(2) + E*E
      accum(3) = accum(3) + v
      ibin = int(v/delv)
      P(ibin) = P(ibin) + 1
      END
C
      SUBROUTINE output(P,accum,mcs,accept,nbin,delv)
      IMPLICIT NONE
      REAL*8 E2ave, Eave, prob, sigma2, v, vave, P(-100:100),
     +       accum(3), accept, delv
      INTEGER mcs, nbin, ibin
      accept = accept/mcs
      WRITE(*,*) 'acceptance probability =', accept
      vave = accum(3)/mcs
      WRITE(*,*) 'mean velocity =', vave
      Eave = accum(1)/mcs
      WRITE(*,*) 'mean energy =', Eave
      E2ave = accum(2)/mcs
      sigma2 = E2ave - Eave*Eave
      WRITE(*,*) 'sigma_E = ', sqrt(sigma2)
      WRITE(*,*)
      WRITE(*,'(2A13)') ' v ', 'P(v)'
      WRITE(*,*)
      v = -nbin*delv
      DO ibin = -nbin, nbin
         IF(p(ibin) .GT. 0) THEN
            prob = p(ibin)/mcs
            WRITE(*,'(2F13.6)') v, prob
         END IF
         v = v + delv
      END DO
      END
