PROGRAM path_integral
DIM x(-1 to 501),psi2(-200 to 200),psi2t(-200 to 200)
DIM dEP(-1 to 501),dEL(-1 to 501),dER(-1 to 501)
OPEN #1: name "qmcdata3",access output,create new
OPEN #2: name "psi^2,3",access output,create new
CALL parameters(c1,c2,N,dxmax,x,nmcs,nout,nequil,dEP,dEL,dER)
FOR iout = 1 to nout+ nequil
    CALL mc(c1,c2,N,dxmax,x,psi2,nmcs,accept,dEP,dEL,dER)
    CALL average(c1,c2,N,dxmax,accept,nmcs,iout,nout,nequil,psi2,psi2t,E,E2,#1,#2)
NEXT iout
CLOSE #1
END

SUB parameters(c1,c2,N,dxmax,x(),nmcs,nout,nequil,dEP(),dEL(),dER())
    RANDOMIZE
    LET ep = 0.15                 ! total time = inverse temp = N*ep
    LET N = 100
    LET dxmax = 1.4               ! twice maximum move
    LET c1 = 1/(2*ep)
    LET c2 = ep/2
    LET nmcs = 50000      ! total steps
    LET nout = 50
    LET nequil = 20
    FOR i = 1 to N
        LET x(i) = (rnd - 0.5)
        IF abs[x(i) - x(i-1)] > 0.2 then
           LET x(i) = (x(i) + x(i-1))/2
        END IF
        LET dEL(i) = [(x(i) - x(i-1))^2]*c1
        LET dER(i) = [(x(i) - x(i+1))^2]*c1
        LET dEP(i) = [exp(-2*x(i)) - 2*exp(-x(i))]*c2
    NEXT i
    LET x(N+1) = x(1)             ! pbc
    LET x(0) = x(N)
    PRINT "ep,N,nmcs,nequil = ",ep,N,nmcs,nequil
END SUB

SUB mc(c1,c2,N,dxmax,x(),psi2(),nmcs,accept,dEP(),dEL(),dER())
    FOR imcs = 1 to nmcs
        LET i = 1 + int(rnd*N)
        LET xtry = x(i) + (rnd - 0.5)*dxmax
        LET dERT = [(xtry - x(i+1))^2]*c1
        LET dELT = [(xtry - x(i-1))^2]*c1
        LET dEPT = [exp(-2*xtry) - 2*exp(-xtry)]*c2
        LET dE = dERT + dELT - dER(i) - dEL(i) + dEPT - dEP(i)
        IF exp(-de) > rnd then
           LET x(i) = xtry
           LET dEP(i) = dEPT
           LET dEL(i) = dELT
           LET dER(i) = dERT
           LET dER(i-1) = dELT
           LET dEL(i+1) = dERT
           IF (i = 1) then         ! pbc
              LET x(N+1) = xtry
           ELSE IF (i = N) then
              LET x(0) = xtry
           END IF
           LET accept = accept + 1
        END IF
        LET ibin = 10*x(i)  ! 10 bins 
        LET psi2(ibin) = psi2(ibin) + 1     ! accumulate probability
    NEXT imcs
END SUB

SUB AVERAGE(c1,c2,N,dxmax,accept,nmcs,iout,nout,nequil,psi2(),psi2t(),E,E2,#1,#2)
    LET accept = accept/(nmcs)
    LET cexact = 1/sqr(pi)
    LET ener = 0
    FOR ibin = -200 to 200
        IF psi2(ibin) > 0 then
           LET psi2(ibin) = psi2(ibin)/nmcs
           LET x = ibin*0.1
           LET ener = ener + psi2(ibin)*(exp(-2*x)-2*exp(-x))/2
           LET ener = ener + x*(-2*exp(-2*x) + 2*exp(-x))/4  ! V(x) + .5*x*dV/dx
        END IF
    NEXT ibin
    PRINT #1: "energy = " , ener
    PRINT #1: "acceptance ratio = " , accept
    PRINT "energy = " , ener
    PRINT "acceptance ratio =", accept
    LET accept = 0
    IF iout > nequil then
       LET e = e + ener
       LET e2 = e2 + ener*ener
       LET eave = e/(iout - nequil)
       LET e2ave = e2/(iout - nequil)
       LET sig = sqr[abs((e2ave - eave*eave)/(iout-nequil))]
       PRINT "eave,e2ave,sig = ", eave,e2ave,sig
       PRINT #1: "eave,e2ave,sig = ", eave,e2ave,sig
       FOR ibin = -200 to 200
           LET psi2t(ibin) = psi2t(ibin) + psi2(ibin)
           LET psi2(ibin) = 0
       NEXT ibin
    END IF
    IF iout = nout + nequil then
       FOR ibin = -200 to 200
           IF psi2t(ibin) > 0 then
              LET psi2t(ibin) = psi2t(ibin)*10/nout
              PRINT #2: str$(ibin/10) & chr$(9) & str$(psi2t(ibin))
           END IF
       NEXT ibin
       FOR ibin = -200 to 200
           IF psi2t(ibin) > 0 then
              LET x = ibin*0.1
              LET d2 = sqr(psi2t(ibin+1)) + sqr(psi2t(ibin-1))
              LET d2 = d2 - 2*sqr(psi2t(ibin))*100
              LET pd2p = d2*sqr(psi2t(ibin))
              LET V = [exp(-2*x) - 2*exp(-x)]/2
              LET pVp = psi2t(ibin)*V
              LET pHp = -pd2p/2 + pVp
              LET enert = enert + pHp/10    ! using running average
              LET enertv = enertv + psi2t(ibin)*(exp(-2*x)-2*exp(-x))/2   ! virial
              LET enertv = enertv + (x*(-2*exp(-2*x) + 2*exp(-x)))/4/10  ! V(x) + .5*x*dV/dx
           END IF
       NEXT ibin
       PRINT #1: "enert,enertv = ",enert,enertv
       PRINT "enert,enertv = ",enert,enertv
    END IF
END SUB