subroutine WHAM_CALC(islice,*)
! Weighed Histogram Analysis Method (WHAM) code
! Written by A. Liwo based on the work of Kumar et al., 
! J.Comput.Chem., 13, 1011 (1992)
!
! 2/1/05 Multiple temperatures allowed.
! 2/2/05 Free energies calculated directly from data points
!  acc. to Eq. (21) of Kumar et al.; final histograms also
!  constructed based on this equation.
! 2/12/05 Multiple parameter sets included
!
! 2/2/05 Parallel version
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "DIMENSIONS.FREE"
      integer MaxBinRms,MaxBinRgy
      parameter (MaxBinRms=100,MaxBinRgy=100)
      integer MaxHdim
c      parameter (MaxHdim=200000)
      parameter (MaxHdim=200)
      integer maxinde
      parameter (maxinde=200)
#ifdef MPI
      include "mpif.h"
      include "COMMON.MPI"
      integer ierror,errcode,status(MPI_STATUS_SIZE) 
#endif
      include "COMMON.CONTROL"
      include "COMMON.IOUNITS"
      include "COMMON.FREE"
      include "COMMON.ENERGIES"
      include "COMMON.FFIELD"
      include "COMMON.SBRIDGE"
      include "COMMON.PROT"
      include "COMMON.ENEPS"
      integer MaxPoint,MaxPointProc
      parameter (MaxPoint=MaxStr,
     & MaxPointProc=MaxStr_Proc)
      double precision finorm_max,potfac,entmin,entmax,expfac,vf
      double precision entfac_min,entfac_min_t
      parameter (finorm_max=1.0d0)
      integer islice
      integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
      integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
     &  nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
      integer htot(0:MaxHdim),histent(0:2000)
      double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
      double precision energia(0:max_ene)
#ifdef MPI
      integer tmax_t,upindE_p
      double precision fi_p(MaxR,MaxT_h,Max_Parm),
     &  fi_p_min(MaxR,MaxT_h,Max_Parm)
      double precision sumW_p(0:Max_GridT,Max_Parm),
     & sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm),
     & sumQ_p(MaxQ1,0:Max_GridT,Max_Parm),
     & sumQsq_p(MaxQ1,0:Max_GridT,Max_Parm),
     & sumEQ_p(MaxQ1,0:Max_GridT,Max_Parm),
     & sumEprim_p(MaxQ1,0:Max_GridT,Max_Parm), 
     & sumEbis_p(0:Max_GridT,Max_Parm)
      double precision hfin_p(0:MaxHdim,maxT_h),
     & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
     & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
      double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
      double precision potEmin_t_all(maxT_h,Max_Parm),entmin_p,entmax_p
      integer histent_p(0:2000)
      logical lprint /.true./
#endif
      double precision rgymin,rmsmin,rgymax,rmsmax
      double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
     & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
     & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
     & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
     & weight,econstr
      double precision fi(MaxR,maxT_h,Max_Parm),
     & fi_min(MaxR,maxT_h,Max_Parm),
     & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
     & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
     & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
     & potEmin_all(maxT_h,Max_Parm),potEmin,potEmin_min,ent,
     & hfin_ent(0:MaxHdim),vmax,aux
      double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
     &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,
     &  eplus,eminus,logfac,tanhT,tt
      double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
     &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
     &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor

      integer ind_point(maxpoint),upindE,indE
      character*16 plik
      character*1 licz1
      character*2 licz2
      character*3 licz3
      character*128 nazwa
      character*30 frm_write
      integer ilen
      external ilen
      write(licz2,'(bz,i2.2)') islice
      nbin1 = 1.0d0/delta
      write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
     &  i2/80(1h-)//)') islice
      write (iout,*) "delta",delta," nbin1",nbin1
      write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
      call flush(iout)
      dmin=0.0d0
      tmax=0
      do i=1,nParmset
        do j=1,nT_h(i)
          potEmin_all(j,i)=1.0d10
        enddo
      enddo
      rgymin=1.0d10
      rmsmin=1.0d10
      rgymax=0.0d0
      rmsmax=0.0d0
      do t=0,MaxN
        htot(t)=0
      enddo
#ifdef MPI
      do i=1,scount(me1)
#else
      do i=1,ntot(islice)
#endif
        if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
        if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
        if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
        if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
        ind_point(i)=0
        do j=nQ,1,-1
          ind=(q(j,i)-dmin+1.0d-8)/delta
          if (j.eq.1) then
            ind_point(i)=ind_point(i)+ind
          else 
            ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
          endif
          if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
            write (iout,*) "Error - index exceeds range for point",i,
     &      " q=",q(j,i)," ind",ind_point(i)
#ifdef MPI 
            write (iout,*) "Processor",me1
            call flush(iout)
            call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
#endif
            stop
          endif
        enddo ! j
        if (ind_point(i).gt.tmax) tmax=ind_point(i)
        htot(ind_point(i))=htot(ind_point(i))+1
#ifdef DEBUG
        write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
     &   " htot",htot(ind_point(i))
        call flush(iout)
#endif
      enddo ! i
      call flush(iout)

      nbin=nbin1**nQ-1
      write (iout,'(a)') "Numbers of counts in Q bins"
      do t=0,tmax
        if (htot(t).gt.0) then
        write (iout,'(i15,$)') t
        liczba=t
        do j=1,nQ
          jj = mod(liczba,nbin1)
          liczba=liczba/nbin1
          write (iout,'(i5,$)') jj
        enddo
        write (iout,'(i8)') htot(t)
        endif
      enddo
      do iparm=1,nParmSet
      write (iout,'(a,i3)') "Number of data points for parameter set",
     & iparm
      write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
     & ib=1,nT_h(iparm))
      write (iout,'(i8)') stot(islice)
      write (iout,'(a)')
      enddo
      call flush(iout)

#ifdef MPI
      call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
     &  WHAM_COMM,IERROR)
      tmax=tmax_t
      call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
     &  MPI_MIN,WHAM_COMM,IERROR)
      call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
     &  MPI_MAX,WHAM_COMM,IERROR)
      call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
     &  MPI_MIN,WHAM_COMM,IERROR)
      call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
     &  MPI_MAX,WHAM_COMM,IERROR)
      rgymin=rgymin_t
      rgymax=rgymax_t
      rmsmin=rmsmin_t
      rmsmax=rmsmax_t
#endif
      rmsmin=deltrms*dint(rmsmin/deltrms)
      rmsmax=deltrms*dint(rmsmax/deltrms)
      rgymin=deltrms*dint(rgymin/deltrgy)
      rgymax=deltrms*dint(rgymax/deltrgy)
      nbin_rms=(rmsmax-rmsmin)/deltrms
      nbin_rgy=(rgymax-rgymin)/deltrgy
      write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
     & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
      nFi=0
      do i=1,nParmSet
        do j=1,nT_h(i)
          nFi=nFi+nR(j,i)
        enddo
      enddo
      write (iout,*) "nFi",nFi
! Compute the Boltzmann factor corresponing to restrain potentials in different
! simulations.
#ifdef MPI
      do i=1,scount(me1)
#else
      do i=1,ntot(islice)
#endif
c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
        do iparm=1,nParmSet
#ifdef DEBUG
          write (iout,'(2i5,21f8.2)') i,iparm,
     &     (enetb(k,i,iparm),k=1,21)
#endif
          call restore_parm(iparm)
#ifdef DEBUG
          write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
     &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
     &      wtor_d,wsccor,wbond
#endif
          do ib=1,nT_h(iparm)
            if (rescale_mode.eq.1) then
              quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
              quotl=1.0d0
              kfacl=1.0d0
              do l=1,5
                quotl1=quotl
                quotl=quotl*quot
                kfacl=kfacl*kfac
                fT(l)=kfacl/(kfacl-1.0d0+quotl)
              enddo
#if defined(FUNCTH)
              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
#elif defined(FUNCT)
              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
#else
              ft(6)=1.0d0
#endif
            else if (rescale_mode.eq.2) then
              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
              quotl=1.0d0
              do l=1,5
                quotl=quotl*quot
                fT(l)=1.12692801104297249644d0/
     &             dlog(dexp(quotl)+dexp(-quotl))
              enddo
#if defined(FUNCTH)
              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
#elif defined(FUNCT)
              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
#else
              ft(6)=1.0d0
#endif
c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
            else if (rescale_mode.eq.0) then
              do l=1,6
                fT(l)=1.0d0
              enddo
            else
              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
     &          rescale_mode
              call flush(iout)
              return1
            endif
            evdw=enetb(1,i,iparm)
            evdw_t=enetb(21,i,iparm)
#ifdef SCP14
            evdw2_14=enetb(17,i,iparm)
            evdw2=enetb(2,i,iparm)+evdw2_14
#else
            evdw2=enetb(2,i,iparm)
            evdw2_14=0.0d0
#endif
#ifdef SPLITELE
            ees=enetb(3,i,iparm)
            evdw1=enetb(16,i,iparm)
#else
            ees=enetb(3,i,iparm)
            evdw1=0.0d0
#endif
            ecorr=enetb(4,i,iparm)
            ecorr5=enetb(5,i,iparm)
            ecorr6=enetb(6,i,iparm)
            eel_loc=enetb(7,i,iparm)
            eello_turn3=enetb(8,i,iparm)
            eello_turn4=enetb(9,i,iparm)
            eturn6=enetb(10,i,iparm)
            ebe=enetb(11,i,iparm)
            escloc=enetb(12,i,iparm)
            etors=enetb(13,i,iparm)
            etors_d=enetb(14,i,iparm)
            ehpb=enetb(15,i,iparm)
            estr=enetb(18,i,iparm)
            esccor=enetb(19,i,iparm)
            edihcnstr=enetb(20,i,iparm)
#ifdef DEBUG
            write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
     &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
     &       etors,etors_d,eello_turn3,eello_turn4,esccor
#endif

#ifdef SPLITELE
            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
     &      +wvdwpp*evdw1
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
     &      +ft(2)*wturn3*eello_turn3
     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
     &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
     &      +wbond*estr
#else
            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
     &      +ft(1)*welec*(ees+evdw1)
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
     &      +ft(2)*wturn3*eello_turn3
     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
     &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
     &      +wbond*estr
#endif
#ifdef DEBUG
            write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
     &        etot
#endif
#ifdef DEBUG
            if (iparm.eq.1 .and. ib.eq.1) then
            write (iout,*)"Conformation",i
            energia(0)=etot
            do k=1,max_ene
              energia(k)=enetb(k,i,iparm)
            enddo
            call enerprint(energia(0),fT)
            endif
#endif
            do kk=1,nR(ib,iparm)
              Econstr=0.0d0
              do j=1,nQ
                dd = q(j,i)
                Econstr=Econstr+Kh(j,kk,ib,iparm)
     &           *(dd-q0(j,kk,ib,iparm))**2
              enddo
              v(i,kk,ib,iparm)=
     &          -beta_h(ib,iparm)*(etot+Econstr)
#ifdef DEBUG
              write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
     &         etot,v(i,kk,ib,iparm)
#endif
            enddo ! kk
          enddo   ! ib
        enddo     ! iparm
      enddo       ! i
! Simple iteration to calculate free energies corresponding to all simulation
! runs.
      do iter=1,maxit 
        
! Compute new free-energy values corresponding to the righ-hand side of the 
! equation and their derivatives.
        write (iout,*) "------------------------fi"
        entfac_min=1.0d10
#ifdef MPI
        do t=1,scount(me1)
#else
        do t=1,ntot(islice)
#endif
          vmax=-1.0d+20
          do i=1,nParmSet
            do k=1,nT_h(i)
              do l=1,nR(k,i)
                vf=v(t,l,k,i)+f(l,k,i)
                if (vf.gt.vmax) vmax=vf
              enddo
            enddo
          enddo        
          denom=0.0d0
          do i=1,nParmSet
            do k=1,nT_h(i)
              do l=1,nR(k,i)
                aux=f(l,k,i)+v(t,l,k,i)-vmax
                if (aux.gt.-200.0d0)
     &          denom=denom+snk(l,k,i,islice)*dexp(aux)
              enddo
            enddo
          enddo
          entfac(t)=-dlog(denom)-vmax
          if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
#ifdef DEBUG
          write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
#endif
        enddo
c#ifdef MPI
c        write (iout,*) "entfac_min before AllReduce",entfac_min
c        call MPI_AllReduce(entfac_min,entfac_min_t,1,
c     &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
c        entfac_min=entfac_min_t
c        write (iout,*) "entfac_min after AllReduce",entfac_min
c#endif
c#ifdef MPI
c        do t=1,scount(me)
c          entfac(t)=entfac(t)-entfac_min
c        enddo
c#else
c        do t=1,ntot(islice)
c          entfac(t)=entfac(t)-entfac_min
c        enddo
c#endif
        do iparm=1,nParmSet
          do iib=1,nT_h(iparm)
            do ii=1,nR(iib,iparm)
#ifdef MPI
              fi_p_min(ii,iib,iparm)=-1.0d10
              do t=1,scount(me)
                aux=v(t,ii,iib,iparm)+entfac(t)
                if (aux.gt.fi_p_min(ii,iib,iparm))
     &                                   fi_p_min(ii,iib,iparm)=aux
              enddo
#else
              do t=1,ntot(islice)
                aux=v(t,ii,iib,iparm)+entfac(t)
                if (aux.gt.fi_min(ii,iib,iparm)) 
     &                                     fi_min(ii,iib,iparm)=aux
              enddo
#endif
            enddo ! ii
          enddo ! iib
        enddo ! iparm
#ifdef MPI
#ifdef DEBUG
        write (iout,*) "fi_min before AllReduce"
        do i=1,nParmSet
          do j=1,nT_h(i)
            write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
          enddo
        enddo
#endif
        call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
     &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
#ifdef DEBUG
        write (iout,*) "fi_min after AllReduce"
        do i=1,nParmSet
          do j=1,nT_h(i)
            write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
          enddo
        enddo
#endif
#endif
        do iparm=1,nParmSet
          do iib=1,nT_h(iparm)
            do ii=1,nR(iib,iparm)
#ifdef MPI
              fi_p(ii,iib,iparm)=0.0d0
              do t=1,scount(me)
                fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
     &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
#ifdef DEBUG
              write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
     &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
     &         fi_p(ii,iib,iparm)
#endif
              enddo
#else
              fi(ii,iib,iparm)=0.0d0
              do t=1,ntot(islice)
                fi(ii,iib,iparm)=fi(ii,iib,iparm)
     &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
              enddo
#endif
            enddo ! ii
          enddo ! iib
        enddo ! iparm

#ifdef MPI
#ifdef DEBUG
        write (iout,*) "fi before MPI_Reduce me",me,' master',master
        do iparm=1,nParmSet
          do ib=1,nT_h(iparm)
            write (iout,*) "iparm",iparm," ib",ib
            write (iout,*) "beta=",beta_h(ib,iparm)
            write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
          enddo
        enddo
#endif
#ifdef DEBUG
        write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
     &   maxR*MaxT_h*nParmSet
        write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
     &   " WHAM_COMM",WHAM_COMM
#endif
        call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
     &   MPI_DOUBLE_PRECISION,
     &   MPI_SUM,Master,WHAM_COMM,IERROR)
#ifdef DEBUG
        write (iout,*) "fi after MPI_Reduce nparmset",nparmset
        do iparm=1,nParmSet
          write (iout,*) "iparm",iparm
          do ib=1,nT_h(iparm)
            write (iout,*) "beta=",beta_h(ib,iparm)
            write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
          enddo
        enddo
#endif
        if (me1.eq.Master) then
#endif
        avefi=0.0d0
        do iparm=1,nParmSet
          do ib=1,nT_h(iparm)
            do i=1,nR(ib,iparm)
              fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
              avefi=avefi+fi(i,ib,iparm)
            enddo
          enddo
        enddo
        avefi=avefi/nFi
        do iparm=1,nParmSet
          write (iout,*) "Parameter set",iparm
          do ib =1,nT_h(iparm)
            write (iout,*) "beta=",beta_h(ib,iparm)
            do i=1,nR(ib,iparm)
              fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
            enddo
            write (iout,'(6f15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
            write (iout,'(6f15.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
          enddo
        enddo

! Compute the norm of free-energy increments.
        finorm=0.0d0
        do iparm=1,nParmSet
          do ib=1,nT_h(iparm)
            do i=1,nR(ib,iparm)
              finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
              f(i,ib,iparm)=fi(i,ib,iparm)
            enddo  
          enddo
        enddo

        write (iout,*) 'Iteration',iter,' finorm',finorm

#ifdef MPI
        endif
        call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
     &   MPI_DOUBLE_PRECISION,Master,
     &   WHAM_COMM,IERROR)
        call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
     &   WHAM_COMM,IERROR)
#endif
! Exit, if the increment norm is smaller than pre-assigned tolerance.
        if (finorm.lt.fimin) then
          write (iout,*) 'Iteration converged'
          goto 20
        endif

      enddo ! iter

   20 continue
! Now, put together the histograms from all simulations, in order to get the
! unbiased total histogram.

C Determine the minimum free energies
#ifdef MPI
      do i=1,scount(me1)
#else
      do i=1,ntot(islice)
#endif
c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
        do iparm=1,nParmSet
#ifdef DEBUG
          write (iout,'(2i5,21f8.2)') i,iparm,
     &     (enetb(k,i,iparm),k=1,21)
#endif
          call restore_parm(iparm)
#ifdef DEBUG
          write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
     &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
     &      wtor_d,wsccor,wbond
#endif
          do ib=1,nT_h(iparm)
            if (rescale_mode.eq.1) then
              quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
              quotl=1.0d0
              kfacl=1.0d0
              do l=1,5
                quotl1=quotl
                quotl=quotl*quot
                kfacl=kfacl*kfac
                fT(l)=kfacl/(kfacl-1.0d0+quotl)
              enddo
#if defined(FUNCTH)
              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
#elif defined(FUNCT)
              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
#else
              ft(6)=1.0d0
#endif
            else if (rescale_mode.eq.2) then
              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
              quotl=1.0d0
              do l=1,5
                quotl=quotl*quot
                fT(l)=1.12692801104297249644d0/
     &             dlog(dexp(quotl)+dexp(-quotl))
              enddo
#if defined(FUNCTH)
              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
#elif defined(FUNCT)
              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
#else
              ft(6)=1.0d0
#endif
c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
            else if (rescale_mode.eq.0) then
              do l=1,6
                fT(l)=1.0d0
              enddo
            else
              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
     &          rescale_mode
              call flush(iout)
              return1
            endif
            evdw=enetb(1,i,iparm)
            evdw_t=enetb(21,i,iparm)
#ifdef SCP14
            evdw2_14=enetb(17,i,iparm)
            evdw2=enetb(2,i,iparm)+evdw2_14
#else
            evdw2=enetb(2,i,iparm)
            evdw2_14=0.0d0
#endif
#ifdef SPLITELE
            ees=enetb(3,i,iparm)
            evdw1=enetb(16,i,iparm)
#else
            ees=enetb(3,i,iparm)
            evdw1=0.0d0
#endif
            ecorr=enetb(4,i,iparm)
            ecorr5=enetb(5,i,iparm)
            ecorr6=enetb(6,i,iparm)
            eel_loc=enetb(7,i,iparm)
            eello_turn3=enetb(8,i,iparm)
            eello_turn4=enetb(9,i,iparm)
            eturn6=enetb(10,i,iparm)
            ebe=enetb(11,i,iparm)
            escloc=enetb(12,i,iparm)
            etors=enetb(13,i,iparm)
            etors_d=enetb(14,i,iparm)
            ehpb=enetb(15,i,iparm)
            estr=enetb(18,i,iparm)
            esccor=enetb(19,i,iparm)
            edihcnstr=enetb(20,i,iparm)
#ifdef DEBUG
            write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
     &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
     &       etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr
#endif

#ifdef SPLITELE
            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
     &      +wvdwpp*evdw1
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
     &      +ft(2)*wturn3*eello_turn3
     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
     &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
     &      +wbond*estr
#else
            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
     &      +ft(1)*welec*(ees+evdw1)
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
     &      +ft(2)*wturn3*eello_turn3
     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
     &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
     &      +wbond*estr
#endif
c            write (iout,*) "i",i," ib",ib,
c     &      " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
c     &      " entfac",entfac(i)
            etot=etot-entfac(i)/beta_h(ib,iparm)
            if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
c            write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
          enddo   ! ib
        enddo     ! iparm
      enddo       ! i
#ifdef DEBUG
      write (iout,*) "The potEmin array before reduction"
      do i=1,nParmSet
        write (iout,*) "Parameter set",i
        do j=1,nT_h(i)
          write (iout,*) j,PotEmin_all(j,i)
        enddo
      enddo
      write (iout,*) "potEmin_min",potEmin_min
#endif
#ifdef MPI
C Determine the minimum energes for all parameter sets and temperatures
      call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
     &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
      do i=1,nParmSet
        do j=1,nT_h(i)
          potEmin_all(j,i)=potEmin_t_all(j,i)
        enddo
      enddo
#endif
      potEmin_min=potEmin_all(1,1)
      do i=1,nParmSet
        do j=1,nT_h(i)
          if (potEmin_all(j,i).lt.potEmin_min) 
     &             potEmin_min=potEmin_all(j,i)
        enddo
      enddo
#ifdef DEBUG
      write (iout,*) "The potEmin array"
      do i=1,nParmSet
        write (iout,*) "Parameter set",i
        do j=1,nT_h(i)
          write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
     &      PotEmin_all(j,i)
        enddo
      enddo
      write (iout,*) "potEmin_min",potEmin_min
#endif

#ifdef MPI
      do t=0,tmax
        hfin_ent_p(t)=0.0d0
      enddo
#else
      do t=0,tmax
        hfin_ent(t)=0.0d0
      enddo
#endif
      write (iout,*) "--------------hist"
#ifdef MPI
      do iparm=1,nParmSet
        do i=0,nGridT
          sumW_p(i,iparm)=0.0d0
          sumE_p(i,iparm)=0.0d0
          sumEbis_p(i,iparm)=0.0d0
          sumEsq_p(i,iparm)=0.0d0
          do j=1,nQ+2
            sumQ_p(j,i,iparm)=0.0d0
            sumQsq_p(j,i,iparm)=0.0d0
            sumEQ_p(j,i,iparm)=0.0d0
          enddo
        enddo
      enddo
      upindE_p=0
#else
      do iparm=1,nParmSet
        do i=0,nGridT
          sumW(i,iparm)=0.0d0
          sumE(i,iparm)=0.0d0
          sumEbis(i,iparm)=0.0d0
          sumEsq(i,iparm)=0.0d0
          do j=1,nQ+2
            sumQ(j,i,iparm)=0.0d0
            sumQsq(j,i,iparm)=0.0d0
            sumEQ(j,i,iparm)=0.0d0
          enddo
        enddo
      enddo
      upindE=0
#endif
c 8/26/05 entropy distribution
#ifdef MPI
      entmin_p=1.0d10
      entmax_p=-1.0d10
      do t=1,scount(me1)
c        ent=-dlog(entfac(t))
        ent=entfac(t)
        if (ent.lt.entmin_p) entmin_p=ent
        if (ent.gt.entmax_p) entmax_p=ent
      enddo
      write (iout,*) "entmin",entmin_p," entmax",entmax_p
      call flush(iout)
      call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
     &  WHAM_COMM,IERROR)
      call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
     &  WHAM_COMM,IERROR)
      ientmax=entmax-entmin 
      if (ientmax.gt.2000) ientmax=2000
      write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
      call flush(iout)
      do t=1,scount(me1)
c        ient=-dlog(entfac(t))-entmin
        ient=entfac(t)-entmin
        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
      enddo
      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
     &  MPI_SUM,WHAM_COMM,IERROR)
      if (me1.eq.Master) then
        write (iout,*) "Entropy histogram"
        do i=0,ientmax
          write(iout,'(f15.4,i10)') entmin+i,histent(i)
        enddo
      endif
#else
      entmin=1.0d10
      entmax=-1.0d10
      do t=1,ntot(islice)
        ent=entfac(t)
        if (ent.lt.entmin) entmin=ent
        if (ent.gt.entmax) entmax=ent
      enddo
      ientmax=-dlog(entmax)-entmin
      if (ientmax.gt.2000) ientmax=2000
      do t=1,ntot(islice)
        ient=entfac(t)-entmin
        if (ient.le.2000) histent(ient)=histent(ient)+1
      enddo
      write (iout,*) "Entropy histogram"
      do i=0,ientmax
        write(iout,'(2f15.4)') entmin+i,histent(i)
      enddo
#endif
      
#ifdef MPI
c      write (iout,*) "me1",me1," scount",scount(me1)

      do iparm=1,nParmSet

#ifdef MPI
        do ib=1,nT_h(iparm)
          do t=0,tmax
            hfin_p(t,ib)=0.0d0
          enddo
        enddo
        do i=1,maxindE
          histE_p(i)=0.0d0
        enddo
#else
        do ib=1,nT_h(iparm)
          do t=0,tmax
            hfin(t,ib)=0.0d0
          enddo
        enddo
        do i=1,maxindE
          histE(i)=0.0d0
        enddo
#endif
        do ib=1,nT_h(iparm)
          do i=0,MaxBinRms
            do j=0,MaxBinRgy
              hrmsrgy(j,i,ib)=0.0d0
#ifdef MPI
              hrmsrgy_p(j,i,ib)=0.0d0
#endif
            enddo
          enddo
        enddo

        do t=1,scount(me1)
#else
        do t=1,ntot(islice)
#endif
          ind = ind_point(t)
#ifdef MPI
          hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
#else
          hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
#endif
          call restore_parm(iparm)
          evdw=enetb(21,t,iparm)
          evdw_t=enetb(1,t,iparm)
#ifdef SCP14
          evdw2_14=enetb(17,t,iparm)
          evdw2=enetb(2,t,iparm)+evdw2_14
#else
          evdw2=enetb(2,t,iparm)
          evdw2_14=0.0d0
#endif
#ifdef SPLITELE
          ees=enetb(3,t,iparm)
          evdw1=enetb(16,t,iparm)
#else
          ees=enetb(3,t,iparm)
          evdw1=0.0d0
#endif
          ecorr=enetb(4,t,iparm)
          ecorr5=enetb(5,t,iparm)
          ecorr6=enetb(6,t,iparm)
          eel_loc=enetb(7,t,iparm)
          eello_turn3=enetb(8,t,iparm)
          eello_turn4=enetb(9,t,iparm)
          eturn6=enetb(10,t,iparm)
          ebe=enetb(11,t,iparm)
          escloc=enetb(12,t,iparm)
          etors=enetb(13,t,iparm)
          etors_d=enetb(14,t,iparm)
          ehpb=enetb(15,t,iparm)
          estr=enetb(18,t,iparm)
          esccor=enetb(19,t,iparm)
          edihcnstr=enetb(20,t,iparm)
          do k=0,nGridT
            betaT=startGridT+k*delta_T
            temper=betaT
c            fT=T0/betaT
c            ft=2*T0/(T0+betaT)
            if (rescale_mode.eq.1) then
              quot=betaT/T0
              quotl=1.0d0
              kfacl=1.0d0
              do l=1,5
                quotl1=quotl
                quotl=quotl*quot
                kfacl=kfacl*kfac
                denom=kfacl-1.0d0+quotl
                fT(l)=kfacl/denom
                ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
                ftbis(l)=l*kfacl*quotl1*
     &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
              enddo
#if defined(FUNCTH)
              ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
     &                  320.0d0
              ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
     &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
#elif defined(FUNCT)
              fT(6)=betaT/T0
              ftprim(6)=1.0d0/T0
              ftbis(6)=0.0d0
#else
              fT(6)=1.0d0
              ftprim(6)=0.0d0
              ftbis(6)=0.0d0
#endif
            else if (rescale_mode.eq.2) then
              quot=betaT/T0
              quotl=1.0d0
              do l=1,5
                quotl1=quotl
                quotl=quotl*quot
                eplus=dexp(quotl)
                eminus=dexp(-quotl)
                logfac=1.0d0/dlog(eplus+eminus)
                tanhT=(eplus-eminus)/(eplus+eminus)
                fT(l)=1.12692801104297249644d0*logfac
                ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
                ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
     &          2*l*quotl1/T0*logfac*
     &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
     &          +ftprim(l)*tanhT)
              enddo
#if defined(FUNCTH)
              ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
     &                 320.0d0
              ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
     &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
#elif defined(FUNCT)
              fT(6)=betaT/T0
              ftprim(6)=1.0d0/T0
              ftbis(6)=0.0d0
#else
              fT(6)=1.0d0
              ftprim(6)=0.0d0
              ftbis(6)=0.0d0
#endif
            else if (rescale_mode.eq.0) then
              do l=1,5
                fT(l)=1.0d0
                ftprim(l)=0.0d0
              enddo
            else
              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
     &          rescale_mode
              call flush(iout)
              return1
            endif
c            write (iout,*) "ftprim",ftprim
c            write (iout,*) "ftbis",ftbis
            betaT=1.0d0/(1.987D-3*betaT)
            if (betaT.ge.beta_h(1,iparm)) then
              potEmin=potEmin_all(1,iparm)+
     &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/ 
     &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
     &         (1.0/betaT-1.0/beta_h(1,iparm))
#ifdef DEBUG
              write(iout,*) "first",temper,potEmin
#endif
            else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
              potEmin=potEmin_all(nT_h(iparm),iparm)+
     &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
     &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
     &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
#ifdef DEBUG
              write (iout,*) "last",temper,potEmin
#endif
            else
              do l=1,nT_h(iparm)-1
                if (betaT.le.beta_h(l,iparm) .and.
     &              betaT.gt.beta_h(l+1,iparm)) then
                  potEmin=potEmin_all(l,iparm)
#ifdef DEBUG
                  write (iout,*) "l",l,
     &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
     &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
#endif
                  exit
                endif
              enddo
            endif
#ifdef DEBUG
            write (iout,*) "k",k," potEmin",potEmin
#endif
#ifdef SPLITELE
            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
     &      +wvdwpp*evdw1
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
     &      +ft(2)*wturn3*eello_turn3
     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
     &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
     &      +wbond*estr
            eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
     &            +ftprim(1)*wtor*etors+
     &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
     &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
     &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
     &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
     &            ftprim(1)*wsccor*esccor
            ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
     &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
     &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
     &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
     &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
     &            ftbis(1)*wsccor*esccor
#else
            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
     &      +ft(1)*welec*(ees+evdw1)
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
     &      +ft(2)*wturn3*eello_turn3
     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
     &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
     &      +wbond*estr
            eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
     &           +ftprim(1)*wtor*etors+
     &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
     &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
     &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
     &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
     &            ftprim(1)*wsccor*esccor
            ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
     &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
     &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
     &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
     &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
     &            ftprim(1)*wsccor*esccor
#endif
            weight=dexp(-betaT*(etot-potEmin)+entfac(t))
#ifdef DEBUG
            write (iout,*) "iparm",iparm," t",t," temper",temper,
     &        " etot",etot," entfac",entfac(t),
     &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,
     &        " boltz",-betaT*(etot-potEmin)+entfac(t),
     &        " weight",weight," ebis",ebis
#endif
            etot=etot-temper*eprim
#ifdef MPI
            sumW_p(k,iparm)=sumW_p(k,iparm)+weight
            sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
            sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
            sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
            do j=1,nQ+2
              sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
              sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
              sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
     &         +etot*q(j,t)*weight
            enddo
#else
            sumW(k,iparm)=sumW(k,iparm)+weight
            sumE(k,iparm)=sumE(k,iparm)+etot*weight
            sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
            sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
            do j=1,nQ+2
              sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
              sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
              sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
     &         +etot*q(j,t)*weight
            enddo
#endif
          enddo
          indE = aint(potE(t,iparm)-aint(potEmin))
          if (indE.ge.0 .and. indE.le.maxinde) then
            if (indE.gt.upindE_p) upindE_p=indE
            histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
          endif
#ifdef MPI
          do ib=1,nT_h(iparm)
            potEmin=potEmin_all(ib,iparm)
            expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
            hfin_p(ind,ib)=hfin_p(ind,ib)+
     &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
             if (rmsrgymap) then
               indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) 
               indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
               hrmsrgy_p(indrgy,indrms,ib)=
     &           hrmsrgy_p(indrgy,indrms,ib)+expfac
             endif
          enddo
#else
          do ib=1,nT_h(iparm)
            potEmin=potEmin_all(ib,iparm)
            expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
            hfin(ind,ib)=hfin(ind,ib)+
     &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
             if (rmsrgymap) then
               indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
               indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
               hrmsrgy(indrgy,indrms,ib)=
     &           hrmsrgy(indrgy,indrms,ib)+expfac
             endif
          enddo
#endif
        enddo ! t
        do ib=1,nT_h(iparm)
          if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
     &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
          if (rmsrgymap) then
          call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
     &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
     &       WHAM_COMM,IERROR)
          endif
        enddo
        call MPI_Reduce(upindE_p,upindE,1,
     &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
        call MPI_Reduce(histE_p(0),histE(0),maxindE,
     &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)

        if (me1.eq.master) then

        if (histout) then

        write (iout,'(6x,$)')
        write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
     &   ib=1,nT_h(iparm))
        write (iout,*)

        write (iout,'(/a)') 'Final histograms'
        if (histfile) then
          if (nslice.eq.1) then
            if (separate_parset) then
              write(licz3,"(bz,i3.3)") myparm
              histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
            else
              histname=prefix(:ilen(prefix))//'.hist'
            endif
          else
            if (separate_parset) then
              write(licz3,"(bz,i3.3)") myparm
              histname=prefix(:ilen(prefix))//'_par'//licz3//
     &         '_slice_'//licz2//'.hist'
            else
              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
            endif
          endif
#if defined(AIX) || defined(PGI)
          open (ihist,file=histname,position='append')
#else
          open (ihist,file=histname,access='append')
#endif
        endif

        do t=0,tmax
          liczba=t
          sumH=0.0d0
          do ib=1,nT_h(iparm)
            sumH=sumH+hfin(t,ib)
          enddo
          if (sumH.gt.0.0d0) then
            do j=1,nQ
              jj = mod(liczba,nbin1)
              liczba=liczba/nbin1
              write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
              if (histfile) 
     &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
            enddo
            do ib=1,nT_h(iparm)
              write (iout,'(e20.10,$)') hfin(t,ib)
              if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
            enddo
            write (iout,'(i5)') iparm
            if (histfile) write (ihist,'(i5)') iparm
          endif
        enddo

        endif

        if (entfile) then
          if (nslice.eq.1) then
            if (separate_parset) then
              write(licz3,"(bz,i3.3)") myparm
              histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
            else
              histname=prefix(:ilen(prefix))//'.ent'
            endif
          else
            if (separate_parset) then
              write(licz3,"(bz,i3.3)") myparm
              histname=prefix(:ilen(prefix))//'par_'//licz3//
     &           '_slice_'//licz2//'.ent'
            else
              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
            endif
          endif
#if defined(AIX) || defined(PGI)
          open (ihist,file=histname,position='append')
#else
          open (ihist,file=histname,access='append')
#endif
          write (ihist,'(a)') "# Microcanonical entropy"
          do i=0,upindE
            write (ihist,'(f8.0,$)') dint(potEmin)+i
            if (histE(i).gt.0.0e0) then
              write (ihist,'(f15.5,$)') dlog(histE(i))
            else
              write (ihist,'(f15.5,$)') 0.0d0
            endif
          enddo
          write (ihist,*)
          close(ihist)
        endif
        write (iout,*) "Microcanonical entropy"
        do i=0,upindE
          write (iout,'(f8.0,$)') dint(potEmin)+i
          if (histE(i).gt.0.0e0) then
            write (iout,'(f15.5,$)') dlog(histE(i))
          else
            write (iout,'(f15.5,$)') 0.0d0
          endif
          write (iout,*)
        enddo
        if (rmsrgymap) then
          if (nslice.eq.1) then
            if (separate_parset) then
              write(licz3,"(bz,i3.3)") myparm
              histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
            else
              histname=prefix(:ilen(prefix))//'.rmsrgy'
            endif
          else
            if (separate_parset) then
              write(licz3,"(bz,i3.3)") myparm
              histname=prefix(:ilen(prefix))//'_par'//licz3//
     &         '_slice_'//licz2//'.rmsrgy'
            else
             histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
            endif
          endif
#if defined(AIX) || defined(PGI)
          open (ihist,file=histname,position='append')
#else
          open (ihist,file=histname,access='append')
#endif
          do i=0,nbin_rms
            do j=0,nbin_rgy
              write(ihist,'(2f8.2,$)') 
     &          rgymin+deltrgy*j,rmsmin+deltrms*i
              do ib=1,nT_h(iparm)
                if (hrmsrgy(j,i,ib).gt.0.0d0) then
                  write(ihist,'(e14.5,$)') 
     &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
     &              +potEmin
                else
                  write(ihist,'(e14.5,$)') 1.0d6
                endif
              enddo
              write (ihist,'(i2)') iparm
            enddo
          enddo
          close(ihist)
        endif
        endif
      enddo ! iparm
#ifdef MPI
      call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
     &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
      call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
     &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
      call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
     &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
      call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
     &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
      call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
     &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
      call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
     &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
     &   WHAM_COMM,IERROR)
      call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
     &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
     &   WHAM_COMM,IERROR)
      call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
     &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
     &   WHAM_COMM,IERROR)
      if (me.eq.master) then
#endif
      write (iout,'(/a)') 'Thermal characteristics of folding'
      if (nslice.eq.1) then
        nazwa=prefix
      else
        nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
      endif
      iln=ilen(nazwa)
      if (nparmset.eq.1 .and. .not.separate_parset) then
        nazwa=nazwa(:iln)//".thermal"
      else if (nparmset.eq.1 .and. separate_parset) then
        write(licz3,"(bz,i3.3)") myparm
        nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
      endif
      do iparm=1,nParmSet
      if (nparmset.gt.1) then
        write(licz3,"(bz,i3.3)") iparm
        nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
      endif
      open(34,file=nazwa)
      if (separate_parset) then
        write (iout,'(a,i3)') "Parameter set",myparm
      else
        write (iout,'(a,i3)') "Parameter set",iparm
      endif
      do i=0,NGridT
        betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
        if (betaT.ge.beta_h(1,iparm)) then
          potEmin=potEmin_all(1,iparm)
        else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
          potEmin=potEmin_all(nT_h(iparm),iparm)
        else
          do l=1,nT_h(iparm)-1
            if (betaT.le.beta_h(l,iparm) .and.
     &          betaT.gt.beta_h(l+1,iparm)) then
              potEmin=potEmin_all(l,iparm)
              exit
            endif
          enddo
        endif

c        write (iout,*) "i",i," potEmin",potEmin

        sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
        sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
     &    sumW(i,iparm)
        sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
     &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
        do j=1,nQ+2
          sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
          sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
     &     -sumQ(j,i,iparm)**2
          sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
     &     -sumQ(j,i,iparm)*sumE(i,iparm)
        enddo
        sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
     &   (startGridT+i*delta_T))+potEmin
        write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
     &   sumW(i,iparm),sumE(i,iparm)
        write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
        write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
     &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
        write (iout,*) 
        write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
     &   sumW(i,iparm),sumE(i,iparm)
        write(frm_write,'( "(",i3,"e15.5,$)" )' ) nQ+2
        write (34,frm_write) (sumQ(j,i,iparm),j=1,nQ+2)
        write(frm_write,'( "(",i3,"e15.5,$)" )' ) (nQ+2)*2+1
        write (34,frm_write) sumEsq(i,iparm)-sumEbis(i,iparm),
     &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
        write (34,*) 
      enddo
      close(34)
      enddo
      if (histout) then
      do t=0,tmax
        if (hfin_ent(t).gt.0.0d0) then
          liczba=t
          jj = mod(liczba,nbin1)
          write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
     &     hfin_ent(t)
          if (histfile) write (ihist,'(f6.3,e20.10," ent")') 
     &      dmin+(jj+0.5d0)*delta,
     &     hfin_ent(t)
        endif
      enddo
      if (histfile) close(ihist)
      endif

#ifdef ZSCORE
! Write data for zscore
      if (nslice.eq.1) then
        zscname=prefix(:ilen(prefix))//".zsc"
      else
        zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
      endif
#if defined(AIX) || defined(PGI)
      open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
#else
      open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
#endif
      write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
      do iparm=1,nParmSet
        write (izsc,'("NT=",i1)') nT_h(iparm)
      do ib=1,nT_h(iparm)
        write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') 
     &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
        jj = min0(nR(ib,iparm),7)
        write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
        write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
        write (izsc,'("&")')
        if (nR(ib,iparm).gt.7) then
          do ii=8,nR(ib,iparm),9
            jj = min0(nR(ib,iparm),ii+8)
            write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj) 
            write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
            write (izsc,'("&")')
          enddo
        endif
        write (izsc,'("FI=",$)')
        jj=min0(nR(ib,iparm),7)
        write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
        write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
        write (izsc,'("&")')
        if (nR(ib,iparm).gt.7) then
          do ii=8,nR(ib,iparm),9
            jj = min0(nR(ib,iparm),ii+8)
            write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj) 
            if (jj.eq.nR(ib,iparm)) then
              write (izsc,*) 
            else
              write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
              write (izsc,'(t80,"&")')
            endif
          enddo
        endif
        do i=1,nR(ib,iparm)
          write (izsc,'("KH=",$)') 
          write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
          write (izsc,'(" Q0=",$)')
          write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
          write (izsc,*)
        enddo
      enddo
      enddo
      close(izsc)
#endif
#ifdef MPI
      endif
#endif

      return

      end