!c  ATOMO IDROGENOIDE
!!!c
!c  campionamento Metropolis variazionale di psi^2, mossa radiale
!c  H = - nabla^2 - 2Z / r		(uso unita` atomiche Ry)!
!c  psi=exp(-r / a) funz.di prova VMC	(a parametro variazionale)

      program	idrogeno
      implicit	none
      real*8	r

      call start(r)
      call vmc(r)
      stop
      end


!c  start: raccogli la posizione iniziale

      subroutine start(r)
      implicit	 none
      real*8	 r
      character*15 nomefile

      nomefile = 'idrogenoide.dat'
      print*,'Atomo idrogenoide con VMC'
      print*,'Hamiltoniana H = - nabla^2 - 2Z / r'
      print*,'Funzione di prova psi(r) = exp (-r / a)'
10    print*,'distanza iniziale dal nucleo (in raggi di Bohr) >>'
      read*,r
      if (r .eq. 0) goto 10	! no divergenze potenziale
      open(1,file=nomefile,status='unknown')
      write(1,*)' !a      <El>      sigmam'
      return
      end

      
!c  vmc: routine principale

      subroutine vmc(r)
      implicit	 none
      integer	 imcs,seed,M,Z,i
      real*8	 r,a,ai,af,deltaa,delta,el,eloc,elcum,el2cum,acc,rnew,elbar,el2bar,sigma2,tmp
      parameter	 (deltaa=0.01d0)

      print*,'n. step MC >>'
      read*,M
15    print*,'carica del nucleo in unita` e >>'
      read*,Z
      print*,'valore minimo di partenza del parametro variazionale >>'
      read*,ai
!      af = ai + 20 * deltaa
      if (Z .lt. 1) goto 15	!Z=1,2,3,...
      delta = 1.d0		!ampiezza iniziale max. passo Metropolis
      do i = 0, 20! a = ai, af, deltaa	!loop sul parametro variaz.
         a = ai + i*deltaa
         elcum  = 0.d0		!azzero le quantita` cumulative
         el2cum = 0.d0		!
         acc    = 0.d0		!
     	 el     = eloc(r,a,Z)	!inizializzo l'en.locale
         do 200 imcs = 1,M		!loop sugli step MC
            elcum  = elcum + el
            el2cum = el2cum + el*el
	    call walk(delta,el,elcum,el2cum,r,rnew)
	    call metrop(r,rnew,acc,el,a,Z)
            !regolo il passo Metropolis
            tmp = acc/imcs
            if (tmp .gt. 0.6) then
               delta=delta*2
            else if (tmp .lt. 0.4) then
               delta=delta/2
            end if
200      continue
         elbar  = elcum/M
	 el2bar = el2cum/M
	 sigma2 = ( el2bar - elbar * elbar ) / M
	 call output(a,elbar,sigma2,acc,M)
    end do
      close(1)
      stop
      end


!c  walk: fa la mossa da sottoporre a Metropolis

      subroutine walk(delta,el,elcum,el2cum,r,rnew)
      implicit	 none
      real*8	 delta,el,elcum,el2cum,r,rnew,drand48,rnd


20    call random_number(rnd)
      rnew   = abs(r + delta * ( rnd - 0.5d0 ))	! r>0
      if (rnew .eq. 0) goto 20
      return
      end


!c  metrop: accetta o rifiuta la mossa

      subroutine metrop(r,rnew,acc,el,a,Z)
      implicit	 none
      integer	 Z
      real*8	 r,rnew,acc,el,a,eloc,p,drand48,rnd

      !! Metropolis con mossa radiale: w = r^2 * psi^2	!!
      p = ((rnew/r)**2) * (exp(2*(r-rnew)/a))
      call random_number(rnd) !!
      if ( p .gt. rnd) then			!!
         r   = rnew					!!
         acc = acc+1.d0					!!
         el  = eloc(r,a,Z)				!!
      end if						!!
      return
      end


!c  output: scrivi su file e su video i risultati

      subroutine output(a,elbar,sigma2,acc,M)
      implicit	 none
      integer	 M
      real*8	 a,elbar,sigma2,acc

      write(*,50)'a=',a,' -> <E>=',elbar,'sigma=',sqrt(sigma2),' [acc=',acc/M,']'
      write(1,60)a,elbar,sqrt(sigma2)
50    format(a3,f6.2,a8,f9.5,a16,f7.4,a6,f6.4,a1)
60    format(3(1x,f8.5))
      return
      end


!c  eloc: energia locale (1/psi) H psi

      function  eloc(r,a,Z)
      implicit	none
      integer	Z
      real*8	a,r,eloc,T,V,tmp1,tmp2

      ! eloc(a,r) con psi = exp(-r/a)
      tmp1 = 1.d0/a
      tmp2 = 2.d0/r
      T = tmp1 * tmp2 - tmp1 * tmp1
      V = - Z * tmp2
      eloc = T + V
      return
      end