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