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