subroutine make_ensembles(islice,*)
! construct the conformational ensembles at REMD temperatures
      implicit none
      include "DIMENSIONS"
      include "DIMENSIONS.ZSCOPT"
      include "DIMENSIONS.FREE"
#ifdef MPI
      include "mpif.h"
      include "COMMON.MPI"
      integer ierror,errcode,status(MPI_STATUS_SIZE) 
#endif
      include "COMMON.IOUNITS"
      include "COMMON.CONTROL"
      include "COMMON.HOMOLOGY"
      include "COMMON.FREE"
      include "COMMON.ENERGIES"
      include "COMMON.FFIELD"
      include "COMMON.INTERACT"
      include "COMMON.SBRIDGE"
      include "COMMON.CHAIN"
      include "COMMON.PROTFILES"
      include "COMMON.PROT"
      real*4 csingle(3,maxres2)
      double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
     &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
      double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
     &      escloc,eliptran,esaxs,
     &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
     &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt,
     &      ehomology_constr,edfadis,edfator,edfanei,edfabet
      integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
      double precision qfree,sumprob,eini,efree,rmsdev
      character*80 bxname
      character*2 licz1,licz2
      character*3 licz3,licz4
      character*5 ctemper
      integer ilen
      external ilen
      real*4 Fdimless(MaxStr),Fdimless_(MaxStr)
      double precision enepot(MaxStr)
      integer iperm(MaxStr)
      integer islice
 
#ifdef MPI
      if (me.eq.Master) then
#endif
      write (licz2,'(bz,i2.2)') islice
      if (nslice.eq.1) then
        if (.not.separate_parset) then
          bxname = prefix(:ilen(prefix))//".bx"
        else
          write (licz3,'(bz,i3.3)') myparm
          bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
        endif
      else
        if (.not.separate_parset) then
          bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
        else
          write (licz3,'(bz,i3.3)') myparm
          bxname = prefix(:ilen(prefix))//"par_"//licz3//
     &      "_slice_"//licz2//".bx"
        endif
      endif
      open (ientout,file=bxname,status="unknown",
     &  form="unformatted",access="direct",recl=lenrec1)
#ifdef MPI
      endif
#endif
      do iparm=1,nParmSet
        if (iparm.ne.iparmprint) exit
        call restore_parm(iparm)
        do ib=1,nT_h(iparm)
#ifdef DEBUG
          write (iout,*) "iparm",iparm," ib",ib
#endif
          temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
c          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
c          quotl=1.0d0
c          kfacl=1.0d0
c          do l=1,5
c            quotl1=quotl
c            quotl=quotl*quot
c            kfacl=kfacl*kfac
c            fT(l)=kfacl/(kfacl-1.0d0+quotl)
c          enddo
            if (rescale_mode.eq.1) then
              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
#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)=quot
#else
              ft(6)=1.0d0
#endif
              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
            else if (rescale_mode.eq.2) then
              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
#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))/3200.d0
#elif defined(FUNCT)
              ft(6)=quot
#else 
              ft(6)=1.0d0
#endif
              quotl=1.0d0
              do l=1,5
                quotl=quotl*quot
                fT(l)=1.12692801104297249644d0/
     &             dlog(dexp(quotl)+dexp(-quotl))
              enddo
c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
            else if (rescale_mode.eq.0) then
              do l=1,5
                fT(l)=0.0d0
              enddo
            else
              write (iout,*)
     &        "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
              call flush(iout)
              return1
            endif
#ifdef MPI
          do i=1,scount(me1)
#else
          do i=1,ntot(islice)
#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)
            eliptran=enetb(22,i,iparm)
            esaxs=enetb(26,i,iparm)
            ehomology_constr=enetb(27,i,iparm)
            edfadis=enetb(28,i,iparm)
            edfator=enetb(29,i,iparm)
            edfanei=enetb(30,i,iparm)
            edfabet=enetb(31,i,iparm)
            if (homol_nset.gt.1)
     &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
#ifdef SPLITELE
          if (shield_mode.gt.0) then
            etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
     &      +ft(1)*welec*ees
     &      +ft(1)*wvdwpp*evdw1
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+nss*ebr+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+wliptran*eliptran+wsaxs*esaxs
     &      +ehomology_constr
     &      +wdfa_dist*edfadis
     &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
          else
            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+nss*ebr+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+wliptran*eliptran+wsaxs*esaxs
     &      +ehomology_constr
     &      +wdfa_dist*edfadis
     &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
          endif
#else
          if (shield_mode.gt.0) then
            etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
     &      +ft(1)*welec*(ees+evdw1)
     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
     &      +wstrain*ehpb+nss*ebr+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+wliptran*eliptran+wsaxs*esaxs
     &      +ehomology_constr
     &      +wdfa_dist*edfadis
     &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
          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+nss*ebr+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+wliptran*eliptran+wsaxs*esaxs
     &      +ehomology_constr
     &      +wdfa_dist*edfadis
     &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
          endif
#endif
#ifdef MPI
            Fdimless_(i)=
     &        beta_h(ib,iparm)*etot-entfac(i)
            potE(i,iparm)=etot
#ifdef DEBUG
            write (iout,*) i,indstart(me)+i-1,ib,
     &       1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
     &       -entfac(i),Fdimless(i)
#endif
#else
            Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
            potE(i,iparm)=etot
#endif
          enddo   ! i
#ifdef MPI
          call MPI_Gatherv(Fdimless_(1),scount(me),
     &     MPI_REAL,Fdimless(1),
     &     scount(0),idispl(0),MPI_REAL,Master,
     &     WHAM_COMM, IERROR)
#ifdef DEBUG
          call MPI_Gatherv(potE(1,iparm),scount(me),
     &     MPI_DOUBLE_PRECISION,potE(1,iparm),
     &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
     &     WHAM_COMM, IERROR)
          call MPI_Gatherv(entfac(1),scount(me),
     &     MPI_DOUBLE_PRECISION,entfac(1),
     &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
     &     WHAM_COMM, IERROR)
#endif
          if (me.eq.Master) then
#ifdef DEBUG
          write (iout,*) "The FDIMLESS array before sorting"
          do i=1,ntot(islice)
            write (iout,*) i,fdimless(i)
          enddo
#endif
#endif
          do i=1,ntot(islice)
            iperm(i)=i
          enddo
          call mysort1(ntot(islice),Fdimless,iperm)
#ifdef DEBUG
          write (iout,*) "The FDIMLESS array after sorting"
          do i=1,ntot(islice)
            write (iout,*) i,iperm(i),fdimless(i)
          enddo
#endif
          qfree=0.0d0
          do i=1,ntot(islice)
            qfree=qfree+exp(-fdimless(i)+fdimless(1))
          enddo
c          write (iout,*) "qfree",qfree
          nlist=1
          sumprob=0.0
          do i=1,min0(ntot(islice),ensembles) 
            sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
#ifdef DEBUG
            write (iout,*) i,ib,beta_h(ib,iparm),
     &       1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
     &       potE(iperm(i),iparm),
     &       -entfac(iperm(i)),fdimless(i),sumprob
#endif
            if (sumprob.gt.0.99d0) goto 122
            nlist=nlist+1
          enddo  
  122     continue
#ifdef MPI
          endif
          call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM, 
     &       IERROR)
          call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
     &       IERROR)
          do i=1,nlist
            ii=iperm(i)
            iproc=0
            do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
              iproc=iproc+1
            enddo
            if (iproc.ge.nprocs) then
              write (iout,*) "Fatal error: processor out of range",iproc
              call flush(iout)
              if (bxfile) then
                close (ientout)
              else
                close (ientout,status="delete")
              endif
              return1
            endif
            ik=ii-indstart(iproc)+1
            if (iproc.ne.Master) then
              if (me.eq.iproc) then
#ifdef DEBUG
                write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
     &           " energy",potE(ik,iparm)
#endif
                call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
     &            Master,i,WHAM_COMM,IERROR)
              else if (me.eq.Master) then
                call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
     &            WHAM_COMM,STATUS,IERROR)
              endif
            else if (me.eq.Master) then
              enepot(i)=potE(ik,iparm)
            endif
          enddo
#else
          do i=1,nlist
            enepot(i)=potE(iperm(i),iparm)
          enddo
#endif
#ifdef MPI
          if (me.eq.Master) then
#endif
          write(licz3,'(bz,i3.3)') iparm
          write(licz2,'(bz,i2.2)') islice
          if (temper.lt.100.0d0) then
            write(ctemper,'(f3.0)') temper
          else if (temper.lt.1000.0) then
            write (ctemper,'(f4.0)') temper
          else
            write (ctemper,'(f5.0)') temper
          endif
          if (nparmset.eq.1) then
            if (separate_parset) then
              write(licz4,'(bz,i3.3)') myparm
              pdbname=prefix(:ilen(prefix))//"_par"//licz4
            else
              pdbname=prefix(:ilen(prefix))
            endif
          else
            pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
          endif
          if (nslice.eq.1) then
            pdbname=pdbname(:ilen(pdbname))//"_T_"//
     &        ctemper(:ilen(ctemper))//"pdb"
          else
            pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
     &        ctemper(:ilen(ctemper))//"pdb" 
          endif
          open(ipdb,file=pdbname)
          write (iout,*) "Before reading nlist",nlist
          do i=1,nlist
            read (ientout,rec=iperm(i))
     &        ((csingle(l,k),l=1,3),k=1,nres),
     &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
     &        nss,(ihpb(k),jhpb(k),k=1,nss),
     &        eini,efree,rmsdev,iscor
            do j=1,2*nres
              do k=1,3
                c(k,j)=csingle(k,j)
              enddo
            enddo
            eini=fdimless(i)
            call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
          enddo
#ifdef MPI
          endif
#endif
        enddo     ! ib
      enddo       ! iparm
      if (bxfile) then
        close(ientout)
      else
        close(ientout,status="delete")
      endif
      return
      end
!--------------------------------------------------
      subroutine mysort1(n, x, ipermut)
      implicit none
      integer i,j,imax,ipm,n
      real x(n)
      integer ipermut(n)
      real xtemp
      do i=1,n
        xtemp=x(i)
        imax=i
        do j=i+1,n
          if (x(j).lt.xtemp) then
            imax=j
            xtemp=x(j)
          endif
        enddo
        x(imax)=x(i)
        x(i)=xtemp
        ipm=ipermut(imax)
        ipermut(imax)=ipermut(i)
        ipermut(i)=ipm
      enddo
      return
      end