PROGRAM dqmc2
! Diffusion quantum MC, use trial wave function to select initial
! configurations, and then use importance sampling.
! one-dimensional harmonic oscillator
! see P. J. Reynolds, J. Tobochnik, and H. Gould,
! "Diffusion Quantum Monte Carlo," Computers in Physics 4, 662 (1990).
DIM conf(10000)
CALL initial(m,nmcs,E0,dt,conf())
CALL diffuse(m,nmcs,E0,dt,conf())
END

SUB initial(m,nmcs,E0,dt,conf())
   ! choose parameters and produce initial configurations
   RANDOMIZE
   DECLARE DEF psit2             ! trial wave function squared
   ! E0 = reference energy = approximation to ground state
   INPUT prompt "guess for ground state energy E0 = ? ": E0
   INPUT prompt "time interval dt = ? ": dt
   INPUT prompt "number of Monte Carlo steps nmcs = ? ": nmcs
   LET m = 100                   ! desired number of configurations
   LET nequil0 = 20              ! equilibration time for initial configuration
   LET ncorrel = 10              ! # MC steps between initial configurations
   LET dx = 0.5
   LET x = 2.0*rnd - 1           ! initial guess for x
   FOR imcs = 1 to nequil0       ! make initial choice of x "more" random
      LET xtry = x + (rnd - 0.5)*dx
      LET p = psit2(xtry)/psit2(x)
      IF rnd <= p then LET x = xtry  ! test for new configuration
   NEXT imcs
   ! choose initial configurations every ncorrel moves
   FOR i = 1 to m
      FOR isep = 1 to ncorrel
         LET xtry = x + (rnd - 0.5)*dx
         LET p = psit2(xtry)/psit2(x)
         IF rnd <= p then LET x = xtry
      NEXT isep
      LET conf(i) = x
   NEXT i
END SUB

SUB diffuse(m,nmcs,E0,dt,conf())
   DECLARE DEF el                ! el local energy
   DECLARE DEF f
   DIM g(10000)                  ! put random numbers into an array
   LET m0 = m
   LET ecum = 0.0
   PRINT "MC steps", "mean ground state energy", "# configurations"
   FOR imcs = 1 to nmcs
      IF mod(imcs,20) = 0 then  ! adjust reference energy
         LET E0 = E0 - 0.1*(m - m0)/(dt*m0)
         LET ecum = ecum + E0
         PRINT imcs,20*ecum/imcs,,m
      END IF
      LET mnow = m              ! number of configurations at this time
      LET i = 1                 ! number of configuration
      CALL gauss(0.5,dt,g(),mnow)    ! produce mnow Gaussian random numbers
      LET iran = 0              ! index for random numbers
      DO
         LET iran = iran + 1
         LET x = conf(i)
         LET conf(i) = x + g(iran)   ! new configuration from Gaussian distribution
         LET conf(i) = conf(i) + 0.5*dt*f(x)   ! add "drift" term
         ! make copies if needed
         LET iw = rnd + exp(dt*(-0.5*(el(conf(i)) + el(conf(i) - g(iran))) + E0))
	 LET iw = int(iw)
         IF iw = 0 then
            ! delete configuration
            LET conf(i) = conf(mnow)
            LET conf(mnow) = conf(m)
            LET mnow = mnow - 1
            LET m = m - 1
            IF m = 0 then
               PRINT "m,mnow = ";m,mnow
               STOP
            END IF
         ELSEIF iw = 1 then
            LET i = i + 1
         ELSE
            FOR icopy = 1 to iw - 1  ! make copies of configuration
               LET m = m + 1
               LET conf(m) = conf(i)
            NEXT icopy
            LET i = i + 1
         END IF
      LOOP until i > mnow
   NEXT imcs
END SUB

DEF psit2(x)
   ! unnormalized trial wavefunction squared
   LET psit2 = exp(-2*x*x)
END DEF

DEF f(x)
   ! F = 2 d/dx ln(psiT), see eq. (9)
   LET f = -4*x
END DEF

DEF el(x)
   ! EL = H*psiT/psiT, first term is v(x)
   LET el = 0.5*x*x + 1 - 2*x*x
END DEF

SUB gauss(D,dt,g(),n)
   LET tpi = 6.283185307
   FOR i = 1 to n step 2
      LET a = sqr(-4.0*D*dt*log(rnd))
      LET theta = tpi*(rnd)
      LET g(i) = a*sin(theta)
      LET g(i+1) = a*cos(theta)
   NEXT i
END SUB