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