      PROGRAM ideal
*     1次元の古典理想気体にたいするデーモン・アルゴリズム
      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                                ! 初期デーモン・エネルギー
      WRITE(*,'(A,$)') 'number of MC steps per particle = '
      READ(*,*) mcs
      WRITE(*,'(A,$)') 'maximum change in velocity = '
      READ(*,*) dvmax
*     エネルギーを各粒子に等しく分配する
      vinitial = sqrt(2*E/N)                ! 質量は 1
*     各粒子の初速度は全て同じ
      DO i = 1, N
         v(i) = vinitial
      END DO
*     和の初期化
      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          ! 速度の試行変化
         ipart = int(rnd(0)*N + 1)          ! ランダムに粒子を選択する
         vtrial = v(ipart) + dv             ! 試行速度
*        試行のエネルギー変化
         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
*     粒子当りのモンテカルロ・ステップごとにデータを積算する
      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                    ! 平均のデーモン・エネルギー
      norm = norm/N
      accept_prob = accept*norm             ! 受容確率
*     １粒子当りの平均値
      Esave = Ecum*norm                     ! １粒子当りの平均エネルギー
      WRITE(*,*) 'mean demon energy =', Edave
      WRITE(*,*) 'mean system energy per particle =', Esave
      WRITE(*,*) 'acceptance probability =', accept_prob
      END
