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