   program hos
   implicit none
   integer, parameter :: i4=selected_int_kind(9)
   integer, parameter :: i8=selected_int_kind(15)
   integer, parameter :: r8=selected_real_kind(15,9)
   integer, parameter :: nwalkmax=1000000
   real(kind=r8), dimension(nwalkmax) :: xo,xn
   real(kind=r8) :: omega,dt,etrial,xnew,prob,psiold,psinew
   integer(kind=i4) :: idmc
   integer(kind=i4) :: istep
   integer(kind=i4) :: nwalk0,nwalk,nw,neq,nav,nstep,nb
   integer(kind=i4) :: i,j,k,n
   real(kind=r8) :: esum,e2sum,eblock,e2block,wsum,wblock
   real(kind=r8) :: eloc,sigma,weight,ene,ene2,err
   real(kind=r8),allocatable :: ebl(:),ebl2(:),err2(:)
   integer(kind=i4) :: iwt
   real(kind=r8) :: rgauss,csi
   integer(kind=i8) :: irn
   real(kind=r8) :: alpha
   integer(kind=i4) :: nmax,nmin,nacc
   integer(kind=i4) :: idx

   open(unit=5,file='harmonic.in',status='old')
   read(5,*) idmc   ! 0 VMC 1 DMC
   read(5,*) irn    ! number used to initialize the random number generator
   read(5,*) omega  ! strenght of the harmonic potential
   read(5,*) alpha  ! alpha parameter of the trial wave function, psi=exp(-alpha*x**2)
   read(5,*) dt     ! time step
   read(5,*) etrial ! trial energy
   read(5,*) nwalk0 ! number of walkers
   read(5,*) neq    ! number of blocks for equilibration
   read(5,*) nav    ! number of blocks for average
   read(5,*) nstep  ! number of steps per block

   write(6,'(''One dimensional harmonic oscillator, DMC'')')
   write(6,'(''Random number seed = '',t40,i20)') irn
   write(6,'(''Omega strenght of the potential = '',t50,f10.5)') omega
   write(6,'(''Alpha parameter of the trial wave function = '',t50,f10.5)') alpha
   write(6,'(''Time step = '',t50,f10.5)') dt
   write(6,'(''Trial energy = '',t50,f10.5)') etrial
   write(6,'(''Walkers = '',t50,i10)') nwalk0
   write(6,'(''Number of blocks for equilibration = '',t50,i10)') neq
   write(6,'(''Number of blocks for statistics = '',t50,i10)') nav
   write(6,'(''Number of steps for each block = '',t50,i10)') nstep
   write(6,*) 
   write(6,'(''iteration  energy block   errror         total energy    error'')')

   sigma=sqrt(dt)
   allocate(ebl(nav+neq),ebl2(nav+neq),err2(nav+neq))
   if (nwalk0.ge.nwalkmax) then
      write(6,*) 'Walkers number larger than maximum'
      stop
   endif

! Walkers initializations: give a random position to each walker
   do i=1,nwalk0
      call ran(csi,irn)
      xo(i)=1.0_r8/omega*(0.5-csi)
   enddo

! Initialize all the quantities
   nwalk=nwalk0
   esum=0.0_r8
   e2sum=0.0_r8
   wsum=0.0_r8
   nb=0
   nacc=0
   istep=0
! Loop over equilibration and averaging blocks
   do i=1,nav+neq
      if (i.eq.neq+1) then
         esum=0.0_r8
         e2sum=0.0_r8
         wsum=0.0_r8
         nb=0
         nacc=0
         write(6,*) 'Equilibration done!'
         write(6,*)
      endif
      eblock=0.0_r8
      e2block=0.0_r8
      wblock=0.0_r8
      do j=1,nstep
         istep=istep+1
         nw=0
! Loop over the walkers
         do n=1,nwalk
            if (idmc.eq.0) then ! Variational Monte Carlo algotithm
                nw=nwalk
                xnew=xo(n)+sigma*rgauss(irn)                                         !Gaussian shift of the coordinate
                psiold=exp(-0.5_r8*alpha*xo(n)**2)                                   !Psi_T(x)
                psinew=exp(-0.5_r8*alpha*xnew**2)                                    !Psi_T(x')
                prob=psinew**2/psiold**2                                             !Ratio of |Psi_T(x')|**2/|Psi_T(x)|**2
                call ran(csi,irn)                                                    !Sample the random number for the acceptance
                if (csi.lt.prob) then                                                !The proposed move is accepted
                   eloc=-0.5_r8*((alpha*xnew)**2-alpha)+0.5_r8*omega**2*xnew**2      !Compute E_L(x')
                   xn(n)=xnew
                   nacc=nacc+1
                else                                                                 !The proposed move is rejected
                   eloc=-0.5_r8*((alpha*xo(n))**2-alpha)+0.5_r8*omega**2*xo(n)**2    !Compute E_L(x)
                   xn(n)=xo(n)
                endif 
                weight=1.0_r8                                                        !The weight in VMC is always 1

            elseif (idmc.eq.1) then                                                  !Psitrial weight importance sampling algorithm 
                xnew=xo(n)+sigma*rgauss(irn)                                         !Gaussian shift of the coordinates  
                eloc=-0.5_r8*((alpha*xnew)**2-alpha)+0.5_r8*omega**2*xnew**2         !Compute E_L(x')
                psiold=exp(-0.5_r8*alpha*xo(n)**2)                                   !Psi_T(x)
                psinew=exp(-0.5_r8*alpha*xnew**2)                                    !Psi_T(x')
                weight=exp(-dt*(0.5_r8*omega**2*xnew**2-etrial))*psinew/psiold       !Importance sampled psitrial weight
                call ran(csi,irn)                                                    !Sample random number for the branching
                iwt=weight+csi                                                       !Integer number of copies of the walker
                do k=1,iwt
                   nw=nw+1
                   xn(nw)=xnew
                enddo

            elseif (idmc.eq.2) then                                                  !Elocal importance sampling with drift
                xnew=xo(n)-sigma**2*alpha*xo(n)+sigma*rgauss(irn)                    !Drifted gaussian shift of the coordinates
                eloc=-0.5_r8*((alpha*xnew)**2-alpha)+0.5_r8*omega**2*xnew**2         !Compute E_L(x') 
                weight=exp(-dt*(eloc-etrial))                                        !The weight is given by the local energy
                call ran(csi,irn)                                                    !Sample random number for the branching
                iwt=weight+csi                                                       !Integer number of copies of the walker
                do k=1,iwt
                   nw=nw+1
                   xn(nw)=xnew
                enddo
            endif
            esum=esum+eloc                                                           !Accumulate E_L
            e2sum=e2sum+eloc**2                                                      !Accumulate E_L**2
            eblock=eblock+eloc                                                       !Accumulate E_L for the block
            e2block=e2block+eloc**2                                                  !Accumulate E_L**2 for the block   
            wsum=wsum+weight                                                         !Accumulate denominator
            wblock=wblock+weight                                                     !Accumulate denominator for the block
         enddo

!Copy the new array into the old one
         do n=1,nw                                                                   
            xo(n)=xn(n)
         enddo
         nwalk=nw
         nb=nb+1
      enddo
!Compute all the expectation values and errors
      ebl(i)=eblock/wblock                                                           
      ebl2(i)=e2block/wblock
      err2(i)=sqrt(abs(ebl2(i)-ebl(i)**2)/max(1,nstep))
      ene=esum/wsum
      ene2=e2sum/wsum
      err=sqrt(abs(ene2-ene**2)/max(1,nb))

      write(6,'(i7,5f15.10,i5)') nb,ebl(i),err2(i),ene,err,real(nacc/wsum),nwalk     !Write the output
      if  (i.ge.neq+1) then
          write(1006,'(i7,4f15.10,i6,2f12.5,3i10)') nb,ebl(i),err2(i),ene,err
      endif
    enddo

    end program

!Box-Mueller method for gaussianly distributed random number
    function rgauss(irn)                                                             
    implicit none
    integer, parameter :: i8=selected_int_kind(15)
    integer, parameter :: r8=selected_real_kind(15,9)
    real(kind=r8) :: pi
    real(kind=r8) :: rgauss,x1,x2
    integer(kind=i8) :: irn
    call ran(x1,irn)
    call ran(x2,irn)
    pi=4.0_r8*atan(1.0_r8)
    rgauss=sqrt(-2.0_r8*log(x1))*cos(2.0_r8*pi*x2)
    return
    end

!Random number generator for number uniformly distributed between 0 and 1    
    subroutine ran(rn,irn)                                                           
    implicit none
    integer, parameter :: i4=selected_int_kind(9)
    integer, parameter :: i8=selected_int_kind(15)
    integer, parameter :: r8=selected_real_kind(15,9)
    integer(kind=i8),  parameter :: mask24 = ishft(1_i8,24)-1
    integer(kind=i8),  parameter :: mask48 = ishft(1_i8,48_i8)-1_i8
    real(kind=r8),  parameter :: twom48=2.0_r8**(-48)
    integer(kind=i8),  parameter :: mult1 = 44485709377909_i8
    integer(kind=i8),  parameter :: m11 = iand(mult1,mask24)
    integer(kind=i8),  parameter :: m12 = iand(ishft(mult1,-24),mask24)
    integer(kind=i8),  parameter :: iadd1 = 96309754297_i8
    integer(kind=i8) :: irn
    real(kind=r8) :: rn
    integer(kind=i8) :: is1,is2
    is2 = iand(ishft(irn,-24),mask24)
    is1 = iand(irn,mask24)
    irn = iand(ishft(iand(is1*m12+is2*m11,mask24),24)+is1*m11+iadd1,mask48)
    rn = ior(irn,1_i8)*twom48
    return
    end subroutine ran

