subroutine omove
       implicit none
       include 'pimc.par'
       integer il,l,js,jsm,jsp,i,is
       integer lm,jm,jk,step,iflag,nj,mi
       integer mover(mmovers)
       real*8 phonon,gaussn,eps,xold(mdim,mmovers),pold,pnew,xm,delta
       real*8 ronf
       real*4 etime,tim(2)
       parameter (eps=1.d-32)
       external gaussn
       include 'pimc.cm'
       data iflag/0/
	
       if (iflag.eq.0) then
        iflag=1
        if (mod(nslices,nmovers).ne.0) then
         do lm=1,nslices
          nmovers=nmovers-1
          if (mod(nslices,nmovers).eq.0) goto 10
         enddo
10       write(*,*) 'number of movers modified to ',nmovers
        endif
        nrap=nslices/nmovers
       endif
c do one block of MC dynamics
       do step=1,nstep
c loops recursively over time slices
        do lm=1,nrap
         do nj=1,np
          js=lm
          do jm=1,nmovers
           mover(jm)=iwrap(js+(jm-1)*nrap)
c          write (*,*) nslices,nmovers,nrap,lm,jm,mover(jm)
          enddo
          do jm=1,nmovers
           js=mover(jm)
           jsm=iwrap(js-1)
           jsp=iwrap(js+1)
           sold(jm) = spring(js,nj)
           soldm(jm) = spring(jsm,nj)
           spring(js,nj)=0.d0
           spring(jsm,nj)=0.d0
           do l=1,ndim
            xold(l,jm) = x(l,nj,js)
c comput the trial free particle term
            xm=(x(l,nj,jsm)+x(l,nj,jsp))/2.d0
            x(l,nj,js)=xm+sigma*gaussn()
            spring(js,nj) = spring(js,nj)
     &                    +cke*(x(l,nj,js)-x(l,nj,jsp))**2
            spring(jsm,nj) = spring(jsm,nj)
     &                    +cke*(x(l,nj,js)-x(l,nj,jsm))**2
           enddo
          enddo
          delta=0.0
c computes the trial action
          if (alpha.ne.0.0) then
           pold=0.d0
           pnew=0.d0
           js=lm
           do jm=1,nmovers
            js=mover(jm)
c first comovers
            do jk=jm+1,nmovers
             is=mover(jk)
             il=abs(js-is)
c            write (*,*) js,is,il
             pold=pold+phonon(il,xold(1,jm),xold(1,jk),ndim)
             pnew=pnew+phonon(il,x(1,nj,js),x(1,nj,is),ndim)
            enddo
c then immobiles
            do jk=1,nmovers
             is=mover(jk)
             do l=1,nrap-1
              is=iwrap(is+1)
              il=abs(js-is)
c             write (*,*) js,is,il
              pold=pold+phonon(il,xold(1,jm),x(1,nj,is),ndim)
              pnew=pnew+phonon(il,x(1,nj,js),x(1,nj,is),ndim)
             enddo
            enddo
           enddo
           delta=-pnew+pold
          endif
          ntry=ntry+1
          call rmar(ronf)
          ronf=dlog(max(eps,ronf))
c         write (*,*) delta,ronf
	  if ( delta .ge. ronf ) then
c    ** .. accept **
           nacc = nacc + 1
c update action
           action = action-delta
          else
c    ** .. reject and replace **
           do jm=1,nmovers
            js=mover(jm)
            do l=1,ndim
             x(l,nj,js) = xold(l,jm)
            enddo
            spring(js,nj) = sold(jm)
            spring(jsm,nj) = soldm(jm)
           enddo
          endif
         enddo            ! over nj
        enddo        ! over nrap
c computing averages
        if (mod(step,nanal).eq.0) then
         tt=etime(tim)
         call anal(1,step)
         ttanal=ttanal+etime(tim)-tt
        endif
        if (mod(step,nspill).eq.0)
     &      write (31) (((x(l,nj,js),l=1,ndim),nj=1,np),js=1,nslices)
       enddo         ! over steps

       return
       end