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