c....... simple DMC code for fermions
c........(c) L. Mitas, July 1998, NATO ASI
c....... does two electron systems, eg, He,H-,Li+ etc
c....... ulec() is a uniform random number generator (0,1) 
c....... He triplet -2.175229
c....... He sp triplet -2.133
c
c.......................................initialize etc
	implicit double precision (a-h,o-z)
	parameter (mwalk=2000,mpar=2,mdim=3,mblm=1000)
	dimension r(mwalk,mpar,mdim),rnew(mwalk,mpar,mdim),
     &  	  dr(mwalk,mpar,mdim),eloc(mwalk),elocn(mwalk),
     &            eblock(mblm),eblock2(mblm),eloch(mwalk)
     &           ,drn(mwalk,mpar,mdim),wf(mwalk),wfn(mwalk),
     &           phi1(mpar),phi1d(mpar),phi1dd(mpar),
     &           phi2(mpar),phi2d(mpar),phi2dd(mpar) 
     &           ,drm(mpar)
c
	read (1,*) z,nwmean,esteig,tau,mblock,mstep,iskip
c
	iswit=0
	twopi=2.*3.14159265358979
	al=z
	al2=al/2.
	alal=al*al
	nwalk=nwmean
	rspan=3.
	eestim=esteig
	et=eestim
c
	write (*,1001) z,tau,nwmean
	write (*,1002) mblock,mstep,iskip
        write (2,1001) z,tau,nwmean
        write (2,1002) mblock,mstep,iskip
c
c.......................................generate walkers
c
	do id=1,mdim
	do ip=1,mpar
	do m=1,nwalk
	r(m,ip,id)=rspan*(2.*ulec()-1.)
c	if (m.eq.1) write (*,*) r(m,ip,id)
	enddo
	enddo
	enddo
c
c.......................................global loop over time evolution
c					blocks and inside block steps
c
	do ib=1,mblock
	eblock(ib)=0.
	eblock2(ib)=0.
	do is=1,mstep
c
c....... f200=(2-Zr)e(-Zr/2)
c
c.........................................move=diffuse and drift		
c
		do id=1,mdim
		do ip=1,mpar
		do m=1,nwalk
		gauss=dsqrt(-2.*dlog(ulec()))*dcos(twopi*ulec())
		chi=dsqrt(tau)*gauss
 		rnew(m,ip,id)=r(m,ip,id)+chi+tau*dr(m,ip,id)
		enddo
		enddo
		enddo	
c........................................calculate wf, drift, kinetic
c					 potential and local energy

		icross=0
		do m=1,nwalk
		r12=dsqrt((rnew(m,1,1)-rnew(m,2,1))**2+
     &		         (rnew(m,1,2)-rnew(m,2,2))**2+
     &			 (rnew(m,1,3)-rnew(m,2,3))**2)
		epot=1./r12
c........................................calc orbitals and derivatives
		   do ip=1,mpar
		   r1=dsqrt(rnew(m,ip,1)**2+rnew(m,ip,2)**2+
     &                rnew(m,ip,3)**2)
		   epot=epot-z/r1
                   phi1(ip)=exp(-al*r1)
                   phi1d(ip)=-al*phi1(ip)/r1
                   phi1dd(ip)=alal*phi1(ip)+2.*phi1d(ip)
c
                   ex1=exp(-al2*r1)
                   phi2(ip)=(2.-al*r1)*ex1
                   phi2d(ip)=-al*(2-al2*r1)*ex1/r1
		   phi2dd(ip)=al*al2*(3.-al2*r1)*ex1+2.*phi2d(ip)
		   enddo
c
		wfn(m)=phi1(1)*phi2(2)-phi1(2)*phi2(1)
		ekin=-.5*(phi1dd(1)*phi2(2)-phi1(2)*phi2dd(1)
     &                   +phi1(1)*phi2dd(2)-phi1dd(2)*phi2(1))/wfn(m)
		elocn(m)=ekin+epot
c........................................preventive cut-off of drift
c                                        (not crucial)
		drm(1)=0.
		drm(2)=0.
		do id=1,mdim
		drn(m,1,id)=rnew(m,1,id)*(phi1d(1)*phi2(2)-
     &                      phi1(2)*phi2d(1))
     &                      /wfn(m)
		drm(1)=drm(1)+drn(m,1,id)**2
		drn(m,2,id)=rnew(m,2,id)*(phi1(1)*phi2d(2)-
     &                      phi1d(2)*phi2(1))
     &                      /wfn(m)
		drm(2)=drm(2)+drn(m,2,id)**2
		enddo
c
		do ip=1,mpar
		if (drm(ip).gt.1./tau)		then
		do id=1,mdim
		drn(m,ip,id)=drn(m,ip,id)/sqrt(drm(ip))
		enddo
						endif
		enddo
c.........................................check crossing the node
c
		if (wfn(m)*wf(m).lt.0.)		then
		icross=icross+1
		elocn(m)=eloc(m)
		wfn(m)=wf(m)
		do ip=1,mpar
		do id=1,mdim
		drn(m,ip,id)=dr(m,ip,id)
		rnew(m,ip,id)=r(m,ip,id)
		enddo
		enddo
						endif
		enddo
c
c.........................................branch and update the list
c
		mn=0
		taueff=tau*float(nwalk-icross)/float(nwalk)
		do m=1,nwalk
		p=dexp(-taueff*(eloc(m)+elocn(m)-2.*et)/2.)
		eta=ulec()
		mnew=p+eta
			do m1=1,mnew
			mn=mn+1
			if (mn.gt.mwalk)		then
			write (*,*) ' population too large !!!'
							stop
							endif
			eloch(mn)=elocn(m)
			wf(mn)=wfn(m)
				do ip=1,mpar
				do id=1,mdim
				r(mn,ip,id)=rnew(m,ip,id)
				dr(mn,ip,id)=drn(m,ip,id)
				enddo
				enddo	
			enddo
		enddo
		nwalk=mn
		ampl=.01/tau
		etupd=-dlog(dfloat(nwalk)/dfloat(nwmean))
		et=eestim+ampl*etupd
c	
c.........................................accumulate energies for averages
c					  and dispersions
c
		rnorm=1./dfloat(nwalk)
		ecurr=0.
		e2curr=0.
		do m=1,nwalk
		eloc(m)=eloch(m)
		ecurr=ecurr+rnorm*(eloc(m)-esteig)
		e2curr=e2curr+rnorm*((eloc(m)-esteig)**2)
		enddo
		ecurr=ecurr+esteig
c
		eblock(ib)=eblock(ib)+(ecurr-esteig)/dfloat(mstep)
		eblock2(ib)=eblock2(ib)+e2curr/dfloat(mstep)
c
	enddo
		eblock(ib)=eblock(ib)+esteig
c......................................................end of inside the block
c......................................................write something out
c
		eacc=0.
		disacc=0.
		error=0.
		if(ib.gt.iskip) then
			do jb=iskip+1,ib
			eacc=eacc+eblock(jb)/dfloat(ib-iskip)
			disacc=disacc+eblock2(jb)/dfloat(ib-iskip)
			enddo
				if(ib.gt.iskip+1)	then
				   do jb=iskip+1,ib
				   error=error+(eblock(jb)-eacc)**2
				   enddo
			error=dsqrt(error/dfloat((ib-iskip)*(ib-iskip-1)))
							endif
				endif
c
	if (ib.gt.iskip)	then
				eestim=eacc
				else
				eestim=(eblock(ib)+eestim)/2.
				endif
	write (*,1000) ib,eblock(ib),eacc,error,sqrt(disacc),nwalk,et
	write (2,1000) ib,eblock(ib),eacc,error,sqrt(disacc),nwalk,et
1000	format(1x,'block=',i2,'  ebl= ',f9.5,'  eave= ',f9.5,
     &  '  err= ',f9.5,'  disp= ', f9.5,'  nwalk=',i4,' et=',f9.5)
c
	enddo
c......................................................done
1001    format ('System: 2e Z= ',f5.2,'  tau= ',f6.3, ' walkers ',i4)
1002    format ('Blocks=', i3, ' steps/bl =',i4,' skip ',i4,' blocks',/)
	end	
c
	include 'ulec.f'