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