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