************************************************************************
C                        METROP_LAYR1.F
************************************************************************
C Implements the multilevel metropolis algorithm.
C Uses 6 term polynomial expansion of the particle interaction.
C Uses perm,mgauss,samdist,ran2
************************************************************************
      subroutine metrop(x,y,z,xp,yp,zp,r1,rgrlat,rgroff,
     x     dutfix,ddutfix,u00,ugrlat,ugroff,image,
     x     idum,nd,ngrlat,ngroff,icount2,icount3,
     x     icount4,nacc)
      include 'pimc_sqr3sc.dcl'

      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension xp(2*m,nat),yp(2*m,nat),zp(2*m,nat)
      dimension t(nperm),sum(0:nperm) 
      dimension t1(nperm),sum1(0:nperm)
      dimension xbar0(4),ybar0(4),zbar0(4)
      dimension xpbar0(4),ypbar0(4),zpbar0(4)
      dimension r1(nr1),u00(nr1)
      dimension rgrlat(nr1),rgroff(nr1)
      dimension ugrlat(nr1),ugroff(nr1)
      dimension dutfix(lev,nr1),ddutfix(lev,nr1)
	
      integer matoms(nperm),iatom(nperm,4),iatomp(nperm,4)
      integer matoms1(nperm),iatom1(nperm,4),iatomp1(nperm,4)
      integer index1(2*m),index2(3*m)
      integer image(2*m,nat),isave1(4),isave2(4)

      dimension tarray(2)
      
*///////////////////////////////////////////////////////////////////////
C Useful constants.  
*///////////////////////////////////////////////////////////////////////

      cutoff=exp(-argcut)
      pi=4.0*atan(1.0)
      con=1./(4.*alambda*float(nbis)*tau)
      icount=0
      icntret=0
      icntmax=1
      time_mg=0.
      time_sd=0.
      do i=1,2*m
         index2(i)=i
      end do
      do i=1,m
         index1(i)=i
         index1(i+m)=i
         index2(i+2*m)=i
      end do

*///////////////////////////////////////////////////////////////////////
C Loop over time slices.  Perm() creates permutation table.  Mgauss() 
C proposes new particle positions and calculates the ratio of transition 
C amplitudes.  Samdist calculates probability amplitudes for a given 
C configuration.
*///////////////////////////////////////////////////////////////////////

      icheckit=0
      do islice=1,m
         icntret=icntret+icount
         icount=0
         ichol=0
         ipc=1
         call perm(x,y,z,t,matoms,iatom,iatomp,sum,con,
     x        islice,iptot)
         xnormi=sum(iptot)
            
C Loop over many particles per slice

         do imov=1,nmov               ! Upper limit here is arbitrary
            ireject=0

C Choose a permutation

            xr=ran2(idum)       ! Choose permutation     
            chi=sum(iptot)*xr
            do iperm=1,iptot
               if(sum(iperm-1).lt.chi.and.chi.le.sum(iperm)) then
                  ipc=iperm
                  goto 30
               end if
            end do
 30         continue

            sumtmp=0.0
            do ii=1,matoms(ipc)                 
               ia=iatom(ipc,ii)
               iap=iatomp(ipc,ii)
               iarg=islice+nbis
                                ! Record permutation
               xp(islice,ia)=x(islice,ia)
               yp(islice,ia)=y(islice,ia)
               zp(islice,ia)=z(islice,ia)
               xp(iarg,iap)=x(iarg,ia)
               yp(iarg,iap)=y(iarg,ia)
               zp(iarg,iap)=z(iarg,ia)
                                ! Calculate identity permutation
               call dstpbc(x(iarg,ia),x(islice,ia),xdim,tt1,sgn)
               call dstpbc(y(iarg,ia),y(islice,ia),ydim,tt2,sgn)
               tt3=z(iarg,ia)-z(islice,ia)
               dst=tt1*tt1+tt2*tt2+tt3*tt3
               sumtmp=sumtmp+con*dst
c               tidp=exp(-con*dst)
            end do
            tidp=exp(-sumtmp)
            
C Begin bisection routine
            delsumold=log(t(ipc)/tidp)
            natp=matoms(ipc)
c            if(natp.gt.1) print*,natp,t(ipc),tidp,delsumold
            do k=1,lev
               delsum=0.0
               tsum=0.
               iimax=2**(k-1)
               inc=nbis/iimax
               nsl1=inc
               nsl2=0
               nsld=(nsl1-nsl2)/2
               nsla=(nsl1+nsl2)/2
               sig2=float(nsld)*alambda*tau
               do ii1=1,iimax
                  do ll=1,natp
                     ia=iatom(ipc,ll)
                     iap=iatomp(ipc,ll)
                     is1=islice+nsl1
                     is2=islice+nsl2
                     call dstpbc(xp(is1,iap),xp(is2,iap),
     x                    xdim,dstx,sgnx)
                     call dstpbc(yp(is1,iap),yp(is2,iap),
     x                    ydim,dsty,sgny)
                     xpbar0(ll)=xp(is1,iap)-sgnx*dstx/2.
                     ypbar0(ll)=yp(is1,iap)-sgny*dsty/2.
                     zpbar0(ll)=(zp(is2,iap)+zp(is1,iap))/2.
                     if(xpbar0(ll).gt.xdim) xpbar0(ll)=xpbar0(ll)-xdim
                     if(xpbar0(ll).lt.0.) xpbar0(ll)=xpbar0(ll)+xdim
                     if(ypbar0(ll).gt.ydim) ypbar0(ll)=ypbar0(ll)-ydim
                     if(ypbar0(ll).lt.0.) ypbar0(ll)=ypbar0(ll)+ydim

                     call dstpbc(x(is1,ia),x(is2,ia),xdim,dstx,sgnx)
                     call dstpbc(y(is1,ia),y(is2,ia),ydim,dsty,sgny)
                     xbar0(ll)=x(is1,ia)-sgnx*dstx/2.
                     ybar0(ll)=y(is1,ia)-sgny*dsty/2.
                     zbar0(ll)=(z(is2,ia)+z(is1,ia))/2.
                     if(xbar0(ll).gt.xdim) xbar0(ll)=xbar0(ll)-xdim
                     if(xbar0(ll).lt.0.) xbar0(ll)=xbar0(ll)+xdim
                     if(ybar0(ll).gt.ydim) ybar0(ll)=ybar0(ll)-ydim
                     if(ybar0(ll).lt.0.) ybar0(ll)=ybar0(ll)+ydim
                  end do
                  move=1
c                  start1=etime(tarray)
                  call mgauss(xp,yp,zp,xpbar0,ypbar0,zpbar0,
     x                 r1,dutfix,ddutfix,
     x                 iatomp,idum,natp,sig2,nd,ipc,islice,nsla,k,
     x                 move,tpexp,detp)
c                  finish1=etime(tarray)
c                  time_mg=time_mg+finish1-start1
                  if(move.gt.1) then
                     ireject=1
                     goto 998
                  end if
                  move=0
                  call mgauss(x,y,z,xbar0,ybar0,zbar0,
     x                 r1,dutfix,ddutfix,iatom,idum,
     x                 natp,sig2,nd,ipc,islice,nsla,k,move,texp,det)
                  if(move.gt.1) then
                     ireject=1
                     goto 998
                  end if
                  detexp=-.5*(log(det)-log(detp))
                  tsum=tsum+texp-tpexp+detexp
                  move=1
c                  start2=etime(tarray)
                  call samdist(xp,yp,zp,r1,rgrlat,rgroff,
     x                 u00,ugrlat,ugroff,iatomp,nsl1,
     x                 nsl2,k,islice,nd,ngrlat,ngroff,ipc,
     x                 natp,move,xnum)
c                  finish2=etime(tarray)
c                  time_sd=time_sd+finish2-start2
                  if(move.gt.1) then
                     ireject=1
                     goto 998
                  end if
                  move=0
                  call samdist(x,y,z,r1,rgrlat,rgroff,
     x                 u00,ugrlat,ugroff,iatom,
     x                 nsl1,nsl2,k,islice,nd,
     x                 ngrlat,ngroff,ipc,natp,move,denom)
c                  print*,texp,tpexp,xnum,denom
                  arg=-xnum+denom
                  delsum=delsum+arg
                  nsl2=nsl2+inc
                  nsl1=nsl1+inc
                  nsld=(nsl1-nsl2)/2
                  nsla=(nsl1+nsl2)/2
               end do           ! End bisection loop
               expsum=tsum+delsum-delsumold
c               print*,expsum
c               print*

               if(expsum.lt.0.)then ! Conditionally accept configs.
                  accept=exp(expsum)
                  xr=ran2(idum)
                  if(xr.gt.accept) then
                     ireject=1
                     goto 998
                  end if
               endif
               delsumold=delsum
            end do              ! End loop over nbis levels
            
C Final test:  check acceptence of original permutation.  Not necessary for
C identity permutations.
            
            if(matoms(ipc).gt.1) then 
               call perm(xp,yp,zp,t1,matoms1,iatom1,iatomp1,sum1,
     x              con,islice,iptot1)
               xnormp=sum1(iptot1)
c               print*,xnormi,xnormp,xnormi/xnormp
            else
               xnormp=xnormi
            end if

c            ratio=xnormp/xnormi
            ratio=xnormi/xnormp
            xr=ran2(idum)
            if(xr.gt.ratio) then 
               ireject=1
               goto 998
            end if

            nacc=nacc+1
            do jj=1,matoms(ipc)                ! Update positions
               ia=iatom(ipc,jj)
               iap=iatomp(ipc,jj)
               do ii=islice+nbis+1,islice+m-1
                  xp(ii,iap)=x(ii,ia)
                  yp(ii,iap)=y(ii,ia)
                  zp(ii,iap)=z(ii,ia)
               end do
            end do
            do ii=islice,islice+m-1
               ii1=index1(ii)
               ii2=index2(ii+m)
               do jj=1,matoms(ipc)
                  ia=iatom(ipc,jj)
                  iap=iatomp(ipc,jj)
                  do kk=1,nat
                     if(image(ii1,kk).eq.ia) then
                        isave1(jj)=iap
                        isave2(jj)=kk
                        goto 50
                     end if
                  end do
 50               continue
               end do
               do jj=1,matoms(ipc)
                  image(ii1,isave2(jj))=isave1(jj)
                  image(ii1+m,isave1(jj))=isave2(jj)
               end do
               do kk=1,nat
                  ik=image(ii,kk)
                  xp(ii2,ik)=xp(ii,kk)
                  yp(ii2,ik)=yp(ii,kk)
                  zp(ii2,ik)=zp(ii,kk)
               end do
            end do
            do ii=1,m
               do jj=1,nat
                  x(ii,jj)=xp(ii,jj)
                  x(ii+m,jj)=xp(ii+m,jj)
                  y(ii,jj)=yp(ii,jj)
                  y(ii+m,jj)=yp(ii+m,jj)
                  z(ii,jj)=zp(ii,jj)
                  z(ii+m,jj)=zp(ii+m,jj)
               end do
            end do
                                        ! Update permutation stuff
            isaveit=matoms(ipc)
            if(isaveit.gt.1) then
               if(isaveit.eq.2) icount2=icount2+1
               if(isaveit.eq.3) icount3=icount3+1
               if(isaveit.eq.4) icount4=icount4+1
               xnormi=xnormp
               iptot=iptot1
               sum(0)=sum1(0)   
               do iperm=1,iptot
                  t(iperm)=t1(iperm)
                  matoms(iperm)=matoms1(iperm)
                  sum(iperm)=sum1(iperm)
                  do jj=1,matoms(iperm)
                     iatom(iperm,jj)=iatom1(iperm,jj)
                     iatomp(iperm,jj)=iatomp1(iperm,jj)
                  end do
               end do
            end if

 998        continue 

            if(ireject.eq.1) then            ! Handle rejection
c               if(matoms(ipc).gt.1.and.t(ipc).gt.tidp) then
c                  print*,'move rejected at level',k
c               end if
c               if(matoms(ipc).eq.2) icount2=icount2+1
c               if(matoms(ipc).eq.3) icount3=icount3+1
c               if(matoms(ipc).eq.4) icount4=icount4+1
               do ia=0,nbis   
                  iarg=islice+ia
                  do ib=1,matoms(ipc)
                     ic=iatom(ipc,ib)
                     xp(iarg,ic)=x(iarg,ic)
                     yp(iarg,ic)=y(iarg,ic)
                     zp(iarg,ic)=z(iarg,ic)
                  end do
               end do
            end if

         end do                 ! End loop imov
      end do                    ! End loop over time slices
c      print*,time_mg,time_sd
 999  continue
      return
      end
