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