PROGRAM qmwalk
DIM x(2000),psi(-500 to 500)
CALL initial(x(),N,N0,Vref,binsize,ds,dt,mcs,nequil,Esum)
FOR i = 1 to nequil          ! equilibration
    CALL walk(x(),N,N0,Vref,Vave,ds,dt)
NEXT i
FOR imcs = 1 to mcs
    CALL walk(x(),N,N0,Vref,Vave,ds,dt)
    CALL data(x(),psi(),N,Vref,Vave,Esum,imcs,binsize)
NEXT imcs
CALL plot_psi(psi(),binsize)
END

SUB initial(x(),N,N0,Vref,binsize,ds,dt,mcs,nequil,Esum)
    DECLARE DEF V
    RANDOMIZE
    INPUT prompt "desired number of walkers = ": N0
    INPUT prompt "step length = ": ds
    LET dt = ds*ds           ! time step
    INPUT prompt "number of Monte Carlo steps per walker = ": mcs
    LET nequil = round(0.4*mcs)
    ! choose initial number of walkers equal to desired number
    LET N = N0
    ! choose initial positions of walkers at random
    LET width = 1            ! initial width of region for walkers
    LET Vref = 0
    FOR i = 1 to N
        LET x(i) = (2*rnd - 1)*width
        LET Vref = Vref + V(x(i))
    NEXT i
    LET Vref = Vref/N
    LET binsize = 2*ds       ! binsize for psi array
    LET Esum = 0
    CLEAR
    PRINT "MC steps","N","E","Vref"
    PRINT
END SUB

SUB walk(x(),N,N0,Vref,Vave,ds,dt)
    DECLARE DEF V
    LET Nin = N              ! # walkers at beginning of trial
    LET Vsum = 0
    FOR i = Nin to 1 step -1
        IF rnd < 0.5 then
           LET x(i) = x(i) + ds
        ELSE
           LET x(i) = x(i) - ds
        END IF
        LET potential = V(x(i))   ! potential at x
        LET dV = potential - Vref
        IF dV < 0 then       ! check to add walker
           IF rnd < -dV*dt then
              LET N = N + 1
              LET x(N) = x(i)     ! new walker
              ! factor of 2 since two walkers at x
              LET Vsum = Vsum + 2*potential
           ELSE
              LET Vsum = Vsum + potential   ! only do old walker
           END IF
        ELSE
           IF rnd < dV*dt then    ! check to remove walker
              LET x(i) = x(N)
              LET N = N - 1
           ELSE
              LET Vsum = Vsum + potential
           END IF
        END IF
    NEXT i
    LET Vave = Vsum/N        ! mean potential
    LET Vref = Vave - (N - N0)/(N0*dt)      ! new reference energy
END SUB

SUB data(x(),psi(),N,Vref,Vave,Esum,imcs,binsize)
    ! accumulate data
    LET Esum = Esum + Vave
    ! bin walkers
    FOR i = 1 to N
        LET bin = round(x(i)/binsize)
        LET psi(bin) = psi(bin) + 1
    NEXT i
    IF mod(imcs,10) = 0 then
       PRINT imcs,
       PRINT using "####": N;
       PRINT,
       PRINT using "--.###": Esum/imcs;
       PRINT,
       PRINT using "--.###": Vref
    END IF
END SUB

SUB plot_psi(psi(),binsize)
    PRINT "Hit any key to see |psi|"
    DO
    LOOP until key input
    GET KEY k
    CLEAR
    LET i = 0
    LET pmax = psi(i)*psi(i)
    LET sum = pmax
    DO
       LET i = i + 1
       LET P = psi(i)*psi(i) + psi(-i)*psi(-i)
       IF P > pmax then LET pmax = P
       LET sum = sum + P
    LOOP until P = 0
    LET imax = i
    LET norm = sqr(sum*binsize)
    LET ymax = sqr(pmax)/norm
    SET WINDOW -imax,imax,-0.1*ymax,1.1*ymax
    PLOT LINES: -imax,0;imax,0
    FOR i = -imax to imax
        PLOT LINES: i-1,psi(i-1)/norm;i,psi(i)/norm
    NEXT i
END SUB

DEF V(x)
    LET V = 0.5*x*x
END DEF