PROGRAM qmwalk
      IMPLICIT NONE
      REAL*8 Esum, Vave, Vref, binsize, ds, dt, psi(-500:500), x(2000)
      INTEGER IR_, N, N0, i, imcs, mcs, nequil
      DATA psi, x/3001*0.0D0/
      CALL GWopen(IR_, 0)
      CALL initial(x,N,N0,Vref,binsize,ds,dt,mcs,nequil,Esum)
      DO i = 1, nequil                      ! 平衡化
         CALL walk(x,N,N0,Vref,Vave,ds,dt)
      END DO
      DO imcs = 1, mcs
         CALL walk(x,N,N0,Vref,Vave,ds,dt)
         CALL mydata(x,psi,N,Vref,Vave,Esum,imcs,binsize)
      END DO
      CALL plot_psi(psi,binsize)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(x,N,N0,Vref,binsize,ds,dt,mcs,nequil,Esum)
      IMPLICIT NONE
      REAL*8 width, x(2000), Vref, binsize, ds, dt, Esum, V
      REAL dmy, rnd
      INTEGER IR_, i, N, N0, mcs, nequil
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'desired number of walkers = '
      READ(*,*) N0
      WRITE(*,'(A,$)') 'step length = '
      READ(*,*) ds
      dt = ds*ds                            ! 時間刻み
      WRITE(*,'(A,$)') 'number of Monte Carlo steps per walker = '
      READ(*,*) mcs
      nequil = int(0.4D0*mcs)
*     初期の粒子数は任意に選ぶ
      N = N0
*     粒子の初期位置はランダムに選ぶ
      width = 1                             ! 粒子の領域の初期幅
      Vref = 0
      DO i = 1, N
         x(i) = (2*rnd(0) - 1)*width
         Vref = Vref + V(x(i))
      END DO
      Vref = Vref/N
      binsize = 2*ds                        ! 配列 psi に対応する区分け用の箱の幅
      Esum = 0
      CALL GWclear(IR_, -1)
      WRITE(*,'(4A12)') 'MC steps', 'N', 'E', 'Vref'
      WRITE(*,*)
      END
C
      SUBROUTINE walk(x,N,N0,Vref,Vave,ds,dt)
      IMPLICIT NONE
      REAL*8 Nin, Vsum, dV, potential, x(2000), Vref, Vave, ds, dt, V
      REAL rnd
      INTEGER i, N, N0
      Nin = N                               ! 試行の初めにおける粒子数
      Vsum = 0
      DO i = Nin, 1, -1
         IF(rnd(0) .LT. 0.5D0) THEN
            x(i) = x(i) + ds
         ELSE
            x(i) = x(i) - ds
         END IF
         potential = V(x(i))                ! x でのポテンシャル
         dV = potential - Vref
         IF(dV .LT. 0) THEN
            IF(rnd(0) .LT. -dV*dt) THEN
               N = N + 1
               x(N) = x(i)                  ! 新しい粒子
*              係数 2 は x に 2 個の粒子があるため
               Vsum = Vsum + 2*potential
            ELSE
               Vsum = Vsum + potential      ! 古い粒子のみ
            END IF
         ELSE
            IF(rnd(0) .LT. dV*dt) THEN
               x(i) = x(N)
               N = N - 1
            ELSE
               Vsum = Vsum + potential
            END IF
         END IF
      END DO
      Vave = Vsum/N                         ! 平均ポテンシャル
      Vref = Vave - (N - N0)/(N0*dt)        ! 新しい参照エネルギー
      END
C
      SUBROUTINE mydata(x,psi,N,Vref,Vave,Esum,imcs,binsize)
      IMPLICIT NONE
      REAL*8 x(2000), psi(-500:500), Vref, Vave, Esum, binsize
      INTEGER bin, i, N, imcs
*     データを積算する
      Esum = Esum + Vave
*     粒子を'箱'に入れる
      DO i = 1, N
         bin = int(x(i)/binsize)
         psi(bin) = psi(bin) + 1
      END DO
      IF(mod(imcs,10) .EQ. 0) THEN
         WRITE(*, '(I12,I12,F12.3,F12.3)') 
     +            imcs, N, Esum/imcs, Vref
      END IF
      END
C
      SUBROUTINE plot_psi(psi,binsize)
      IMPLICIT NONE
      REAL*8 P, norm, pmax, sum, ymax, psi(-500:500), binsize
      INTEGER IR_, i, imax
      LOGICAL Ltmp1_
      i = 0
      pmax = psi(i)*psi(i)
      sum = pmax
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         i = i + 1
         P = psi(i)*psi(i) + psi(-i)*psi(-i)
         IF(P .GT. pmax) pmax = P
         sum = sum + P
         Ltmp1_ = P .NE. 0
      END DO
      imax = i
      norm = sqrt(sum*binsize)
      ymax = sqrt(pmax)/norm
      CALL GWindow(IR_, REAL(-imax), REAL(-0.1D0*ymax), REAL(imax)
     +, REAL(1.1D0*ymax))
      CALL GWline(IR_, REAL(-imax), 0.0, REAL(imax), 0.0)
      DO i = -imax, imax
         CALL GWline(IR_, REAL(i-1), REAL(psi(i-1)/norm), REAL(i)
     +   , REAL(psi(i)/norm))
      END DO
      END
C
      FUNCTION V(x)
      IMPLICIT NONE
      REAL*8 V, x
      V = 0.5D0*x*x
      END