Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement