Brock University - Department of Chemistry

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

University Page

Faculty Page

Chemistry Page

This page is: http://chemiris.labs.brocku.ca/~chemweb/faculty/rothstein/jce/hatom.html