subroutine muca_delta(remd_t_bath,remd_ene,i,iex,delta)
      implicit real*8 (a-h,o-z)     
      include 'DIMENSIONS'
      include 'COMMON.MUCA'
      include 'COMMON.MD'
      double precision remd_t_bath(maxprocs)
      double precision remd_ene(maxprocs)
      double precision muca_ene
      double precision betai,betaiex,delta

      betai=1.0/(Rb*remd_t_bath(i))
      betaiex=1.0/(Rb*remd_t_bath(iex))
      
      delta=betai*(muca_ene(remd_ene(iex),i,remd_t_bath)-
     &                muca_ene(remd_ene(i),i,remd_t_bath))
     &          -betaiex*(muca_ene(remd_ene(iex),iex,remd_t_bath)-
     &                muca_ene(remd_ene(i),iex,remd_t_bath))

      return
      end
      
      double precision function muca_ene(energy,i,remd_t_bath)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.MUCA'
      include 'COMMON.MD'
      double precision y,yp,energy
      double precision remd_t_bath(maxprocs)
      integer i
 
      if (energy.lt.elowi(i)) then
        call splint(emuca,nemuca,nemuca2,nmuca,elowi(i),y,yp)
        muca_ene=remd_t_bath(i)*Rb*(yp*(energy-elowi(i))+y)
      elseif (energy.gt.ehighi(i)) then
        call splint(emuca,nemuca,nemuca2,nmuca,ehighi(i),y,yp)
        muca_ene=remd_t_bath(i)*Rb*(yp*(energy-ehighi(i))+y)
      else
        call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
        muca_ene=remd_t_bath(i)*Rb*y
      endif
      return
      end
      
      subroutine read_muca
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.MUCA'
      include 'COMMON.CONTROL'
      include 'COMMON.MD'
      include 'COMMON.REMD'
      include 'COMMON.SETUP'
      include 'COMMON.IOUNITS'
      double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
      imtime=0
      do i=1,4*maxres
        hist(i)=0
      enddo
      if (modecalc.eq.14.and..not.remd_tlist) then
                print *,"MUCAREMD works only with TLIST"
                stop
      endif
      open(89,file='muca.input')      
      read(89,*) 
      read(89,*) 
      if (modecalc.eq.14) then
        read(89,*) (elowi(i),ehighi(i),i=1,nrep)
       if (remd_mlist) then
        k=0
        do i=1,nrep
         do j=1,remd_m(i)
          i2rep(k)=i
          k=k+1
         enddo
        enddo
        elow=elowi(i2rep(me))
        ehigh=ehighi(i2rep(me))
        elowi(me+1)=elow
        ehighi(me+1)=ehigh        
       else 
        elow=elowi(me+1)
        ehigh=ehighi(me+1)
       endif
      else
        read(89,*) elow,ehigh
        elowi(1)=elow
        ehighi(1)=ehigh
      endif
      i=0
      do while(.true.)
       i=i+1
       read(89,*,end=100) emuca(i),nemuca(i)
cd      nemuca(i)=nemuca(i)*remd_t(me+1)*Rb
      enddo
 100  continue
      nmuca=i-1
      hbin=emuca(nmuca)-emuca(nmuca-1)
      write (iout,*) 'hbin',hbin      
      write (iout,*) me,'elow,ehigh',elow,ehigh
      yp1=0
      ypn=0
      call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
      factor_min=0.0d0
      factor_min=muca_factor(ehigh)
      call print_muca
      return
      end


      subroutine print_muca
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.MUCA'
      include 'COMMON.CONTROL'
      include 'COMMON.MD'
      include 'COMMON.REMD'
      include 'COMMON.SETUP'
      include 'COMMON.IOUNITS'
      double precision yp1,ypn,yp,x,muca_factor,y,muca_ene
      double precision dummy(maxprocs)

      if (remd_mlist) then
           k=0
           do i=1,nrep
            do j=1,remd_m(i)
             i2rep(k)=i
             k=k+1
            enddo
           enddo
      endif

      do i=1,nmuca
c       print *,'nemuca ',emuca(i),nemuca(i)
       do j=0,4
        x=emuca(i)+hbin/5*j
        if (modecalc.eq.14) then
         if (remd_mlist) then
          yp=muca_factor(x)*remd_t(i2rep(me))*Rb
          dummy(me+1)=remd_t(i2rep(me))
          y=muca_ene(x,me+1,dummy)
         else
          yp=muca_factor(x)*remd_t(me+1)*Rb
          y=muca_ene(x,me+1,remd_t)
         endif
         write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
     &      'muca factor ',x,yp,' muca ene',y
        else
         yp=muca_factor(x)*t_bath*Rb
         dummy(1)=t_bath
         y=muca_ene(x,1,dummy)
         write (iout,'(i4,i12,a12,2f15.5,a10,f15.5)') me,imtime,
     &      'muca factor ',x,yp,' muca ene',y         
        endif
       enddo
      enddo
      if(mucadyn.gt.0) then
       do i=1,nmuca
         write(iout,'(a13,i8,2f12.5)') 'nemuca after ',
     &             imtime,emuca(i),nemuca(i)
       enddo
      endif
      return
      end

      subroutine muca_update(energy)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.MUCA'
      include 'COMMON.CONTROL'
      include 'COMMON.MD'
      include 'COMMON.REMD'
      include 'COMMON.SETUP'
      include 'COMMON.IOUNITS'
      double precision energy
      double precision yp1,ypn
      integer k
      logical lnotend

      k=int((energy-emuca(1))/hbin)+1
      
      IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN
       if(energy.ge.ehigh) 
     &        write (iout,*) 'MUCA reject',energy,emuca(k)
       if(energy.ge.ehigh.and.(energy-ehigh).lt.hbin) then
         write (iout,*) 'MUCA ehigh',energy,emuca(k)
         do i=k,nmuca
           hist(i)=hist(i)+1
         enddo
       endif
       if(k.gt.0.and.energy.lt.ehigh) hist(k)=hist(k)+1       
      ELSE
       if(k.gt.0.and.k.lt.4*maxres) hist(k)=hist(k)+1       
      ENDIF
      if(mod(imtime,mucadyn).eq.0) then

         do i=1,nmuca
          IF(muca_smooth.eq.2.or.muca_smooth.eq.3) THEN
           nemuca(i)=nemuca(i)+dlog(hist(i)+1)
          ELSE
           if (hist(i).gt.0) hist(i)=dlog(hist(i))         
           nemuca(i)=nemuca(i)+hist(i)
          ENDIF
          hist(i)=0
          write(iout,'(a24,i8,2f12.5)')'nemuca before smoothing ',
     &          imtime,emuca(i),nemuca(i)
         enddo


         lnotend=.true.
         ismooth=0
         ist=2
         ien=nmuca-1
        IF(muca_smooth.eq.1.or.muca_smooth.eq.3) THEN         
c         lnotend=.false.         
c         do i=1,nmuca-1
c           do j=i+1,nmuca
c            if(nemuca(j).lt.nemuca(i)) lnotend=.true.
c           enddo
c         enddo         
         do while(lnotend)
          ismooth=ismooth+1
          write (iout,*) 'MUCA update smoothing',ist,ien
          do i=ist,ien
           nemuca(i)=(nemuca(i-1)+nemuca(i)+nemuca(i+1))/3
          enddo
          lnotend=.false.
          ist=0
          ien=0
          do i=1,nmuca-1
           do j=i+1,nmuca
            if(nemuca(j).lt.nemuca(i)) then 
              lnotend=.true.
              if(ist.eq.0) ist=i-1
              if(ien.lt.j+1) ien=j+1
            endif
           enddo
          enddo
         enddo
        ENDIF 

         write (iout,*) 'MUCA update ',imtime,' smooth= ',ismooth
         yp1=0
         ypn=0
         call spline(emuca,nemuca,nmuca,yp1,ypn,nemuca2)
         call print_muca
         
      endif
      return
      end
      
      double precision function muca_factor(energy)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.MUCA'
      double precision y,yp,energy
      
      if (energy.lt.elow) then
        call splint(emuca,nemuca,nemuca2,nmuca,elow,y,yp)
      elseif (energy.gt.ehigh) then
        call splint(emuca,nemuca,nemuca2,nmuca,ehigh,y,yp)
      else
        call splint(emuca,nemuca,nemuca2,nmuca,energy,y,yp)
      endif
      
      if(yp.ge.factor_min) then
       muca_factor=yp
      else
       muca_factor=factor_min
      endif
cd      print *,'energy, muca_factor',energy,muca_factor
      return
      end
      

      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL*8 yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      REAL*8 p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then 
      y2(1)=0. 
      u(1)=0.
      else 
         y2(1)=-0.5
         u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do i=2,n-1 
         sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
         p=sig*y2(i-1)+2.
         y2(i)=(sig-1.)/p
         u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
     *        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
      enddo 
      if (ypn.gt..99e30) then 
         qn=0.
         un=0.
      else 
         qn=0.5
         un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do k=n-1,1,-1 
         y2(k)=y2(k)*y2(k+1)+u(k) 
      enddo 
      return
      END 


      SUBROUTINE splint(xa,ya,y2a,n,x,y,yp)
      INTEGER n
      REAL*8 x,y,xa(n),y2a(n),ya(n),yp
      INTEGER k,khi,klo
      REAL*8 a,b,h
      klo=1 
      khi=n
 1    if (khi-klo.gt.1) then
         k=(khi+klo)/2
         if (xa(k).gt.x) then
            khi=k
         else
            klo=k
         endif
         goto 1
      endif 
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint' 
      a=(xa(khi)-x)/h 
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     *     ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
      yp=-ya(klo)/h+ya(khi)/h-3*(a**2)*y2a(klo)*h/6.
     +     +(3*(b**2)-1)*y2a(khi)*h/6.
      return
      END