c.......Diffusion Monte Carlo codes for hydrogen atom
        program hatom
c.......simple diffusion Monte Carlo program for the H atom
        implicit double precision(a-h,o-z)
        integer *4 is
        dimension r(2000,3),rp(2000,3),b(2000),energy(2000)
c.......random number seed
        is=54833
        print*,"random number seed=",is
c.......initialize random number generator
        call SRAND4(is)
c.......itf is number of iterations per block
        itf=500
        pi=4.*atan(1.)
c.......tau is time step
        tau=0.01d0
        sq=sqrt(tau)
c.......m is target  ensemble size
        m=1000
c.......mp is current ensemble size
        mp=m
c.......generate random ensemble
        do i=1,m
          do j=1,3
             r(i,j)=(drand4()-0.5d0)*.6d0
          enddo
        enddo
c.......nbk is number of blocks
        nbk=11
c.......loop over blocks
        do 10 ibk=1,nbk
        enrg=0
c.......loop over iterations
        do 22 it=1,itf
c.......loop over walkers
          do i=1,mp
c.......dis is distance from nucleus
          dis=sqrt(r(i,1)**2+r(i,2)**2+r(i,3)**2)
          el=-1/dis
c........do 'move'
            do j=1,3
c........bm is standard normal pseudo-random variate
            bm=sqrt(-2*dlog(drand4()))*cos(pi*drand4())
            r(i,j)=r(i,j)+sq*bm
            enddo
          dis=sqrt(r(i,1)**2+r(i,2)**2+r(i,3)**2)
          el=-1/dis+el
          b(i)=exp(-el/2*tau)
c.........end loop over walkers
          enddo
c.......compute average branching factor
        bbar=0
          do i=1,mp
          bbar=bbar+b(i)
          enddo
        bbar=bbar/mp
c.........do branching
          l=0
c.........an ensemble of size l will result from branching
c.............each walker will be located at rp
          do i=1,mp
          nrep=int((b(i)/mp)*(m/bbar)+drand4())
          if (nrep.eq.0) go to 55
            do k=1,nrep
            l=l+1
               do j=1,3
               rp(l,j)=r(i,j)
               enddo
            enddo
 55        continue
          enddo
c.......put surviving configurations back into r matrix
        mp=l
c.......en is energy for current iteration
        en=-dlog(bbar)/tau
        enrg=enrg+en
           do i=1,mp
              do j=1,3
              r(i,j)=rp(i,j)
              enddo
           enddo
c.......end loop over iterations
 22     continue
c.......energy is average block energy
        energy(ibk)=enrg/itf
        print*,"block ave energy",energy(ibk)
c.......end loop over blocks
 10      continue
c.......compute grand mean and standard error of block energies
        se=0
        se2=0
c.......disallow first block for equilibration
        neq=1
        n=nbk-neq
        do i=neq+1,nbk
        se=se+energy(i)
        se2=se2+energy(i)**2
        enddo
        average=se/n
        se=sqrt((se2-se**2/n)/(n-1)/n)
        print*,"# blocks & grand mean",n,average
        print*,"standard error",se
        end


c.......Diffusion Monte Carlo codes for hydrogen molecule and hydrogen molecular ion

        program h2h2p
c.......simple diffusion Monte Carlo program for H2
c.......to do H2+ refer to c**** lines of code
        implicit double precision(a-h,o-z)
        integer *4 is
        dimension r(3000,2,3),rp(3000,2,3),tot(3000),b(3000)
c.......random number seed
        is=64801
        print*,"random number seed=",is
c.......initialize random number generator
        call SRAND4(is)
c.......itf is number of iterations per block
        itf=500
        pi=4.*atan(1.)
c.......tau is time step
        tau=.01d0
        sq=sqrt(tau)
c.......m is target ensemble size
        m=1000
c.......mp is current ensemble size
        mp=m
c.......ne is number of electrons
        ne=2
c*******replace above with ne=1
c.......hrnn is half the internuclear separation
        hrnn=0.7
c*******replace above with hrnn=1.d0
c.......generate random ensemble
        do j=1,3
         do k=1,ne
          do i=1,m
           r(i,k,j)=(drand4()-0.5d0)*1.d0
          enddo
         enddo
        enddo
c.......nbk is number of blocks
        nbk=11
c.......loop over blocks
        do 10 ibk=1,nbk
        enrg=0
c.......loop over iterations
        do 22 it=1,itf
c.......loop over walkers
          do i=1,mp
c.......disAi and disBi are distances of electron i
c...........from nuclei A and B, respectively
c.......dis12 is interelectron distance
          disA1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)+hrnn)**2)
          disB1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)-hrnn)**2)
c*******omit following 6 lines for H2+
          disA2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)+hrnn)**2)
          disB2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)-hrnn)**2)
          dis12=sqrt((r(i,1,1)-r(i,2,1))**2+(r(i,1,2)-r(i,2,2))**2
     + +(r(i,1,3)-r(i,2,3))**2)
          pot=-1.d0/disA1-1.d0/disB1-1.d0/disA2-1.d0/disB2
     + +1.d0/dis12
c********replace above with pot=-1.d0/disA1-1.d0/disB1
c......do 'move'
             do j=1,3
               do k=1,ne
c......bm is standard normal pseudo-random variable
               bm=sqrt(-2*dlog(drand4()))*cos(pi*drand4())
               r(i,k,j)=r(i,k,j)+sq*bm
               enddo
             enddo
          disA1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)+hrnn)**2)
          disB1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)-hrnn)**2)
c*******omit following 6 lines for H2+
          disA2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)+hrnn)**2)
          disB2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)-hrnn)**2)
          dis12=sqrt((r(i,1,1)-r(i,2,1))**2+(r(i,1,2)-r(i,2,2))**2
     + +(r(i,1,3)-r(i,2,3))**2)
          pot=-1.d0/disA1-1.d0/disB1-1.d0/disA2-1.d0/disB2
     +  +1.d0/dis12+pot
c********replace above with pot=-1.d0/disA1-1.d0/disB1+pot
          b(i)=exp(-pot/2*tau)
c......end loop over walkers
          enddo
c.......compute average branching factor
        bbar=0
          do i=1,mp
          bbar=bbar+b(i)
          enddo
        bbar=bbar/mp
c.......do branching
        l=0
c.......an ensemble of size l will result from branching
c.........each walker will be located at rp
          do i=1,mp
          nrep=int((b(i)/mp)*(m/bbar)+drand4())
          if (nrep.eq.0) go to 55
            do kk=1,nrep
            l=l+1
              do j=1,3
               do k=1,ne
               rp(l,k,j)=r(i,k,j)
               enddo
              enddo
            enddo
 55     continue
          enddo
c.......put surviving configurations back into r matrix
        mp=l
c.......en is energy for current iteration
        en=-dlog(bbar)/tau
        enrg=enrg+en
         do j=1,3
           do k=1,ne
             do i=1,mp
              r(i,k,j)=rp(i,k,j)
             enddo
           enddo
         enddo
c.......end loop over iterations
 22     continue
c.......tot is average total energy of current block
        tot(ibk)=enrg/itf+.5/hrnn
        print*,"block ave total energy",tot(ibk)
c.......end loop over blocks
 10     continue
c.......compute grand mean and standard error of block energies
        se=0
        se2=0
c.......disallow first block for equilibration
        neq=1
        n=nbk-neq
        do i=neq+1,nbk
        se=se+tot(i)
        se2=se2+tot(i)**2
        enddo
        average=se/n
        se=sqrt((se2-se**2/n)/(n-1)/n)
        print*,"# blocks & grand mean",n,average
        print*,"standard error",se
        end