1      program he
2C$Id: mc.f,v 1.2 1995-02-02 23:24:16 d3g681 Exp $
3      implicit double precision (a-h, o-z)
4      include 'msgtypesf.h'
5      parameter (maxpt=512)
6      dimension x(6,maxpt), psix(maxpt), ex(maxpt)
7      data esq/0.0d0/
8c
9c     Monte Carlo Test Program ... evaluates the expectation
10c     value of the energy for a simple wavefunction for helium.
11c
12c     Initalize the message passing environment
13c
14      call pbeginf
15      call evon
16c
17c     process zero reads in input and then broadcasts to others
18c
19      if (nodeid().eq.0) then
20        write(6,1)
211       format(' He atom variational monte carlo '//
22     &         ' Input neq and nstep '$)
23c       call flush(6)
24        read(5,*) neq, nstep
25        nstep = ( (nstep+99)/100 ) * 100
26        neq = ( (neq+99)/100 ) * 100
27      endif
28      call brdcst(1+MSGINT, neq, mitob(1), 0)
29      call brdcst(1+MSGINT, nstep, mitob(1), 0)
30c
31c     Divide up the work ... each process does a subset of the points
32c     or configurations ... make sure that total is indeed maxpt
33c
34      nproc = nnodes()
35      npoint = maxpt / nproc
36      nremain = maxpt - npoint*nproc
37      if (nodeid().lt.nremain) npoint = npoint + 1
38c
39c     Actually do the work. Routine init generates the intial points
40c     in a 2x2x2 cube. Equilibriate for neq moves then compute averages
41c     for nstep moves.
42c
43      rjunk = timer()
44      call init(npoint, x, psix)
45      call mover(x, psix, ex, e, esq, npoint, neq)
46      call mover(x, psix, ex, e, esq, npoint, nstep)
47c
48c     Sum the results from all the processors
49c
50      call dgop(2+MSGDBL, e, 1, '+')
51      call dgop(3+MSGDBL, esq, 1, '+')
52c
53c     Write out the results before terminating
54c
55      if (nodeid().eq.0) then
56         e = e / dble(nproc)
57         esq = esq / dble(nproc)
58         err = sqrt((esq-e*e) / dble(nproc*nstep/100))
59         used = timer()
60         write(6,2) e, err, used, nproc
61 2       format(' energy =',f10.6,' +/-',f9.6,': used =',f7.2,' secs',
62     $        ': nproc =',i3)
63      endif
64c
65      call pend
66      call fexit
67c
68      end
69      subroutine mover(x, psix, ex, e, esq, npoint, nstep)
70      implicit double precision (a-h, o-z)
71      dimension x(6,npoint), psix(npoint), ex(npoint)
72c
73c     move the set of points nstep times accumulating averages
74c     for the energy and square of the energy
75c
76      e = 0.0d0
77      esq = 0.0d0
78      eb = 0.0d0
79      do 10 istep = 1, nstep
80c
81c     sample a new set of points
82c
83         do 20 ipoint = 1, npoint
84            call sample(x(1, ipoint), psix(ipoint), ex(ipoint))
85 20      continue
86c
87c     accumulate average(s)
88c
89         do 30 ipoint = 1, npoint
90            eb = eb + ex(ipoint)
91 30      continue
92c
93c     block averages every 100 moves to reduce statistical correlation
94c
95         if (mod(istep,100).eq.0) then
96            eb = eb / dble(npoint*100)
97            e = e + eb
98            esq = esq + eb*eb
99            eb = 0.0d0
100         endif
101 10   continue
102c
103c     compute final averages
104c
105      e = e / dble(nstep/100)
106      esq = esq / dble(nstep/100)
107c
108      end
109      subroutine sample(x, psix, ex)
110      implicit double precision (a-h, o-z)
111      dimension x(6), xnew(6)
112c
113c     sample a new point ... i.e. move current point by a
114c     random amount and accept the move according to the
115c     ratio of the square of the wavefunction at the new
116c     point and the old point.
117c
118c     generate trial point and info at the new point
119c
120      do 10 i = 1,6
121         xnew(i) = x(i) + (drand48(0)-0.5d0)*0.3d0
122 10   continue
123      call rvals(xnew, r1, r2, r12, r1dr2)
124      psinew = psi(r1, r2, r12)
125c
126c     accept or reject the move
127c
128      prob = min((psinew / psix)**2, 1.0d0)
129      if (prob .gt. drand48(0)) then
130         do 20 i = 1,6
131            x(i) = xnew(i)
132 20      continue
133         psix = psinew
134      else
135         call rvals(x, r1, r2, r12, r1dr2)
136      endif
137      ex = elocal(r1, r2, r12, r1dr2)
138c
139      end
140      subroutine rvals(x, r1, r2, r12, r1dr2)
141      implicit double precision (a-h, o-z)
142      dimension x(6)
143c
144c     compute required distances etc.
145c
146      r1 = dsqrt(x(1)**2 + x(2)**2 + x(3)**2)
147      r2 = dsqrt(x(4)**2 + x(5)**2 + x(6)**2)
148      r12 = dsqrt((x(1)-x(4))**2 + (x(2)-x(5))**2 + (x(3)-x(6))**2)
149      r1dr2 = x(1)*x(4) + x(2)*x(5) + x(3)*x(6)
150c
151      end
152      double precision function psi(r1, r2, r12)
153      implicit double precision (a-h, o-z)
154c
155c     compute value of the trial wavefunction
156c
157      psi = dexp(-2.0d0*(r1+r2)) * (1.0d0 + 0.5d0*r12)
158c
159      end
160      double precision function elocal(r1, r2, r12, r1dr2)
161      implicit double precision (a-h, o-z)
162c
163c     compute local energy = (H psi) / psi
164c
165      f = r12*(1.0d0 + 0.5d0*r12)
166      g = 0.5d0*r12 + r1 +r2 - r1dr2*(1.0d0/r1 + 1.0d0/r2)
167      elocal = -4.0d0 + g / f
168c
169      end
170      subroutine init(npoint, x, psix)
171      implicit double precision (a-h, o-z)
172      dimension x(6,npoint), psix(npoint)
173c
174c     distribute points in a 2x2x2 cube.
175c
176c     in parallel version make random no. seed depend on the
177c     process number ... for a production code should be more
178c     sophisticated than this.
179c
180      call srand48(2*nodeid()+1)
181c
182      do 10 ipoint = 1,npoint
183         do 20 i = 1,6
184            x(i,ipoint) = (drand48(0) - 0.5d0) * 2.0d0
185 20      continue
186         call rvals(x(1,ipoint), r1, r2, r12, r1dr2)
187         psix(ipoint) = psi(r1, r2, r12)
188 10   continue
189c
190      end
191