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