source/cluster/cluster.f90

   1       program cluster
   2 !
   3 ! Program to cluster united-residue MCM results.
   4 !
   5       use clust_data
   6       use probability
   7       use tracking
   8       use hc_
   9       use io_clust
  10 !#define CLUSTER
  11       use io_units
  12       use io_base, only: permut
  13       use geometry_data, only: nres,theta,phi,alph,omeg,&
  14                          c,cref
  15       use energy_data, only: nnt,nct
  16       use control_data, only: symetr,outpdb,outmol2,titel,&
  17                           iopt,print_dist !,MaxProcs
  18       use control, only: tcpu,initialize
  19
  20       use wham_data, only: punch_dist
  21       use io_wham, only: parmread
  22       use work_part
  23 !      include 'DIMENSIONS'
  24 !      include 'sizesclu.dat'
  25 #ifdef MPI
  26       use mpi_data
  27       implicit none
  28       include "mpif.h"
  29       integer :: IERROR,ERRCODE !STATUS(MPI_STATUS_SIZE)
  30 #else
  31       implicit none
  32 !      include "COMMON.MPI"
  33 #endif
  34 !      include 'COMMON.TIME1'
  35 !      include 'COMMON.INTERACT'
  36 !      include 'COMMON.NAMES'
  37 !      include 'COMMON.GEO'
  38 !      include 'COMMON.HEADER'
  39 !      include 'COMMON.CONTROL'
  40 !      include 'COMMON.CHAIN'
  41 !      include 'COMMON.VAR'
  42 !      include 'COMMON.CLUSTER'
  43 !      include 'COMMON.IOUNITS'
  44 !      include 'COMMON.FREE'
  45
  46       logical,dimension(:),allocatable :: printang !(max_cut)
  47       integer,dimension(:),allocatable :: printpdb !(max_cut)
  48       integer,dimension(:),allocatable :: printmol2 !(max_cut)
  49       character(len=240) lineh
  50       REAL(kind=4),dimension(:),allocatable :: CRIT,MEMBR !(maxconf)
  51       REAL(kind=4),dimension(:),allocatable :: CRITVAL !(maxconf-1)
  52       INTEGER,dimension(:),allocatable :: IA,IB !(maxconf)
  53       INTEGER,dimension(:,:),allocatable :: ICLASS !(maxconf,maxconf-1)
  54       INTEGER,dimension(:),allocatable :: HVALS !(maxconf-1)
  55       INTEGER,dimension(:),allocatable :: IORDER,HEIGHT !(maxconf-1)
  56       integer,dimension(:),allocatable :: nn !(maxconf)
  57       integer :: ndis
  58       real(kind=4),dimension(:),allocatable :: DISNN !(maxconf)
  59       LOGICAL,dimension(:),allocatable :: FLAG !(maxconf)
  60       integer :: i,j,k,l,m,n,len,lev,idum,ii,ind,jj,icut,ncon,&
  61         it,ncon_work,ind1,kkk
  62       real(kind=8) :: t1,t2,difconf
  63
  64       real(kind=8),dimension(:),allocatable :: varia !(maxvar)
  65       real(kind=8),dimension(:),allocatable :: list_conf_ !(maxvar)
  66       real(kind=8) :: hrtime,mintime,sectime
  67       logical :: eof
  68       external :: difconf
  69 !el
  70       real(kind=4),dimension(:),allocatable :: diss_ !(maxdist)
  71       integer,dimension(:),allocatable :: scount_ !(maxdist)
  72 #ifdef MPI
  73       call MPI_Init( IERROR )
  74       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
  75       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
  76       Master = 0
  77       if (ierror.gt.0) then
  78         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
  79         call mpi_finalize(ierror)
  80         stop
  81       endif
  82 !EL /04/2017  No longer needed, dynamic allocation of arrays
  83 !      if (nprocs.gt.MaxProcs+1) then
  84 !        write (2,*) "Error - too many processors",&
  85 !         nprocs,MaxProcs+1
  86 !        write (2,*) "Increase MaxProcs and recompile"
  87 !        call MPI_Finalize(IERROR)
  88 !        stop
  89 !      endif
  90 #endif
  91 !elwrite(iout,*) "before parmread"
  92       allocate(printang(max_cut))
  93       allocate(printpdb(max_cut))
  94       allocate(printmol2(max_cut))
  95       call initialize
  96 !elwrite(iout,*) "before parmread"
  97       call openunits
  98 !elwrite(iout,*) "before parmread"
  99       call parmread
 100       call read_control
 101 !elwrite(iout,*) "after read control"
 102       call molread
 103 !      if (refstr) call read_ref_structure(*30)
 104       do i=1,nres
 105         phi(i)=0.0D0
 106         theta(i)=0.0D0
 107         alph(i)=0.0D0
 108         omeg(i)=0.0D0
 109       enddo
 110 !      write (iout,*) "Before permut"
 111 !       write (iout,*) "symetr", symetr
 112 !      call flush(iout)
 113       call permut(symetr)
 114 !      write (iout,*) "after permut"
 115 !      call flush(iout)
 116       print *,'MAIN: nnt=',nnt,' nct=',nct
 117
 118       DO I=1,NCUT
 119         PRINTANG(I)=.FALSE.
 120         PRINTPDB(I)=0
 121         printmol2(i)=0
 122         IF (RCUTOFF(I).LT.0.0) THEN
 123           RCUTOFF(I)=ABS(RCUTOFF(I))
 124           PRINTANG(I)=.TRUE.
 125           PRINTPDB(I)=outpdb
 126           printmol2(i)=outmol2
 127         ENDIF
 128       ENDDO
 129       write (iout,*) 'Number of cutoffs:',NCUT
 130       write (iout,*) 'Cutoff values:'
 131       DO ICUT=1,NCUT
 132         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),&
 133           printpdb(icut),printmol2(icut)
 134       ENDDO
 135       DO I=1,NRES-3
 136         MULT(I)=1
 137       ENDDO
 138       allocate(list_conf(maxconf))
 139       do i=1,maxconf
 140         list_conf(i)=i
 141       enddo
 142       call read_coords(ncon,*20)
 143
 144       allocate(list_conf_(maxconf))
 145       do i=1,maxconf
 146         list_conf_(i)=list_conf(i)
 147       enddo
 148       deallocate(list_conf)
 149       allocate(list_conf(ncon))
 150       do i=1,ncon
 151         list_conf(i)=list_conf_(i)
 152       enddo
 153       deallocate(list_conf_)
 154
 155 !el      call alloc_clust_arrays(ncon)
 156
 157       write (iout,*) 'from read_coords: ncon',ncon
 158
 159       write (iout,*) "nT",nT
 160       do iT=1,nT
 161       write (iout,*) "iT",iT
 162 #ifdef MPI
 163       call work_partition(.true.,ncon)
 164 #endif
 165 !elwrite(iout,*)"after work partition, ncon_work=", ncon_work,ncon
 166
 167       call probabl(iT,ncon_work,ncon,*20)
 168
 169 !elwrite(iout,*)"after probabl, ncon_work=", ncon_work,ncon
 170
 171       if (ncon_work.lt.2) then
 172         write (iout,*) "Too few conformations; clustering skipped"
 173         exit
 174       endif
 175 #ifdef MPI
 176       ndis=ncon_work*(ncon_work-1)/2
 177       call work_partition(.true.,ndis)
 178 #endif
 179 !el      call alloc_clust_arrays(ncon_work)
 180       allocate(ICC(ncon_work))
 181       allocate(DISS(maxdist))
 182
 183       DO I=1,NCON_work
 184         ICC(I)=I
 185       ENDDO
 186       WRITE (iout,'(A80)') TITEL
 187       t1=tcpu()
 188 !
 189 ! CALCULATE DISTANCES
 190 !
 191       call daread_ccoords(1,ncon_work)
 192       ind1=0
 193       DO I=1,NCON_work-1
 194         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 195         do k=1,2*nres
 196           do l=1,3
 197             c(l,k)=allcart(l,k,i)
 198           enddo
 199         enddo
 200         kkk=1
 201         do k=1,nres
 202           do l=1,3
 203             cref(l,k,kkk)=c(l,k)
 204           enddo
 205         enddo
 206         DO J=I+1,NCON_work
 207           IND=IOFFSET(NCON_work,I,J)
 208 #ifdef MPI
 209           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
 210 #endif
 211           ind1=ind1+1
 212           DISS(IND1)=DIFCONF(I,J)
 213 !          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 214 #ifdef MPI
 215           endif
 216 #endif
 217         ENDDO
 218       ENDDO
 219       t2=tcpu()
 220       WRITE (iout,'(/a,1pe14.5,a/)') &
 221        'Time for distance calculation:',T2-T1,' sec.'
 222       t1=tcpu()
 223       PRINT '(a)','End of distance computation'
 224 !el---------------
 225       allocate(diss_(maxdist))
 226       allocate(scount_(0:nprocs))
 227
 228       do i=1,maxdist
 229         diss_(i)=diss(i)
 230       enddo
 231       do i=0,nprocs
 232         scount_(i)=scount(i)
 233       enddo
 234 !el-----------
 235 #ifdef MPI
 236       call MPI_Gatherv(diss_(1),scount_(me),MPI_REAL,diss(1),&
 237            scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 238       if (me.eq.master) then
 239 #endif
 240       deallocate(diss_)
 241       deallocate(scount_)
 242
 243       open(80,file='/tmp/distance',form='unformatted')
 244       do i=1,ndis
 245         write(80) diss(i)
 246       enddo
 247       if (punch_dist) then
 248         do i=1,ncon_work-1
 249           do j=i+1,ncon_work
 250             IND=IOFFSET(NCON,I,J)
 251             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),&
 252               energy(j)-energy(i)
 253           enddo
 254         enddo
 255       endif
 256 !
 257 ! Print out the RMS deviation matrix.
 258 !
 259       if (print_dist) CALL DISTOUT(NCON_work)
 260 !
 261 !  call hierarchical clustering HC from F. Murtagh
 262 !
 263       N=NCON_work
 264       LEN = (N*(N-1))/2
 265       write(iout,*) "-------------------------------------------"
 266       write(iout,*) "HIERARCHICAL CLUSTERING using"
 267       if (iopt.eq.1) then
 268         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 269       elseif (iopt.eq.2) then
 270         write(iout,*) "SINGLE LINK METHOD"
 271       elseif (iopt.eq.3) then
 272         write(iout,*) "COMPLETE LINK METHOD"
 273       elseif (iopt.eq.4) then
 274         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 275       elseif (iopt.eq.5) then
 276         write(iout,*) "MCQUITTY'S METHOD"
 277       elseif (iopt.eq.6) then
 278         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 279       elseif (iopt.eq.7) then
 280         write(iout,*) "CENTROID METHOD"
 281       else
 282         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 283         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 284         stop
 285       endif
 286       write(iout,*)
 287       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 288       write(iout,*) "February 1986"
 289       write(iout,*) "References:"
 290       write(iout,*) "1. Multidimensional clustering algorithms"
 291       write(iout,*) "   Fionn Murtagh"
 292       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 293       write(iout,*) "2. Multivariate data analysis"
 294       write(iout,*) "   Fionn Murtagh and Andre Heck"
 295       write(iout,*) "   Kluwer Academic Publishers, 1987"
 296       write(iout,*) "-------------------------------------------"
 297       write(iout,*)
 298
 299 #ifdef DEBUG
 300       write (iout,*) "The TOTFREE array"
 301       do i=1,ncon_work
 302         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
 303       enddo
 304 #endif
 305       allocate(CRIT(N),MEMBR(N)) !(maxconf)
 306       allocate(CRITVAL(N-1)) !(maxconf-1)
 307       allocate(IA(N),IB(N))
 308       allocate(ICLASS(N,N-1)) !(maxconf,maxconf-1)
 309       allocate(HVALS(N-1)) !(maxconf-1)
 310       allocate(IORDER(N-1),HEIGHT(N-1)) !(maxconf-1)
 311       allocate(nn(N)) !(maxconf)
 312       allocate(DISNN(N)) !(maxconf)
 313       allocate(FLAG(N)) !(maxconf)
 314
 315       call flush(iout)
 316       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 317       LEV = N-1
 318       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 319       if (lev.lt.2) then
 320         write (iout,*) "Too few conformations to cluster."
 321         goto 192
 322       endif
 323       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 324 !      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 325
 326       allocate(licz(maxgr))
 327       allocate(iass(maxgr))
 328       allocate(nconf(maxgr,maxingr))
 329       allocate(totfree_gr(maxgr))
 330
 331       do i=1,maxgr
 332         licz(i)=0
 333       enddo
 334       icut=1
 335       i=1
 336       NGR=i+1
 337       do j=1,n
 338         licz(iclass(j,i))=licz(iclass(j,i))+1
 339         nconf(iclass(j,i),licz(iclass(j,i)))=j
 340 !        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 341 !     &    nconf(iclass(j,i),licz(iclass(j,i)))
 342       enddo
 343       do i=1,lev-1
 344
 345          idum=lev-i
 346          DO L=1,LEV
 347             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 348          ENDDO
 349  190     IDUM=L
 350          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),&
 351           " icut",icut," cutoff",rcutoff(icut)
 352          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 353           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 354           write (iout,'(a,f8.2)') 'Maximum distance found:',&
 355                     CRITVAL(IDUM)
 356           CALL SRTCLUST(ICUT,ncon_work,iT)
 357           CALL TRACK(ICUT)
 358           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 359           icut=icut+1
 360           if (icut.gt.ncut) goto 191
 361          ENDIF
 362          NGR=i+1
 363          do l=1,maxgr
 364           licz(l)=0
 365          enddo
 366          do j=1,n
 367          enddo
 368          do j=1,n
 369           licz(iclass(j,i))=licz(iclass(j,i))+1
 370           nconf(iclass(j,i),licz(iclass(j,i)))=j
 371 !d        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),&
 372 !d         nconf(iclass(j,i),licz(iclass(j,i)))
 373 !d          print *,j,iclass(j,i),
 374 !d     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 375          enddo
 376       enddo
 377  191  continue
 378 !
 379       if (plot_tree) then
 380         CALL WRITRACK
 381         CALL PLOTREE
 382       endif
 383 !
 384       t2=tcpu()
 385       WRITE (iout,'(/a,1pe14.5,a/)') &
 386        'Total time for clustering:',T2-T1,' sec.'
 387 #ifdef MPI
 388       endif
 389 #endif
 390  192  continue
 391       enddo
 392 !
 393       close(icbase,status="delete")
 394 #ifdef MPI
 395 !el      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 396       call MPI_Finalize(IERROR)
 397 #endif
 398       stop '********** Program terminated normally.'
 399    20 write (iout,*) "Error reading coordinates"
 400 #ifdef MPI
 401 !el      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 402       call MPI_Finalize(IERROR)
 403 #endif
 404       stop
 405    30 write (iout,*) "Error reading reference structure"
 406 #ifdef MPI
 407 !el      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 408       call MPI_Finalize(IERROR)
 409 #endif
 410       stop
 411       end program cluster
 412 !---------------------------------------------------------------------------
 413 !
 414 !---------------------------------------------------------------------------
 415       real(kind=8) function difconf(icon,jcon)
 416
 417       use clust_data
 418
 419       use io_units, only: iout
 420       use io_base, only: permut
 421       use geometry_data, only: nres,c,cref,tabperm
 422       use energy_data, only: nct,nnt
 423       use control_data, only: symetr,lside,nend,nstart
 424       use regularize_, only: fitsq
 425       implicit none
 426 !      include 'DIMENSIONS'
 427 !      include 'sizesclu.dat'
 428 !      include 'COMMON.CONTROL'
 429 !      include 'COMMON.CLUSTER'
 430 !      include 'COMMON.CHAIN'
 431 !      include 'COMMON.INTERACT'
 432 !      include 'COMMON.VAR'
 433 !      include 'COMMON.IOUNITS'
 434       logical :: non_conv
 435       real(kind=8) :: przes(3),obrot(3,3)
 436       real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
 437       integer :: i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
 438       integer :: iaperm,ibezperm,run
 439       real(kind=8) :: rms,rmsmina
 440 !      write (iout,*) "tu dochodze"
 441       rmsmina=10d10
 442       nperm=1
 443       do i=1,symetr
 444       nperm=i*nperm
 445       enddo
 446 !      write (iout,*) "nperm",nperm
 447       call permut(symetr)
 448 !      write (iout,*) "tabperm", tabperm(1,1)
 449       do kkk=1,nperm
 450       if (lside) then
 451         ii=0
 452         chalen=int((nend-nstart+2)/symetr)
 453         do run=1,symetr
 454          do i=nstart,(nstart+chalen-1)
 455           zzz=tabperm(kkk,run)
 456 !          write (iout,*) "tutaj",zzz
 457           ii=ii+1
 458           iaperm=(zzz-1)*chalen+i
 459           ibezperm=(run-1)*chalen+i
 460           do j=1,3
 461             xx(j,ii)=allcart(j,iaperm,jcon)
 462             yy(j,ii)=cref(j,ibezperm,kkk)
 463           enddo
 464          enddo
 465         enddo
 466         do run=1,symetr
 467          do i=nstart,(nstart+chalen-1)
 468           zzz=tabperm(kkk,run)
 469           ii=ii+1
 470           iaperm=(zzz-1)*chalen+i
 471           ibezperm=(run-1)*chalen+i
 472 !          if (itype(i).ne.10) then
 473             ii=ii+1
 474             do j=1,3
 475               xx(j,ii)=allcart(j,iaperm+nres,jcon)
 476               yy(j,ii)=cref(j,ibezperm+nres,kkk)
 477             enddo
 478            enddo
 479 !          endif
 480         enddo
 481         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 482       else
 483         chalen=int((nct-nnt+2)/symetr)
 484         do run=1,symetr
 485          do i=nnt,(nnt+chalen-1)
 486           zzz=tabperm(kkk,run)
 487 !           write (iout,*) "tu szukaj", zzz,run,kkk
 488           iaperm=(zzz-1)*chalen+i
 489           ibezperm=(run-1)*chalen+i
 490 !        do i=nnt,nct
 491           do j=1,3
 492             c(j,i)=allcart(j,iaperm,jcon)
 493           enddo
 494          enddo
 495         enddo
 496         call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,&
 497              przes,&
 498              obrot,non_conv)
 499       endif
 500       if (rms.lt.0.0) then
 501         print *,'error, rms^2 = ',rms,icon,jcon
 502         stop
 503       endif
 504       if (non_conv) print *,non_conv,icon,jcon
 505       if (rmsmina.gt.rms) rmsmina=rms
 506       enddo
 507       difconf=dsqrt(rmsmina)
 508       return
 509       end function difconf
 510 !------------------------------------------------------------------------------
 511       subroutine distout(ncon)
 512
 513       use clust_data
 514       use hc_, only:ioffset
 515       use io_units, only: iout
 516       implicit none
 517 !      include 'DIMENSIONS'
 518 !      include 'sizesclu.dat'
 519       integer :: ncon
 520       integer,parameter :: ncol=10
 521 !      include 'COMMON.IOUNITS'
 522 !      include 'COMMON.CLUSTER'
 523       integer :: i,j,k,jlim,jlim1,nlim,ind
 524       real(kind=4) :: b(ncol)
 525
 526       write (iout,'(a)') 'The distance matrix'
 527       do 1 i=1,ncon,ncol
 528       nlim=min0(i+ncol-1,ncon)
 529       write (iout,1000) (k,k=i,nlim)
 530       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 531  1000 format (/8x,10(i4,3x))
 532  1020 format (/1x,80(1h-)/)
 533       do 2 j=i,ncon
 534       jlim=min0(j,nlim)
 535       if (jlim.eq.j) then
 536         b(jlim-i+1)=0.0d0
 537         jlim1=jlim-1
 538       else
 539         jlim1=jlim
 540       endif
 541       do 3 k=i,jlim1
 542        if (j.lt.k) then
 543           IND=IOFFSET(NCON,j,k)
 544        else
 545           IND=IOFFSET(NCON,k,j)
 546        endif
 547     3  b(k-i+1)=diss(IND)
 548       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 549     2 continue
 550     1 continue
 551  1010 format (i5,3x,10(f6.2,1x))
 552       return
 553       end subroutine distout
 554 !------------------------------------------------------------------------------
 555 ! srtclust.f
 556 !------------------------------------------------------------------------------
 557       SUBROUTINE SRTCLUST(ICUT,NCON,IB)
 558
 559       use clust_data
 560       use io_units, only: iout
 561 !      implicit real*8 (a-h,o-z)
 562 !     include 'DIMENSIONS'
 563 !     include 'sizesclu.dat'
 564 !     include 'COMMON.CLUSTER'
 565 !     include 'COMMON.FREE'
 566 !     include 'COMMON.IOUNITS'
 567       implicit none
 568       real(kind=8),dimension(:),allocatable :: prob !(maxgr)
 569       real(kind=8) :: emin,ene,en1,sumprob
 570       integer :: igr,i,ii,li1,li2,ligr,ico,jco,ind1,ind2
 571       integer :: jgr,li,nco,ib,ncon,icut
 572 !
 573 ! Compute free energies of clusters
 574 !
 575       allocate(prob(maxgr))
 576
 577       do igr=1,ngr
 578       emin=totfree(nconf(igr,1))
 579       totfree_gr(igr)=1.0d0
 580       do i=2,licz(igr)
 581         ii=nconf(igr,i)
 582         totfree_gr(igr)=totfree_gr(igr)+dexp(-totfree(ii)+emin)
 583       enddo
 584 !      write (iout,*) "igr",igr," totfree",emin,
 585 !     &    " totfree_gr",totfree_gr(igr)
 586       totfree_gr(igr)=emin-dlog(totfree_gr(igr))
 587 !      write (iout,*) igr," efree",totfree_gr(igr)/beta_h(ib)
 588       enddo
 589 !
 590 !  SORT CONFORMATIONS IN GROUPS ACC. TO ENERGY
 591 !
 592       DO 16 IGR=1,NGR
 593       LIGR=LICZ(IGR)
 594       DO 17 ICO=1,LIGR-1
 595       IND1=NCONF(IGR,ICO)
 596       ENE=totfree(IND1)
 597       DO 18 JCO=ICO+1,LIGR
 598       IND2=NCONF(IGR,JCO)
 599       EN1=totfree(IND2)
 600       IF (EN1.LT.ENE) THEN
 601         NCONF(IGR,ICO)=IND2
 602         NCONF(IGR,JCO)=IND1
 603         IND1=IND2
 604         ENE=EN1
 605       ENDIF
 606    18 CONTINUE
 607    17 CONTINUE
 608    16 CONTINUE
 609 !
 610 !  SORT GROUPS
 611 !
 612       DO 71 IGR=1,NGR
 613       ENE=totfree_gr(IGR)
 614       DO 72 JGR=IGR+1,NGR
 615       EN1=totfree_gr(JGR)
 616       IF (EN1.LT.ENE) THEN
 617         LI1=LICZ(IGR)
 618         LI2=LICZ(JGR)
 619         LI=MAX0(LI1,LI2)
 620         DO 73 I=1,LI
 621         NCO=NCONF(IGR,I)
 622         NCONF(IGR,I)=NCONF(JGR,I)
 623         NCONF(JGR,I)=NCO
 624    73   CONTINUE
 625         totfree_gr(igr)=en1
 626         totfree_gr(jgr)=ene
 627         ENE=EN1
 628         LICZ(IGR)=LI2
 629         LICZ(JGR)=LI1
 630       ENDIF
 631    72 CONTINUE
 632    71 CONTINUE
 633       write (iout,'("Free energies and probabilities of clusters at",f6.1," K")')&
 634            1.0d0/(1.987d-3*beta_h(ib)) !'
 635       prob(1)=1.0d0
 636       sumprob=1.0d0
 637       do i=2,ngr
 638         prob(i)=dexp(-(totfree_gr(i)-totfree_gr(1)))
 639         sumprob=sumprob+prob(i)
 640       enddo
 641       do i=1,ngr
 642         prob(i)=prob(i)/sumprob
 643       enddo
 644       sumprob=0.0d0
 645       write (iout,'("clust   efree    prob sumprob")')
 646       do i=1,ngr
 647         sumprob=sumprob+prob(i)
 648         write (iout,'(i5,f8.1,2f8.5)') i,totfree_gr(i)/beta_h(ib),&
 649           prob(i),sumprob
 650       enddo
 651       DO 81 IGR=1,NGR
 652       LI=LICZ(IGR)
 653       DO 82 I=1,LI
 654    82 IASS(NCONF(IGR,I))=IGR
 655    81 CONTINUE
 656       if (lgrp) then
 657         do i=1,ncon
 658           iass_tot(i,icut)=iass(i)
 659 !          write (iout,*) icut,i,iass(i),iass_tot(i,icut)
 660         enddo
 661       endif
 662       RETURN
 663       END SUBROUTINE SRTCLUST
 664 !------------------------------------------------------------------------------
 665 !------------------------------------------------------------------------------