program dqmc2
*	diffusion quantum MC program
*	use trial wave function to select initial
*	configurations, and importance sampling
*	see P. J. Reynolds, J. Tobochnik, and H. Gould,
*	"Diffusion Quantum Monte Carlo," Computers in Physics 4, 662 (1990).
	real*8 conf(10000),E0,dt
	integer*4 m,nmcs,iseed
	call initial(m,nmcs,E0,dt,conf,iseed)
	call diffus(m,nmcs,E0,dt,conf,iseed)
	stop
	end

	subroutine initial(m,nmcs,E0,dt,conf,iseed)
*	choose parameters and generate initial configurations
	real*8 conf(10000),E0,dt
	integer*4 m,nmcs,iseed
	real*8 dx,x,xtry,p,psit2
	integer*4 nequil0,ncorrel,isep,imcs,i
	write(*,*) 'enter E0,nmcs,dt,iseed'
	read(*,*) E0,nmcs,dt,iseed
*	E0 = reference energy = approximation to ground state
*	desired number of configurations
	m = 100
*	equilibration for MC for initial configurations
	nequil0 = 20
	ncorrel = 10
	dx = 0.5e0
	x = 2.0e0*(ran(iseed) - 0.5e0) ! initial guess for x
*	equilibration
	do 10 imcs = 1,nequil0
	xtry = x + (ran(iseed) - 0.5e0)*dx
	p = psit2(xtry)/psit2(x)
	if (p .gt. ran(iseed)) x = xtry
*	MC test for saving configuration
10	continue
*	choose initial configuration every ncorrel moves
	do 20 i = 1,m
	   do 30 isep = 1,ncorrel
	      xtry = x + (ran(iseed) - 0.5e0)*dx
	      p = psit2(xtry)/psit2(x)
	      if (p .gt. ran(iseed)) x = xtry
30	   continue
	   conf(i) = x
20	continue
	end

	subroutine diffus(m,nmcs,E0,dt,conf,iseed)
	real*8 E0,dt,conf(100)
	integer*4 m,nmcs,iseed
	real*8 g(10000),v,ecum,f,el
	integer*4 i,m0,mnow,imcs,iw,ir,icopy
	m0 = m
	ecum = 0.0e0
	do 10 imcs = 1,nmcs
*	   adjust reference energy
	   if (mod(imcs,20) .eq. 0) then
	      E0 = E0 - 0.1e0*(m - m0)/(dt*m0)
	      ecum = ecum + E0
	      write(*,*) imcs,20.e0*ecum/imcs,m
	   end if
	   mnow = m   ! number of configurations at this time
	   i = 1      ! index of configuration
*	   produce mnow Gaussian random numbers
	   call gran(0.5e0,dt,g,mnow,iseed)
	   ir = 0     ! index for random numbers
50	   continue
	   ir = ir + 1
*	   change configuration from Gaussian distribution
	   conf(i) = conf(i) + g(ir) + 0.5*dt*f(conf(i))
*	   make copies if needed
	   iw = int(ran(iseed) +
     &          exp(dt*(-0.5e0*(el(conf(i)) + el(conf(i)-g(ir))) + E0)))
	   if (iw .eq. 0) then
*	      delete configuration
	      conf(i) = conf(mnow)
	      conf(mnow) = conf(m)
	      mnow = mnow - 1
	      m = m - 1
	      if (m .eq. 0) then
	         write(*,*) 'm,mnow = ',m,mnow
	         stop
	      end if
	   elseif (iw .eq. 1) then
	      i = i + 1
	   else
	      do 20 icopy = 1,iw - 1
*	         make copies of configuration
	         m = m + 1
	         if (m .gt. 10000) write(*,*) m,i,iw
	         conf(m) = conf(i)
20	      continue
	      i = i + 1
	   end if
	   if (i .le. mnow) goto 50
10	continue
	end

	real*8 function psit2(x)
*	unnormalized trial wavefunction squared
	real*8 x
	psit2 = exp(-2*x*x)
	end

	real*8 function f(x)
*	f = 2 d/dx ln(psiT)
	real*8 x
	f = -4*x
	end

	real*8 function el(x)
*	el = H*psit/psit, first term is V(x)
	real*8 x
	el = 0.5e0*x*x + 1 - 2*x*x
	end

	subroutine gran(D,dt,g,n,iseed)
	real*8 D,dt,g(100000)
	integer*4 iseed,n
	real*8 a,theta,tpi
	integer*4 i
	tpi = 6.283185307e0
	do 10 i = 1,n,2
	   a = sqrt(-4.0*D*dt*log(ran(iseed)))
	   theta = tpi*(ran(iseed))
	   g(i) = a*sin(theta)
	   g(i+1) = a*cos(theta)
10	continue
	end