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