unres_package_Oct_2016 from emilial
[unres4.git] / source / cluster / cluster.f90
diff --git a/source/cluster/cluster.f90 b/source/cluster/cluster.f90
new file mode 100644 (file)
index 0000000..1b6767c
--- /dev/null
@@ -0,0 +1,664 @@
+      program cluster
+!
+! Program to cluster united-residue MCM results.
+!
+      use clust_data
+      use probability
+      use tracking
+      use hc_
+      use io_clust
+!#define CLUSTER
+      use io_units
+      use io_base, only: permut
+      use geometry_data, only: nres,theta,phi,alph,omeg,&
+                         c,cref
+      use energy_data, only: nnt,nct
+      use control_data, only: symetr,outpdb,outmol2,titel,&
+                          iopt,print_dist,MaxProcs
+      use control, only: tcpu,initialize
+
+      use wham_data, only: punch_dist
+      use io_wham, only: parmread
+      use work_part 
+!      include 'DIMENSIONS'
+!      include 'sizesclu.dat'
+#ifdef MPI
+      use mpi_data
+      implicit none
+      include "mpif.h"
+      integer :: IERROR,ERRCODE !STATUS(MPI_STATUS_SIZE)
+#else
+      implicit none
+!      include "COMMON.MPI"
+#endif
+!      include 'COMMON.TIME1'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.GEO'
+!      include 'COMMON.HEADER'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CLUSTER'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.FREE'
+
+      logical,dimension(:),allocatable :: printang !(max_cut)
+      integer,dimension(:),allocatable :: printpdb !(max_cut)
+      integer,dimension(:),allocatable :: printmol2 !(max_cut)
+      character(len=240) lineh
+      REAL(kind=4),dimension(:),allocatable :: CRIT,MEMBR !(maxconf)
+      REAL(kind=4),dimension(:),allocatable :: CRITVAL !(maxconf-1)
+      INTEGER,dimension(:),allocatable :: IA,IB !(maxconf)
+      INTEGER,dimension(:,:),allocatable :: ICLASS !(maxconf,maxconf-1)
+      INTEGER,dimension(:),allocatable :: HVALS !(maxconf-1)
+      INTEGER,dimension(:),allocatable :: IORDER,HEIGHT !(maxconf-1)
+      integer,dimension(:),allocatable :: nn !(maxconf)
+      integer :: ndis
+      real(kind=4),dimension(:),allocatable :: DISNN !(maxconf)
+      LOGICAL,dimension(:),allocatable :: FLAG !(maxconf)
+      integer :: i,j,k,l,m,n,len,lev,idum,ii,ind,jj,icut,ncon,&
+        it,ncon_work,ind1,kkk
+      real(kind=8) :: t1,t2,difconf
+      
+      real(kind=8),dimension(:),allocatable :: varia !(maxvar)
+      real(kind=8),dimension(:),allocatable :: list_conf_ !(maxvar)
+      real(kind=8) :: hrtime,mintime,sectime
+      logical :: eof
+      external :: difconf
+!el
+      real(kind=4),dimension(:),allocatable :: diss_ !(maxdist)
+      integer,dimension(:),allocatable :: scount_ !(maxdist)
+#ifdef MPI
+      call MPI_Init( IERROR )
+      call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
+      call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
+      Master = 0
+      if (ierror.gt.0) then
+        write(iout,*) "SEVERE ERROR - Can't initialize MPI."
+        call mpi_finalize(ierror)
+        stop
+      endif
+      if (nprocs.gt.MaxProcs+1) then
+        write (2,*) "Error - too many processors",&
+         nprocs,MaxProcs+1
+        write (2,*) "Increase MaxProcs and recompile"
+        call MPI_Finalize(IERROR)
+        stop
+      endif
+#endif
+!elwrite(iout,*) "before parmread"
+      allocate(printang(max_cut))
+      allocate(printpdb(max_cut))
+      allocate(printmol2(max_cut))
+      call initialize
+!elwrite(iout,*) "before parmread"
+      call openunits
+!elwrite(iout,*) "before parmread"
+      call parmread
+      call read_control
+!elwrite(iout,*) "after read control"
+      call molread
+!      if (refstr) call read_ref_structure(*30)
+      do i=1,nres
+        phi(i)=0.0D0
+        theta(i)=0.0D0
+        alph(i)=0.0D0
+        omeg(i)=0.0D0
+      enddo
+!      write (iout,*) "Before permut"
+!       write (iout,*) "symetr", symetr
+!      call flush(iout)
+      call permut(symetr)
+!      write (iout,*) "after permut"
+!      call flush(iout)
+      print *,'MAIN: nnt=',nnt,' nct=',nct
+
+      DO I=1,NCUT
+        PRINTANG(I)=.FALSE.
+        PRINTPDB(I)=0
+        printmol2(i)=0
+        IF (RCUTOFF(I).LT.0.0) THEN
+          RCUTOFF(I)=ABS(RCUTOFF(I))
+          PRINTANG(I)=.TRUE.
+          PRINTPDB(I)=outpdb
+          printmol2(i)=outmol2
+        ENDIF
+      ENDDO
+      write (iout,*) 'Number of cutoffs:',NCUT
+      write (iout,*) 'Cutoff values:'
+      DO ICUT=1,NCUT
+        WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),&
+          printpdb(icut),printmol2(icut)
+      ENDDO
+      DO I=1,NRES-3  
+        MULT(I)=1
+      ENDDO
+      allocate(list_conf(maxconf))
+      do i=1,maxconf
+        list_conf(i)=i
+      enddo
+      call read_coords(ncon,*20)
+
+      allocate(list_conf_(maxconf))
+      do i=1,maxconf
+        list_conf_(i)=list_conf(i)
+      enddo
+      deallocate(list_conf) 
+      allocate(list_conf(ncon))
+      do i=1,ncon
+        list_conf(i)=list_conf_(i)
+      enddo
+      deallocate(list_conf_) 
+
+!el      call alloc_clust_arrays(ncon)
+      write (iout,*) 'from read_coords: ncon',ncon
+      
+      write (iout,*) "nT",nT
+      do iT=1,nT
+      write (iout,*) "iT",iT
+#ifdef MPI
+      call work_partition(.true.,ncon)
+#endif
+!elwrite(iout,*)"after work partition, ncon_work=", ncon_work,ncon
+
+      call probabl(iT,ncon_work,ncon,*20)
+
+!elwrite(iout,*)"after probabl, ncon_work=", ncon_work,ncon
+
+      if (ncon_work.lt.2) then
+        write (iout,*) "Too few conformations; clustering skipped"
+        exit
+      endif
+#ifdef MPI
+      ndis=ncon_work*(ncon_work-1)/2
+      call work_partition(.true.,ndis)
+#endif
+!el      call alloc_clust_arrays(ncon_work)
+      allocate(ICC(ncon_work))
+      allocate(DISS(maxdist))
+
+      DO I=1,NCON_work
+        ICC(I)=I
+      ENDDO
+      WRITE (iout,'(A80)') TITEL
+      t1=tcpu()
+!
+! CALCULATE DISTANCES
+!
+      call daread_ccoords(1,ncon_work)
+      ind1=0
+      DO I=1,NCON_work-1
+        if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
+        do k=1,2*nres
+          do l=1,3
+            c(l,k)=allcart(l,k,i)
+          enddo 
+        enddo
+        kkk=1
+        do k=1,nres
+          do l=1,3
+            cref(l,k,kkk)=c(l,k)
+          enddo
+        enddo
+        DO J=I+1,NCON_work
+          IND=IOFFSET(NCON_work,I,J)
+#ifdef MPI
+          if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
+#endif
+          ind1=ind1+1
+          DISS(IND1)=DIFCONF(I,J)
+!          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
+#ifdef MPI
+          endif
+#endif
+        ENDDO
+      ENDDO
+      t2=tcpu()
+      WRITE (iout,'(/a,1pe14.5,a/)') &
+       'Time for distance calculation:',T2-T1,' sec.'
+      t1=tcpu()
+      PRINT '(a)','End of distance computation'
+!el---------------
+      allocate(diss_(maxdist))
+      allocate(scount_(0:nprocs))
+      
+      do i=1,maxdist
+        diss_(i)=diss(i)
+      enddo
+      do i=0,nprocs
+        scount_(i)=scount(i)
+      enddo
+!el-----------
+#ifdef MPI
+      call MPI_Gatherv(diss_(1),scount_(me),MPI_REAL,diss(1),&
+           scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
+      if (me.eq.master) then
+#endif
+      deallocate(diss_)
+      deallocate(scount_)
+      open(80,file='/tmp/distance',form='unformatted')
+      do i=1,ndis
+        write(80) diss(i)
+      enddo
+      if (punch_dist) then
+        do i=1,ncon_work-1
+          do j=i+1,ncon_work
+            IND=IOFFSET(NCON,I,J)
+            write (jrms,'(2i5,2f10.5)') i,j,diss(IND),&
+              energy(j)-energy(i)
+          enddo
+        enddo
+      endif
+!
+! Print out the RMS deviation matrix.
+!
+      if (print_dist) CALL DISTOUT(NCON_work)
+!
+!  call hierarchical clustering HC from F. Murtagh
+!
+      N=NCON_work
+      LEN = (N*(N-1))/2
+      write(iout,*) "-------------------------------------------"
+      write(iout,*) "HIERARCHICAL CLUSTERING using"
+      if (iopt.eq.1) then
+        write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
+      elseif (iopt.eq.2) then
+        write(iout,*) "SINGLE LINK METHOD"
+      elseif (iopt.eq.3) then
+        write(iout,*) "COMPLETE LINK METHOD"
+      elseif (iopt.eq.4) then
+        write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
+      elseif (iopt.eq.5) then
+        write(iout,*) "MCQUITTY'S METHOD"
+      elseif (iopt.eq.6) then
+        write(iout,*) "MEDIAN (GOWER'S) METHOD"
+      elseif (iopt.eq.7) then
+        write(iout,*) "CENTROID METHOD"
+      else
+        write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
+        write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
+        stop
+      endif
+      write(iout,*)
+      write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
+      write(iout,*) "February 1986"
+      write(iout,*) "References:"
+      write(iout,*) "1. Multidimensional clustering algorithms"
+      write(iout,*) "   Fionn Murtagh"
+      write(iout,*) "   Vienna : Physica-Verlag, 1985."
+      write(iout,*) "2. Multivariate data analysis"
+      write(iout,*) "   Fionn Murtagh and Andre Heck"
+      write(iout,*) "   Kluwer Academic Publishers, 1987"
+      write(iout,*) "-------------------------------------------"
+      write(iout,*)
+
+#ifdef DEBUG
+      write (iout,*) "The TOTFREE array"
+      do i=1,ncon_work
+        write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
+      enddo
+#endif
+      allocate(CRIT(N),MEMBR(N)) !(maxconf)
+      allocate(CRITVAL(N-1)) !(maxconf-1)
+      allocate(IA(N),IB(N))
+      allocate(ICLASS(N,N-1)) !(maxconf,maxconf-1)
+      allocate(HVALS(N-1)) !(maxconf-1)
+      allocate(IORDER(N-1),HEIGHT(N-1)) !(maxconf-1)
+      allocate(nn(N)) !(maxconf)
+      allocate(DISNN(N)) !(maxconf)
+      allocate(FLAG(N)) !(maxconf)
+      call flush(iout)
+      CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
+      LEV = N-1
+      write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
+      if (lev.lt.2) then
+        write (iout,*) "Too few conformations to cluster."
+        goto 192
+      endif
+      CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
+!      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
+
+      allocate(licz(maxgr))
+      allocate(iass(maxgr))
+      allocate(nconf(maxgr,maxingr))
+      allocate(totfree_gr(maxgr))
+
+      do i=1,maxgr
+        licz(i)=0
+      enddo
+      icut=1
+      i=1
+      NGR=i+1
+      do j=1,n
+        licz(iclass(j,i))=licz(iclass(j,i))+1
+        nconf(iclass(j,i),licz(iclass(j,i)))=j
+!        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
+!     &    nconf(iclass(j,i),licz(iclass(j,i)))
+      enddo        
+      do i=1,lev-1
+
+         idum=lev-i
+         DO L=1,LEV
+            IF (HEIGHT(L).EQ.IDUM) GOTO 190
+         ENDDO
+ 190     IDUM=L
+         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),&
+          " icut",icut," cutoff",rcutoff(icut)
+         IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
+          WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
+          write (iout,'(a,f8.2)') 'Maximum distance found:',&
+                    CRITVAL(IDUM)
+          CALL SRTCLUST(ICUT,ncon_work,iT)
+          CALL TRACK(ICUT)
+          CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
+          icut=icut+1
+          if (icut.gt.ncut) goto 191
+         ENDIF
+         NGR=i+1
+         do l=1,maxgr
+          licz(l)=0
+         enddo
+         do j=1,n
+         enddo
+         do j=1,n
+          licz(iclass(j,i))=licz(iclass(j,i))+1
+          nconf(iclass(j,i),licz(iclass(j,i)))=j
+!d        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),&
+!d         nconf(iclass(j,i),licz(iclass(j,i)))
+!d          print *,j,iclass(j,i),
+!d     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
+         enddo
+      enddo
+ 191  continue
+!
+      if (plot_tree) then
+        CALL WRITRACK
+        CALL PLOTREE
+      endif
+!
+      t2=tcpu()
+      WRITE (iout,'(/a,1pe14.5,a/)') &
+       'Total time for clustering:',T2-T1,' sec.'
+#ifdef MPI
+      endif
+#endif
+ 192  continue
+      enddo
+!
+      close(icbase,status="delete")
+#ifdef MPI
+!el      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
+#endif
+      stop '********** Program terminated normally.'
+   20 write (iout,*) "Error reading coordinates"
+#ifdef MPI
+!el      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
+#endif
+      stop
+   30 write (iout,*) "Error reading reference structure"
+#ifdef MPI
+!el      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
+#endif
+      stop
+      end program cluster
+!---------------------------------------------------------------------------
+!
+!---------------------------------------------------------------------------
+      real(kind=8) function difconf(icon,jcon)
+
+      use clust_data
+
+      use io_units, only: iout
+      use io_base, only: permut
+      use geometry_data, only: nres,c,cref,tabperm
+      use energy_data, only: nct,nnt
+      use control_data, only: symetr,lside,nend,nstart
+      use regularize_, only: fitsq
+      implicit none
+!      include 'DIMENSIONS'
+!      include 'sizesclu.dat'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.CLUSTER'
+!      include 'COMMON.CHAIN' 
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.VAR'
+!      include 'COMMON.IOUNITS'
+      logical :: non_conv
+      real(kind=8) :: przes(3),obrot(3,3)
+      real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
+      integer :: i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
+      integer :: iaperm,ibezperm,run
+      real(kind=8) :: rms,rmsmina
+!      write (iout,*) "tu dochodze"
+      rmsmina=10d10
+      nperm=1
+      do i=1,symetr
+      nperm=i*nperm
+      enddo
+!      write (iout,*) "nperm",nperm
+      call permut(symetr)
+!      write (iout,*) "tabperm", tabperm(1,1)
+      do kkk=1,nperm
+      if (lside) then
+        ii=0
+        chalen=int((nend-nstart+2)/symetr)
+        do run=1,symetr
+         do i=nstart,(nstart+chalen-1)
+          zzz=tabperm(kkk,run)
+!          write (iout,*) "tutaj",zzz
+          ii=ii+1
+          iaperm=(zzz-1)*chalen+i
+          ibezperm=(run-1)*chalen+i
+          do j=1,3
+            xx(j,ii)=allcart(j,iaperm,jcon)
+            yy(j,ii)=cref(j,ibezperm,kkk)
+          enddo
+         enddo
+        enddo
+        do run=1,symetr
+         do i=nstart,(nstart+chalen-1)
+          zzz=tabperm(kkk,run)
+          ii=ii+1
+          iaperm=(zzz-1)*chalen+i
+          ibezperm=(run-1)*chalen+i
+!          if (itype(i).ne.10) then
+            ii=ii+1
+            do j=1,3 
+              xx(j,ii)=allcart(j,iaperm+nres,jcon)
+              yy(j,ii)=cref(j,ibezperm+nres,kkk)
+            enddo
+           enddo
+!          endif
+        enddo
+        call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
+      else
+        chalen=int((nct-nnt+2)/symetr)
+        do run=1,symetr
+         do i=nnt,(nnt+chalen-1)
+          zzz=tabperm(kkk,run)
+!           write (iout,*) "tu szukaj", zzz,run,kkk
+          iaperm=(zzz-1)*chalen+i
+          ibezperm=(run-1)*chalen+i
+!        do i=nnt,nct
+          do j=1,3
+            c(j,i)=allcart(j,iaperm,jcon)
+          enddo
+         enddo
+        enddo
+        call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,&
+             przes,&
+             obrot,non_conv)
+      endif
+      if (rms.lt.0.0) then
+        print *,'error, rms^2 = ',rms,icon,jcon
+        stop
+      endif
+      if (non_conv) print *,non_conv,icon,jcon
+      if (rmsmina.gt.rms) rmsmina=rms
+      enddo
+      difconf=dsqrt(rmsmina)
+      return
+      end function difconf
+!------------------------------------------------------------------------------
+      subroutine distout(ncon)
+
+      use clust_data
+      use hc_, only:ioffset
+      use io_units, only: iout
+      implicit none
+!      include 'DIMENSIONS'
+!      include 'sizesclu.dat'
+      integer :: ncon
+      integer,parameter :: ncol=10
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CLUSTER'
+      integer :: i,j,k,jlim,jlim1,nlim,ind
+      real(kind=4) :: b(ncol)
+
+      write (iout,'(a)') 'The distance matrix'
+      do 1 i=1,ncon,ncol
+      nlim=min0(i+ncol-1,ncon)
+      write (iout,1000) (k,k=i,nlim)
+      write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
+ 1000 format (/8x,10(i4,3x))
+ 1020 format (/1x,80(1h-)/)
+      do 2 j=i,ncon
+      jlim=min0(j,nlim)
+      if (jlim.eq.j) then
+        b(jlim-i+1)=0.0d0
+        jlim1=jlim-1
+      else
+        jlim1=jlim
+      endif
+      do 3 k=i,jlim1
+       if (j.lt.k) then 
+          IND=IOFFSET(NCON,j,k)
+       else
+          IND=IOFFSET(NCON,k,j)
+       endif
+    3  b(k-i+1)=diss(IND)
+      write (iout,1010) j,(b(k),k=1,jlim-i+1)
+    2 continue
+    1 continue
+ 1010 format (i5,3x,10(f6.2,1x))
+      return
+      end subroutine distout
+!------------------------------------------------------------------------------
+! srtclust.f
+!------------------------------------------------------------------------------
+      SUBROUTINE SRTCLUST(ICUT,NCON,IB)
+
+      use clust_data
+      use io_units, only: iout
+!      implicit real*8 (a-h,o-z)
+!     include 'DIMENSIONS'
+!     include 'sizesclu.dat'
+!     include 'COMMON.CLUSTER'
+!     include 'COMMON.FREE'
+!     include 'COMMON.IOUNITS'
+      implicit none
+      real(kind=8),dimension(:),allocatable :: prob !(maxgr)
+      real(kind=8) :: emin,ene,en1,sumprob
+      integer :: igr,i,ii,li1,li2,ligr,ico,jco,ind1,ind2
+      integer :: jgr,li,nco,ib,ncon,icut
+!
+! Compute free energies of clusters
+!
+      allocate(prob(maxgr))
+
+      do igr=1,ngr
+      emin=totfree(nconf(igr,1))
+      totfree_gr(igr)=1.0d0
+      do i=2,licz(igr)
+        ii=nconf(igr,i)
+        totfree_gr(igr)=totfree_gr(igr)+dexp(-totfree(ii)+emin)
+      enddo
+!      write (iout,*) "igr",igr," totfree",emin,
+!     &    " totfree_gr",totfree_gr(igr)
+      totfree_gr(igr)=emin-dlog(totfree_gr(igr))
+!      write (iout,*) igr," efree",totfree_gr(igr)/beta_h(ib)
+      enddo
+!
+!  SORT CONFORMATIONS IN GROUPS ACC. TO ENERGY
+!
+      DO 16 IGR=1,NGR
+      LIGR=LICZ(IGR)
+      DO 17 ICO=1,LIGR-1
+      IND1=NCONF(IGR,ICO)
+      ENE=totfree(IND1)
+      DO 18 JCO=ICO+1,LIGR
+      IND2=NCONF(IGR,JCO)
+      EN1=totfree(IND2)
+      IF (EN1.LT.ENE) THEN
+        NCONF(IGR,ICO)=IND2
+        NCONF(IGR,JCO)=IND1
+        IND1=IND2
+        ENE=EN1
+      ENDIF
+   18 CONTINUE
+   17 CONTINUE
+   16 CONTINUE
+!
+!  SORT GROUPS
+!
+      DO 71 IGR=1,NGR
+      ENE=totfree_gr(IGR)
+      DO 72 JGR=IGR+1,NGR
+      EN1=totfree_gr(JGR)
+      IF (EN1.LT.ENE) THEN
+        LI1=LICZ(IGR)
+        LI2=LICZ(JGR)
+        LI=MAX0(LI1,LI2)
+        DO 73 I=1,LI   
+        NCO=NCONF(IGR,I)
+        NCONF(IGR,I)=NCONF(JGR,I)
+        NCONF(JGR,I)=NCO
+   73   CONTINUE
+        totfree_gr(igr)=en1
+        totfree_gr(jgr)=ene
+        ENE=EN1
+        LICZ(IGR)=LI2
+        LICZ(JGR)=LI1
+      ENDIF
+   72 CONTINUE
+   71 CONTINUE
+      write (iout,'("Free energies and probabilities of clusters at",f6.1," K")')&
+           1.0d0/(1.987d-3*beta_h(ib)) !'
+      prob(1)=1.0d0
+      sumprob=1.0d0
+      do i=2,ngr
+        prob(i)=dexp(-(totfree_gr(i)-totfree_gr(1)))
+        sumprob=sumprob+prob(i)
+      enddo
+      do i=1,ngr
+        prob(i)=prob(i)/sumprob
+      enddo
+      sumprob=0.0d0
+      write (iout,'("clust   efree    prob sumprob")')
+      do i=1,ngr
+        sumprob=sumprob+prob(i)
+        write (iout,'(i5,f8.1,2f8.5)') i,totfree_gr(i)/beta_h(ib),&
+          prob(i),sumprob
+      enddo
+      DO 81 IGR=1,NGR
+      LI=LICZ(IGR)
+      DO 82 I=1,LI 
+   82 IASS(NCONF(IGR,I))=IGR
+   81 CONTINUE
+      if (lgrp) then
+        do i=1,ncon
+          iass_tot(i,icut)=iass(i)
+!          write (iout,*) icut,i,iass(i),iass_tot(i,icut)
+        enddo
+      endif
+      RETURN
+      END SUBROUTINE SRTCLUST
+!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------