subroutine PMFread
c Read the PMFs from wham
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.PMF'
      include 'COMMON.IOUNITS'
      include 'COMMON.MD'
#ifdef LANG0
#ifdef FIVEDIAG
      include 'COMMON.LANGEVIN.lang0.5diag'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
#else
      include 'COMMON.LANGEVIN'
#endif
      include 'COMMON.QRESTR'
      include 'COMMON.PMF'
      include 'COMMON.REMD'
      integer i,iumb,iiset,j,t,nbin
      double precision beta_h,qtemp,htemp
      read(inp,*) delta_q
      write (iout,*)"PMFread: nT",nrep," delta_q",delta_q," nW",nset
      print *,(remd_t(i),i=1,nrep)
      do j=1,nrep
        do iumb=1,nset
          read (inp,*) iiset,beta_h,nbin
          write (iout,*) iiset,beta_h,nbin
          if (iiset.ne.iumb) then
            write(iout,*) "Error: inconsistency in US windows",
     &         iiset,iumb 
            stop
          endif
          if (beta_h.ne.remd_t(j)) then
            write (iout,*) 
     &     "Error replica temperatures do not match PMF temperatures"
            write (iout,*) beta_h,remd_T(j)
            stop
          endif
          do i=1,nbin
            read (inp,*) qtemp,htemp
            t = int(qtemp/delta_q+1.0d-4)
            write (iout,*) qtemp,t,htemp
            if (i.eq.1) tmin(j,iumb)=t
            if (i.eq.nbin) tmax(j,iumb)=t
            PMFtab(t,j,iumb)=dlog(htemp)*Rb*remd_t(j)
          enddo ! i
        enddo ! iumb
      enddo ! j
      do iumb=1,nset
        write (iout,*)"Input PMFs, restraint",iumb," q0",qinfrag(1,iumb)
        write (iout,'(5x,20f10.1)') (remd_T(j),j=1,nrep)
        do i=0,int(1.0/delta_q)
          write (iout,'(f5.2,$)') i*delta_q
          do j=1,nrep
            if (i.lt.tmin(j,iumb).or.i.gt.tmax(j,iumb)) then
              write (iout,'("    ------",$)')
            else
              write (iout,'(f10.3$)') PMFtab(i,j,iumb)
            endif
          enddo
          write (iout,*) 
        enddo ! i
      enddo ! iumb
      return
      end
c------------------------------------------------------------------------
      subroutine PMF_energy(q,ePMF,ePMF_q)
c Calculate the energy and derivative in q due to the biasing PMF
c Caution! Only ONE q is handled, no multi-D q-restraints available!
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.PMF'
#ifdef MPI
      include 'mpif.h'
      include 'COMMON.SETUP'
      integer IERROR,ERRCODE
#endif
      include 'COMMON.IOUNITS'
      include 'COMMON.MD'
      include 'COMMON.QRESTR'
      include 'COMMON.REMD'
      include 'COMMON.PMF'
      integer i,iqmin,iqmax,irep
      double precision q,qmin,qmax,ePMF,ePMF_q
c Determine the location of the q
      irep=1
      do while (dabs(t_bath-remd_t(irep)).gt.1.0d-4)
        irep=irep+1
      enddo
      if (irep.gt.nrep) then
        write (iout,*) "CHUJ NASTAPIL in PMF_energy!!!"
        write (iout,*) 
     &    "Bath temperature does not match that of any replica"
        write (iout,*) "t_bath",t_bath
        write (iout,*) "remd_T",(remd_T(:nrep))
        write (iout,*) "Terminating..."
        call flush(iout)
#ifdef MPI
        call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
#endif
        stop 
      endif 
      iqmin=tmin(irep,iset)
      iqmax=tmax(irep,iset)
      qmin=iqmin*delta_q
      qmax=iqmax*delta_q
#ifdef DEBUG
      write (iout,*) "PMF_energy q",q," qmin",qmin," qmax",qmax,
     & " irep",irep," iset",iset," iqmin",iqmin," iqmax",iqmax
#endif
      if (q.le.qmin) then
        ePMF_q=(PMFtab(iqmin+1,irep,iset)-PMFtab(iqmin,irep,iset))
     &     /delta_q
        ePMF=PMFtab(iqmin,irep,iset)+ePMF_q*(q-qmin)
#ifdef DEBUG
       write (iout,*) "q<=qmin ePMF",ePMF," ePMF_q",ePMF_q
#endif
      else if (q.ge.qmax) then
        ePMF_q=(PMFtab(iqmax,irep,iset)-PMFtab(iqmax-1,irep,iset))/
     &      delta_q 
        ePMF=PMFtab(iqmax,irep,iset)+ePMF_q*(q-qmax)
#ifdef DEBUG
       write (iout,*) "q>=qmax ePMF",ePMF," ePMF_q",ePMF_q
#endif
      else
        do i=iqmin+1,iqmax
          qmax=i*delta_q
          if (q.ge.qmin .and. q.le.qmax) then
            ePMF_q=(PMFtab(i,irep,iset)-PMFtab(i-1,irep,iset))/delta_q
            ePMF=PMFtab(i-1,irep,iset)+ePMF_q*(q-qmin)
#ifdef DEBUG
       write (iout,*) "qmin<q<qmax bin",i," ePMF",ePMF," ePMF_q",ePMF_q
#endif
            exit
          endif
          qmin=qmax
        enddo
      endif
      return
      end