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