program dqmc1
*	diffusion quantum MC without importance sampling
*	see P. J. Reynolds, J. Tobochnik, and H. Gould,
*	"Diffusion Quantum Monte Carlo," Computers in Physics 4, 662 (1990).
	real*8 E0,dt,conf(10000)
	integer*4 m,nmcs,iseed
	call initial(m,nmcs,E0,dt,conf,iseed)
	call mc(m,nmcs,E0,dt,conf,iseed)
	end

	subroutine initial(m,nmcs,E0,dt,conf,iseed)
	integer*4 m,nmcs,i,iseed
	real*8 dx,x,xtry,p,conf(10000),E0,dt
	write(*,*) 'enter E0,nmcs,dt,iseed'
	read(*,*) E0,nmcs,dt,iseed
	m = 100
	do 30 i = 1,m
	   conf(i) = 2.0d0*(ran(iseed) - 0.5)
30	continue
	end

	subroutine mc(m,nmcs,E0,dt,conf,iseed)
	integer*4 iseed,i,m,m0,mnow,nmcs,iw,icopy,ir
	real*8 conf(10000),dt,g(1000),E0,v
	real*8 ecum
	integer*4 imcs
	m0 = m
	ecum = 0.0d0
	do 10 imcs = 1,nmcs
*	   adjust reference energy
	   if (mod(imcs,20) .eq. 0) then
	      E0 = E0 - 0.1*(m - m0)/(dt*m0)
	      ecum = ecum + E0
	      write(*,*) imcs,20.d0*ecum/dfloat(imcs),m
	   end if
	   mnow = m
	   i = 1
*	   generate mnow random numbers distributed according to gaussian
	   call gran(0.5d0,dt,g,mnow,iseed)
	   ir = 0
20	   continue
*	   move coordinates of a configuration randomly
	   ir = ir + 1
	   conf(i) = conf(i) + g(ir)
*	   make copies if needed
	   iw = int(ran(iseed)
     c     + exp(-0.5d0*dt*(v(conf(i)) + v(conf(i) - g(ir)) - 2*E0)))
	   if (iw .eq. 0) then
	      conf(i) = conf(mnow)
	      conf(mnow) = conf(m)
	      mnow = mnow - 1
	      m = m - 1
	   elseif (iw .eq. 1) then
	      i = i + 1
	   else
	      do 30 icopy = 1,iw-1
	         m = m + 1
	         conf(m) = conf(i)
30	      continue
	      i = i + 1
	   end if
	   if (i .le. mnow) goto 20
10	continue
	end

	real*8 function v(x)
	real*8 x
	v = 0.5d0*x*x
	end
 
	subroutine gran(D,dt,g,n,iseed)
*	produces n random numbers distributed according to a gaussian
	real*8 g(1000),D,dt
	integer*4 iseed,n
	real*8 a,theta
	integer*4 ir
	do 10 ir = 1,n,2
	   a = sqrt(-4.0*D*dt*alog(ran(iseed)))
	   theta = 6.2831853D0*(ran(iseed))
	   g(ir) = a*sin(theta)
	   g(ir+1) = a*cos(theta)
10	continue
	end