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