source/cluster/main_clust.F

   1 C
   2 C Program to cluster united-residue MCM results.
   3 C
   4       implicit none
   5       include 'DIMENSIONS'
   6       include 'sizesclu.dat'
   7 #ifdef MPI
   8       include "mpif.h"
   9       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
  10       include "COMMON.MPI"
  11 #endif
  12       include 'COMMON.TIME1'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.GEO'
  16       include 'COMMON.HEADER'
  17       include 'COMMON.CONTROL'
  18       include 'COMMON.CHAIN'
  19       include 'COMMON.VAR'
  20       include 'COMMON.CLUSTER'
  21       include 'COMMON.IOUNITS'
  22       include 'COMMON.FREE'
  23       logical printang(max_cut)
  24       integer printpdb(max_cut)
  25       integer printmol2(max_cut)
  26       character*240 lineh
  27       REAL CRIT(maxconf),MEMBR(maxconf)
  28       REAL CRITVAL(maxconf-1)
  29       INTEGER IA(maxconf),IB(maxconf)
  30       INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
  31       INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
  32       integer nn,ndis
  33       real*4 DISNN
  34       DIMENSION NN(maxconf),DISNN(maxconf)
  35       LOGICAL FLAG(maxconf)
  36       integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon,
  37      & it,ncon_work,ind1,kkk
  38       double precision t1,t2,tcpu,difconf
  39
  40       double precision varia(maxvar)
  41       double precision hrtime,mintime,sectime
  42       logical eof
  43 #ifdef MPI
  44       call MPI_Init( IERROR )
  45       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
  46       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
  47       Master = 0
  48       if (ierror.gt.0) then
  49         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
  50         call mpi_finalize(ierror)
  51         stop
  52       endif
  53       if (nprocs.gt.MaxProcs+1) then
  54         write (2,*) "Error - too many processors",
  55      &   nprocs,MaxProcs+1
  56         write (2,*) "Increase MaxProcs and recompile"
  57         call MPI_Finalize(IERROR)
  58         stop
  59       endif
  60 #endif
  61
  62       call initialize
  63       call openunits
  64       call read_control
  65       call parmread
  66 !      call read_control
  67       call molread
  68 c      if (refstr) call read_ref_structure(*30)
  69       do i=1,nres
  70         phi(i)=0.0D0
  71         theta(i)=0.0D0
  72         alph(i)=0.0D0
  73         omeg(i)=0.0D0
  74       enddo
  75 c      write (iout,*) "Before permut"
  76 c       write (iout,*) "symetr", symetr
  77 c      call flush(iout)
  78       call permut(symetr)
  79 c      write (iout,*) "after permut"
  80 c      call flush(iout)
  81       print *,'MAIN: nnt=',nnt,' nct=',nct
  82
  83         if (nclust.gt.0) then
  84           PRINTANG(1)=.TRUE.
  85           PRINTPDB(1)=outpdb
  86           printmol2(1)=outmol2
  87           ncut=0
  88       else
  89       DO I=1,NCUT
  90         PRINTANG(I)=.FALSE.
  91         PRINTPDB(I)=0
  92         printmol2(i)=0
  93         IF (RCUTOFF(I).LT.0.0) THEN
  94           RCUTOFF(I)=ABS(RCUTOFF(I))
  95           PRINTANG(I)=.TRUE.
  96           PRINTPDB(I)=outpdb
  97           printmol2(i)=outmol2
  98         ENDIF
  99       ENDDO
 100       endif
 101       if (ncut.gt.0) then
 102       write (iout,*) 'Number of cutoffs:',NCUT
 103       write (iout,*) 'Cutoff values:'
 104       DO ICUT=1,NCUT
 105         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
 106      &    printpdb(icut),printmol2(icut)
 107       ENDDO
 108        else if (nclust.gt.0) then
 109        write (iout,'("Number of clusters requested",i5)') nclust
 110        else
 111       if (me.eq.Master)
 112      & write (iout,*) "ERROR: Either nclust or ncut must be >0"
 113        stop
 114        endif
 115       DO I=1,NRES-3
 116         MULT(I)=1
 117       ENDDO
 118       do i=1,maxconf
 119         list_conf(i)=i
 120       enddo
 121       call read_coords(ncon,*20)
 122       write (iout,*) 'from read_coords: ncon',ncon
 123
 124       write (iout,*) "nT",nT
 125       do iT=1,nT
 126       write (iout,*) "iT",iT
 127 #ifdef MPI
 128       call work_partition(.true.,ncon)
 129 #endif
 130
 131       call probabl(iT,ncon_work,ncon,*20)
 132       write(iout,*) "after probabl"
 133       if (ncon_work.lt.2) then
 134         write (iout,*) "Too few conformations; clustering skipped"
 135         exit
 136       endif
 137 #ifdef MPI
 138       ndis=ncon_work*(ncon_work-1)/2
 139       call work_partition(.true.,ndis)
 140 #endif
 141       write(iout,*) "after work partition"
 142
 143       DO I=1,NCON_work
 144         ICC(I)=I
 145       ENDDO
 146       WRITE (iout,'(A80)') TITEL
 147       t1=tcpu()
 148 C
 149 C CALCULATE DISTANCES
 150 C
 151       call daread_ccoords(1,ncon_work)
 152       write(iout,*) "after daread_ccords"
 153       ind1=0
 154       DO I=1,NCON_work-1
 155         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 156         do k=1,2*nres
 157           do l=1,3
 158             c(l,k)=allcart(l,k,i)
 159           enddo
 160         enddo
 161         kkk=1
 162         do k=1,nres
 163           do l=1,3
 164             cref(l,k,kkk)=c(l,k)
 165           enddo
 166         enddo
 167         DO J=I+1,NCON_work
 168           IND=IOFFSET(NCON_work,I,J)
 169 #ifdef MPI
 170           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
 171 #endif
 172           ind1=ind1+1
 173           DISS(IND1)=DIFCONF(I,J)
 174 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 175 #ifdef MPI
 176           endif
 177 #endif
 178         ENDDO
 179       ENDDO
 180       t2=tcpu()
 181       WRITE (iout,'(/a,1pe14.5,a/)')
 182      & 'Time for distance calculation:',T2-T1,' sec.'
 183       t1=tcpu()
 184       PRINT '(a)','End of distance computation'
 185
 186 #ifdef MPI
 187       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
 188      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 189       if (me.eq.master) then
 190 #endif
 191       open(80,file='/tmp/distance',form='unformatted')
 192       do i=1,ndis
 193         write(80) diss(i)
 194       enddo
 195       if (punch_dist) then
 196         do i=1,ncon_work-1
 197           do j=i+1,ncon_work
 198             IND=IOFFSET(NCON,I,J)
 199             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
 200      &        energy(j)-energy(i)
 201           enddo
 202         enddo
 203       endif
 204 C
 205 C Print out the RMS deviation matrix.
 206 C
 207       if (print_dist) CALL DISTOUT(NCON_work)
 208 C
 209 C  call hierarchical clustering HC from F. Murtagh
 210 C
 211       N=NCON_work
 212       LEN = (N*(N-1))/2
 213       write(iout,*) "-------------------------------------------"
 214       write(iout,*) "HIERARCHICAL CLUSTERING using"
 215       if (iopt.eq.1) then
 216         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 217       elseif (iopt.eq.2) then
 218         write(iout,*) "SINGLE LINK METHOD"
 219       elseif (iopt.eq.3) then
 220         write(iout,*) "COMPLETE LINK METHOD"
 221       elseif (iopt.eq.4) then
 222         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 223       elseif (iopt.eq.5) then
 224         write(iout,*) "MCQUITTY'S METHOD"
 225       elseif (iopt.eq.6) then
 226         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 227       elseif (iopt.eq.7) then
 228         write(iout,*) "CENTROID METHOD"
 229       else
 230         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 231         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 232         stop
 233       endif
 234       write(iout,*)
 235       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 236       write(iout,*) "February 1986"
 237       write(iout,*) "References:"
 238       write(iout,*) "1. Multidimensional clustering algorithms"
 239       write(iout,*) "   Fionn Murtagh"
 240       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 241       write(iout,*) "2. Multivariate data analysis"
 242       write(iout,*) "   Fionn Murtagh and Andre Heck"
 243       write(iout,*) "   Kluwer Academic Publishers, 1987"
 244       write(iout,*) "-------------------------------------------"
 245       write(iout,*)
 246
 247 #ifdef DEBUG
 248       write (iout,*) "The TOTFREE array"
 249       do i=1,ncon_work
 250         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
 251       enddo
 252 #endif
 253       call flush(iout)
 254       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 255       LEV = N-1
 256       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 257       if (lev.lt.2) then
 258         write (iout,*) "Too few conformations to cluster."
 259         goto 192
 260       endif
 261       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 262 !c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 263 !c 3/3/16 AL: added explicit number of cluters
 264       if (nclust.gt.0) then
 265         is=nclust-1
 266         ie=nclust-1
 267         icut=1
 268       else
 269         is=1
 270         ie=lev-1
 271       endif
 272
 273       do i=1,maxgr
 274         licz(i)=0
 275       enddo
 276       icut=1
 277       i=is
 278       NGR=is+1
 279       do j=1,n
 280         licz(iclass(j,i))=licz(iclass(j,i))+1
 281         nconf(iclass(j,i),licz(iclass(j,i)))=j
 282 !c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 283 !c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 284       enddo
 285       do i=1,lev-1
 286
 287          idum=lev-i
 288          DO L=1,LEV
 289             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 290          ENDDO
 291  190     IDUM=L
 292          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 293      &    " icut",icut," cutoff",rcutoff(icut)
 294          IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 295           if (nclust.le.0) WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 296           write (iout,'(a,f8.2)') 'Maximum distance found:',
 297      &              CRITVAL(IDUM)
 298           CALL SRTCLUST(ICUT,ncon_work,iT)
 299           CALL TRACK(ICUT)
 300           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 301           icut=icut+1
 302           if (icut.gt.ncut) goto 191
 303          ENDIF
 304          NGR=i+1
 305          do l=1,maxgr
 306           licz(l)=0
 307          enddo
 308          do j=1,n
 309           licz(iclass(j,i))=licz(iclass(j,i))+1
 310           nconf(iclass(j,i),licz(iclass(j,i)))=j
 311 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 312 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 313 cd          print *,j,iclass(j,i),
 314 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 315          enddo
 316       enddo
 317  191  continue
 318 C
 319       if (plot_tree) then
 320         CALL WRITRACK
 321         CALL PLOTREE
 322       endif
 323 C
 324       t2=tcpu()
 325       WRITE (iout,'(/a,1pe14.5,a/)')
 326      & 'Total time for clustering:',T2-T1,' sec.'
 327
 328 #ifdef MPI
 329       endif
 330 #endif
 331  192  continue
 332       enddo
 333 C
 334       close(icbase,status="delete")
 335 #ifdef MPI
 336       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 337 #endif
 338       stop '********** Program terminated normally.'
 339    20 write (iout,*) "Error reading coordinates"
 340 #ifdef MPI
 341       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 342 #endif
 343       stop
 344    30 write (iout,*) "Error reading reference structure"
 345 #ifdef MPI
 346       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 347 #endif
 348       stop
 349       end
 350 c---------------------------------------------------------------------------
 351       double precision function difconf(icon,jcon)
 352       implicit none
 353       include 'DIMENSIONS'
 354       include 'sizesclu.dat'
 355       include 'COMMON.CONTROL'
 356       include 'COMMON.CLUSTER'
 357       include 'COMMON.CHAIN'
 358       include 'COMMON.INTERACT'
 359       include 'COMMON.VAR'
 360       include 'COMMON.IOUNITS'
 361       logical non_conv
 362       double precision przes(3),obrot(3,3)
 363       double precision xx(3,maxres2),yy(3,maxres2)
 364       integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
 365       integer iaperm,ibezperm,run
 366       double precision rms,rmsmina
 367 c      write (iout,*) "tu dochodze"
 368       rmsmina=10d10
 369       nperm=1
 370       do i=1,symetr
 371       nperm=i*nperm
 372       enddo
 373 c      write (iout,*) "nperm",nperm
 374       call permut(symetr)
 375 c      write (iout,*) "tabperm", tabperm(1,1)
 376       do kkk=1,nperm
 377       if (lside) then
 378         ii=0
 379         chalen=int((nend-nstart+2)/symetr)
 380         do run=1,symetr
 381          do i=nstart,(nstart+chalen-1)
 382           zzz=tabperm(kkk,run)
 383 c          write (iout,*) "tutaj",zzz
 384           ii=ii+1
 385           iaperm=(zzz-1)*chalen+i
 386           ibezperm=(run-1)*chalen+i
 387           do j=1,3
 388             xx(j,ii)=allcart(j,iaperm,jcon)
 389             yy(j,ii)=cref(j,ibezperm,kkk)
 390           enddo
 391          enddo
 392         enddo
 393         do run=1,symetr
 394          do i=nstart,(nstart+chalen-1)
 395           zzz=tabperm(kkk,run)
 396           ii=ii+1
 397           iaperm=(zzz-1)*chalen+i
 398           ibezperm=(run-1)*chalen+i
 399 c          if (itype(i).ne.10) then
 400             ii=ii+1
 401             do j=1,3
 402               xx(j,ii)=allcart(j,iaperm+nres,jcon)
 403               yy(j,ii)=cref(j,ibezperm+nres,kkk)
 404             enddo
 405            enddo
 406 c          endif
 407         enddo
 408         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 409       else
 410         chalen=int((nct-nnt+2)/symetr)
 411         do run=1,symetr
 412          do i=nnt,(nnt+chalen-1)
 413           zzz=tabperm(kkk,run)
 414 c           write (iout,*) "tu szukaj", zzz,run,kkk
 415           iaperm=(zzz-1)*chalen+i
 416           ibezperm=(run-1)*chalen+i
 417 c        do i=nnt,nct
 418           do j=1,3
 419             c(j,i)=allcart(j,iaperm,jcon)
 420           enddo
 421          enddo
 422         enddo
 423         call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,
 424      &       przes,
 425      &       obrot,non_conv)
 426       endif
 427       if (rms.lt.0.0) then
 428         print *,'error, rms^2 = ',rms,icon,jcon
 429         stop
 430       endif
 431       if (non_conv) print *,non_conv,icon,jcon
 432       if (rmsmina.gt.rms) rmsmina=rms
 433       enddo
 434       difconf=dsqrt(rmsmina)
 435       RETURN
 436       END
 437 C------------------------------------------------------------------------------
 438       subroutine distout(ncon)
 439       implicit none
 440       include 'DIMENSIONS'
 441       include 'sizesclu.dat'
 442       integer ncol,ncon
 443       parameter (ncol=10)
 444       include 'COMMON.IOUNITS'
 445       include 'COMMON.CLUSTER'
 446       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 447       real*4 b
 448       dimension b(ncol)
 449       write (iout,'(a)') 'The distance matrix'
 450       do 1 i=1,ncon,ncol
 451       nlim=min0(i+ncol-1,ncon)
 452       write (iout,1000) (k,k=i,nlim)
 453       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 454  1000 format (/8x,10(i4,3x))
 455  1020 format (/1x,80(1h-)/)
 456       do 2 j=i,ncon
 457       jlim=min0(j,nlim)
 458       if (jlim.eq.j) then
 459         b(jlim-i+1)=0.0d0
 460         jlim1=jlim-1
 461       else
 462         jlim1=jlim
 463       endif
 464       do 3 k=i,jlim1
 465        if (j.lt.k) then
 466           IND=IOFFSET(NCON,j,k)
 467        else
 468           IND=IOFFSET(NCON,k,j)
 469        endif
 470     3  b(k-i+1)=diss(IND)
 471       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 472     2 continue
 473     1 continue
 474  1010 format (i5,3x,10(f6.2,1x))
 475       return
 476       end