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
 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
 142       DO I=1,NCON_work
 143         ICC(I)=I
 144       ENDDO
 145       WRITE (iout,'(A80)') TITEL
 146       t1=tcpu()
 147 C
 148 C CALCULATE DISTANCES
 149 C
 150       call daread_ccoords(1,ncon_work)
 151       ind1=0
 152       DO I=1,NCON_work-1
 153         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 154         do k=1,2*nres
 155           do l=1,3
 156             c(l,k)=allcart(l,k,i)
 157           enddo
 158         enddo
 159         kkk=1
 160         do k=1,nres
 161           do l=1,3
 162             cref(l,k,kkk)=c(l,k)
 163           enddo
 164         enddo
 165         DO J=I+1,NCON_work
 166           IND=IOFFSET(NCON_work,I,J)
 167 #ifdef MPI
 168           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
 169 #endif
 170           ind1=ind1+1
 171           DISS(IND1)=DIFCONF(I,J)
 172 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 173 #ifdef MPI
 174           endif
 175 #endif
 176         ENDDO
 177       ENDDO
 178       t2=tcpu()
 179       WRITE (iout,'(/a,1pe14.5,a/)')
 180      & 'Time for distance calculation:',T2-T1,' sec.'
 181       t1=tcpu()
 182       PRINT '(a)','End of distance computation'
 183
 184 #ifdef MPI
 185       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
 186      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 187       if (me.eq.master) then
 188 #endif
 189       open(80,file='/tmp/distance',form='unformatted')
 190       do i=1,ndis
 191         write(80) diss(i)
 192       enddo
 193       if (punch_dist) then
 194         do i=1,ncon_work-1
 195           do j=i+1,ncon_work
 196             IND=IOFFSET(NCON,I,J)
 197             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
 198      &        energy(j)-energy(i)
 199           enddo
 200         enddo
 201       endif
 202 C
 203 C Print out the RMS deviation matrix.
 204 C
 205       if (print_dist) CALL DISTOUT(NCON_work)
 206 C
 207 C  call hierarchical clustering HC from F. Murtagh
 208 C
 209       N=NCON_work
 210       LEN = (N*(N-1))/2
 211       write(iout,*) "-------------------------------------------"
 212       write(iout,*) "HIERARCHICAL CLUSTERING using"
 213       if (iopt.eq.1) then
 214         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 215       elseif (iopt.eq.2) then
 216         write(iout,*) "SINGLE LINK METHOD"
 217       elseif (iopt.eq.3) then
 218         write(iout,*) "COMPLETE LINK METHOD"
 219       elseif (iopt.eq.4) then
 220         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 221       elseif (iopt.eq.5) then
 222         write(iout,*) "MCQUITTY'S METHOD"
 223       elseif (iopt.eq.6) then
 224         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 225       elseif (iopt.eq.7) then
 226         write(iout,*) "CENTROID METHOD"
 227       else
 228         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 229         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 230         stop
 231       endif
 232       write(iout,*)
 233       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 234       write(iout,*) "February 1986"
 235       write(iout,*) "References:"
 236       write(iout,*) "1. Multidimensional clustering algorithms"
 237       write(iout,*) "   Fionn Murtagh"
 238       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 239       write(iout,*) "2. Multivariate data analysis"
 240       write(iout,*) "   Fionn Murtagh and Andre Heck"
 241       write(iout,*) "   Kluwer Academic Publishers, 1987"
 242       write(iout,*) "-------------------------------------------"
 243       write(iout,*)
 244
 245 #ifdef DEBUG
 246       write (iout,*) "The TOTFREE array"
 247       do i=1,ncon_work
 248         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
 249       enddo
 250 #endif
 251       call flush(iout)
 252       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 253       LEV = N-1
 254       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 255       if (lev.lt.2) then
 256         write (iout,*) "Too few conformations to cluster."
 257         goto 192
 258       endif
 259       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 260 !c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 261 !c 3/3/16 AL: added explicit number of cluters
 262       if (nclust.gt.0) then
 263         is=nclust-1
 264         ie=nclust-1
 265         icut=1
 266       else
 267         is=1
 268         ie=lev-1
 269       endif
 270
 271       do i=1,maxgr
 272         licz(i)=0
 273       enddo
 274       icut=1
 275       i=is
 276       NGR=is+1
 277       do j=1,n
 278         licz(iclass(j,i))=licz(iclass(j,i))+1
 279         nconf(iclass(j,i),licz(iclass(j,i)))=j
 280 !c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 281 !c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 282       enddo
 283       do i=1,lev-1
 284
 285          idum=lev-i
 286          DO L=1,LEV
 287             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 288          ENDDO
 289  190     IDUM=L
 290          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 291      &    " icut",icut," cutoff",rcutoff(icut)
 292          IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 293           if (nclust.le.0) WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 294           write (iout,'(a,f8.2)') 'Maximum distance found:',
 295      &              CRITVAL(IDUM)
 296           CALL SRTCLUST(ICUT,ncon_work,iT)
 297           CALL TRACK(ICUT)
 298           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 299           icut=icut+1
 300           if (icut.gt.ncut) goto 191
 301          ENDIF
 302          NGR=i+1
 303          do l=1,maxgr
 304           licz(l)=0
 305          enddo
 306          do j=1,n
 307           licz(iclass(j,i))=licz(iclass(j,i))+1
 308           nconf(iclass(j,i),licz(iclass(j,i)))=j
 309 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 310 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 311 cd          print *,j,iclass(j,i),
 312 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 313          enddo
 314       enddo
 315  191  continue
 316 C
 317       if (plot_tree) then
 318         CALL WRITRACK
 319         CALL PLOTREE
 320       endif
 321 C
 322       t2=tcpu()
 323       WRITE (iout,'(/a,1pe14.5,a/)')
 324      & 'Total time for clustering:',T2-T1,' sec.'
 325
 326 #ifdef MPI
 327       endif
 328 #endif
 329  192  continue
 330       enddo
 331 C
 332       close(icbase,status="delete")
 333 #ifdef MPI
 334       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 335 #endif
 336       stop '********** Program terminated normally.'
 337    20 write (iout,*) "Error reading coordinates"
 338 #ifdef MPI
 339       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 340 #endif
 341       stop
 342    30 write (iout,*) "Error reading reference structure"
 343 #ifdef MPI
 344       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 345 #endif
 346       stop
 347       end
 348 c---------------------------------------------------------------------------
 349       double precision function difconf(icon,jcon)
 350       implicit none
 351       include 'DIMENSIONS'
 352       include 'sizesclu.dat'
 353       include 'COMMON.CONTROL'
 354       include 'COMMON.CLUSTER'
 355       include 'COMMON.CHAIN'
 356       include 'COMMON.INTERACT'
 357       include 'COMMON.VAR'
 358       include 'COMMON.IOUNITS'
 359       logical non_conv
 360       double precision przes(3),obrot(3,3)
 361       double precision xx(3,maxres2),yy(3,maxres2)
 362       integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
 363       integer iaperm,ibezperm,run
 364       double precision rms,rmsmina
 365 c      write (iout,*) "tu dochodze"
 366       rmsmina=10d10
 367       nperm=1
 368       do i=1,symetr
 369       nperm=i*nperm
 370       enddo
 371 c      write (iout,*) "nperm",nperm
 372       call permut(symetr)
 373 c      write (iout,*) "tabperm", tabperm(1,1)
 374       do kkk=1,nperm
 375       if (lside) then
 376         ii=0
 377         chalen=int((nend-nstart+2)/symetr)
 378         do run=1,symetr
 379          do i=nstart,(nstart+chalen-1)
 380           zzz=tabperm(kkk,run)
 381 c          write (iout,*) "tutaj",zzz
 382           ii=ii+1
 383           iaperm=(zzz-1)*chalen+i
 384           ibezperm=(run-1)*chalen+i
 385           do j=1,3
 386             xx(j,ii)=allcart(j,iaperm,jcon)
 387             yy(j,ii)=cref(j,ibezperm,kkk)
 388           enddo
 389          enddo
 390         enddo
 391         do run=1,symetr
 392          do i=nstart,(nstart+chalen-1)
 393           zzz=tabperm(kkk,run)
 394           ii=ii+1
 395           iaperm=(zzz-1)*chalen+i
 396           ibezperm=(run-1)*chalen+i
 397 c          if (itype(i).ne.10) then
 398             ii=ii+1
 399             do j=1,3
 400               xx(j,ii)=allcart(j,iaperm+nres,jcon)
 401               yy(j,ii)=cref(j,ibezperm+nres,kkk)
 402             enddo
 403            enddo
 404 c          endif
 405         enddo
 406         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 407       else
 408         chalen=int((nct-nnt+2)/symetr)
 409         do run=1,symetr
 410          do i=nnt,(nnt+chalen-1)
 411           zzz=tabperm(kkk,run)
 412 c           write (iout,*) "tu szukaj", zzz,run,kkk
 413           iaperm=(zzz-1)*chalen+i
 414           ibezperm=(run-1)*chalen+i
 415 c        do i=nnt,nct
 416           do j=1,3
 417             c(j,i)=allcart(j,iaperm,jcon)
 418           enddo
 419          enddo
 420         enddo
 421         call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,
 422      &       przes,
 423      &       obrot,non_conv)
 424       endif
 425       if (rms.lt.0.0) then
 426         print *,'error, rms^2 = ',rms,icon,jcon
 427         stop
 428       endif
 429       if (non_conv) print *,non_conv,icon,jcon
 430       if (rmsmina.gt.rms) rmsmina=rms
 431       enddo
 432       difconf=dsqrt(rmsmina)
 433       RETURN
 434       END
 435 C------------------------------------------------------------------------------
 436       subroutine distout(ncon)
 437       implicit none
 438       include 'DIMENSIONS'
 439       include 'sizesclu.dat'
 440       integer ncol,ncon
 441       parameter (ncol=10)
 442       include 'COMMON.IOUNITS'
 443       include 'COMMON.CLUSTER'
 444       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 445       real*4 b
 446       dimension b(ncol)
 447       write (iout,'(a)') 'The distance matrix'
 448       do 1 i=1,ncon,ncol
 449       nlim=min0(i+ncol-1,ncon)
 450       write (iout,1000) (k,k=i,nlim)
 451       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 452  1000 format (/8x,10(i4,3x))
 453  1020 format (/1x,80(1h-)/)
 454       do 2 j=i,ncon
 455       jlim=min0(j,nlim)
 456       if (jlim.eq.j) then
 457         b(jlim-i+1)=0.0d0
 458         jlim1=jlim-1
 459       else
 460         jlim1=jlim
 461       endif
 462       do 3 k=i,jlim1
 463        if (j.lt.k) then
 464           IND=IOFFSET(NCON,j,k)
 465        else
 466           IND=IOFFSET(NCON,k,j)
 467        endif
 468     3  b(k-i+1)=diss(IND)
 469       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 470     2 continue
 471     1 continue
 472  1010 format (i5,3x,10(f6.2,1x))
 473       return
 474       end