Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

dmc

By: a guest on Jun 25th, 2012  |  syntax: Fortran  |  size: 7.36 KB  |  views: 37  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. c.......Diffusion Monte Carlo codes for hydrogen atom
  2.         program hatom
  3. c.......simple diffusion Monte Carlo program for the H atom
  4.         implicit double precision(a-h,o-z)
  5.         integer *4 is
  6.         dimension r(2000,3),rp(2000,3),b(2000),energy(2000)
  7. c.......random number seed
  8.         is=54833
  9.         print*,"random number seed=",is
  10. c.......initialize random number generator
  11.         call SRAND4(is)
  12. c.......itf is number of iterations per block
  13.         itf=500
  14.         pi=4.*atan(1.)
  15. c.......tau is time step
  16.         tau=0.01d0
  17.         sq=sqrt(tau)
  18. c.......m is target  ensemble size
  19.         m=1000
  20. c.......mp is current ensemble size
  21.         mp=m
  22. c.......generate random ensemble
  23.         do i=1,m
  24.           do j=1,3
  25.              r(i,j)=(drand4()-0.5d0)*.6d0
  26.           enddo
  27.         enddo
  28. c.......nbk is number of blocks
  29.         nbk=11
  30. c.......loop over blocks
  31.         do 10 ibk=1,nbk
  32.         enrg=0
  33. c.......loop over iterations
  34.         do 22 it=1,itf
  35. c.......loop over walkers
  36.           do i=1,mp
  37. c.......dis is distance from nucleus
  38.           dis=sqrt(r(i,1)**2+r(i,2)**2+r(i,3)**2)
  39.           el=-1/dis
  40. c........do 'move'
  41.             do j=1,3
  42. c........bm is standard normal pseudo-random variate
  43.             bm=sqrt(-2*dlog(drand4()))*cos(pi*drand4())
  44.             r(i,j)=r(i,j)+sq*bm
  45.             enddo
  46.           dis=sqrt(r(i,1)**2+r(i,2)**2+r(i,3)**2)
  47.           el=-1/dis+el
  48.           b(i)=exp(-el/2*tau)
  49. c.........end loop over walkers
  50.           enddo
  51. c.......compute average branching factor
  52.         bbar=0
  53.           do i=1,mp
  54.           bbar=bbar+b(i)
  55.           enddo
  56.         bbar=bbar/mp
  57. c.........do branching
  58.           l=0
  59. c.........an ensemble of size l will result from branching
  60. c.............each walker will be located at rp
  61.           do i=1,mp
  62.           nrep=int((b(i)/mp)*(m/bbar)+drand4())
  63.           if (nrep.eq.0) go to 55
  64.             do k=1,nrep
  65.             l=l+1
  66.                do j=1,3
  67.                rp(l,j)=r(i,j)
  68.                enddo
  69.             enddo
  70.  55        continue
  71.           enddo
  72. c.......put surviving configurations back into r matrix
  73.         mp=l
  74. c.......en is energy for current iteration
  75.         en=-dlog(bbar)/tau
  76.         enrg=enrg+en
  77.            do i=1,mp
  78.               do j=1,3
  79.               r(i,j)=rp(i,j)
  80.               enddo
  81.            enddo
  82. c.......end loop over iterations
  83.  22     continue
  84. c.......energy is average block energy
  85.         energy(ibk)=enrg/itf
  86.         print*,"block ave energy",energy(ibk)
  87. c.......end loop over blocks
  88.  10      continue
  89. c.......compute grand mean and standard error of block energies
  90.         se=0
  91.         se2=0
  92. c.......disallow first block for equilibration
  93.         neq=1
  94.         n=nbk-neq
  95.         do i=neq+1,nbk
  96.         se=se+energy(i)
  97.         se2=se2+energy(i)**2
  98.         enddo
  99.         average=se/n
  100.         se=sqrt((se2-se**2/n)/(n-1)/n)
  101.         print*,"# blocks & grand mean",n,average
  102.         print*,"standard error",se
  103.         end
  104.  
  105.  
  106. c.......Diffusion Monte Carlo codes for hydrogen molecule and hydrogen molecular ion
  107.  
  108.         program h2h2p
  109. c.......simple diffusion Monte Carlo program for H2
  110. c.......to do H2+ refer to c**** lines of code
  111.         implicit double precision(a-h,o-z)
  112.         integer *4 is
  113.         dimension r(3000,2,3),rp(3000,2,3),tot(3000),b(3000)
  114. c.......random number seed
  115.         is=64801
  116.         print*,"random number seed=",is
  117. c.......initialize random number generator
  118.         call SRAND4(is)
  119. c.......itf is number of iterations per block
  120.         itf=500
  121.         pi=4.*atan(1.)
  122. c.......tau is time step
  123.         tau=.01d0
  124.         sq=sqrt(tau)
  125. c.......m is target ensemble size
  126.         m=1000
  127. c.......mp is current ensemble size
  128.         mp=m
  129. c.......ne is number of electrons
  130.         ne=2
  131. c*******replace above with ne=1
  132. c.......hrnn is half the internuclear separation
  133.         hrnn=0.7
  134. c*******replace above with hrnn=1.d0
  135. c.......generate random ensemble
  136.         do j=1,3
  137.          do k=1,ne
  138.           do i=1,m
  139.            r(i,k,j)=(drand4()-0.5d0)*1.d0
  140.           enddo
  141.          enddo
  142.         enddo
  143. c.......nbk is number of blocks
  144.         nbk=11
  145. c.......loop over blocks
  146.         do 10 ibk=1,nbk
  147.         enrg=0
  148. c.......loop over iterations
  149.         do 22 it=1,itf
  150. c.......loop over walkers
  151.           do i=1,mp
  152. c.......disAi and disBi are distances of electron i
  153. c...........from nuclei A and B, respectively
  154. c.......dis12 is interelectron distance
  155.           disA1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)+hrnn)**2)
  156.           disB1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)-hrnn)**2)
  157. c*******omit following 6 lines for H2+
  158.           disA2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)+hrnn)**2)
  159.           disB2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)-hrnn)**2)
  160.           dis12=sqrt((r(i,1,1)-r(i,2,1))**2+(r(i,1,2)-r(i,2,2))**2
  161.      + +(r(i,1,3)-r(i,2,3))**2)
  162.           pot=-1.d0/disA1-1.d0/disB1-1.d0/disA2-1.d0/disB2
  163.      + +1.d0/dis12
  164. c********replace above with pot=-1.d0/disA1-1.d0/disB1
  165. c......do 'move'
  166.              do j=1,3
  167.                do k=1,ne
  168. c......bm is standard normal pseudo-random variable
  169.                bm=sqrt(-2*dlog(drand4()))*cos(pi*drand4())
  170.                r(i,k,j)=r(i,k,j)+sq*bm
  171.                enddo
  172.              enddo
  173.           disA1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)+hrnn)**2)
  174.           disB1=sqrt(r(i,1,1)**2+r(i,1,2)**2+(r(i,1,3)-hrnn)**2)
  175. c*******omit following 6 lines for H2+
  176.           disA2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)+hrnn)**2)
  177.           disB2=sqrt(r(i,2,1)**2+r(i,2,2)**2+(r(i,2,3)-hrnn)**2)
  178.           dis12=sqrt((r(i,1,1)-r(i,2,1))**2+(r(i,1,2)-r(i,2,2))**2
  179.      + +(r(i,1,3)-r(i,2,3))**2)
  180.           pot=-1.d0/disA1-1.d0/disB1-1.d0/disA2-1.d0/disB2
  181.      +  +1.d0/dis12+pot
  182. c********replace above with pot=-1.d0/disA1-1.d0/disB1+pot
  183.           b(i)=exp(-pot/2*tau)
  184. c......end loop over walkers
  185.           enddo
  186. c.......compute average branching factor
  187.         bbar=0
  188.           do i=1,mp
  189.           bbar=bbar+b(i)
  190.           enddo
  191.         bbar=bbar/mp
  192. c.......do branching
  193.         l=0
  194. c.......an ensemble of size l will result from branching
  195. c.........each walker will be located at rp
  196.           do i=1,mp
  197.           nrep=int((b(i)/mp)*(m/bbar)+drand4())
  198.           if (nrep.eq.0) go to 55
  199.             do kk=1,nrep
  200.             l=l+1
  201.               do j=1,3
  202.                do k=1,ne
  203.                rp(l,k,j)=r(i,k,j)
  204.                enddo
  205.               enddo
  206.             enddo
  207.  55     continue
  208.           enddo
  209. c.......put surviving configurations back into r matrix
  210.         mp=l
  211. c.......en is energy for current iteration
  212.         en=-dlog(bbar)/tau
  213.         enrg=enrg+en
  214.          do j=1,3
  215.            do k=1,ne
  216.              do i=1,mp
  217.               r(i,k,j)=rp(i,k,j)
  218.              enddo
  219.            enddo
  220.          enddo
  221. c.......end loop over iterations
  222.  22     continue
  223. c.......tot is average total energy of current block
  224.         tot(ibk)=enrg/itf+.5/hrnn
  225.         print*,"block ave total energy",tot(ibk)
  226. c.......end loop over blocks
  227.  10     continue
  228. c.......compute grand mean and standard error of block energies
  229.         se=0
  230.         se2=0
  231. c.......disallow first block for equilibration
  232.         neq=1
  233.         n=nbk-neq
  234.         do i=neq+1,nbk
  235.         se=se+tot(i)
  236.         se2=se2+tot(i)**2
  237.         enddo
  238.         average=se/n
  239.         se=sqrt((se2-se**2/n)/(n-1)/n)
  240.         print*,"# blocks & grand mean",n,average
  241.         print*,"standard error",se
  242.         end