Want more features on Pastebin? Sign Up, it's FREE!

# dmc

By: a guest on Jun 25th, 2012  |  syntax: Fortran  |  size: 7.36 KB  |  views: 38  |  expires: Never
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
clone this paste RAW Paste Data
Top