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