source/cluster/wham/src-HCD-5D/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, is,ie
  38       double precision t1,t2,tcpu,difconf
  39
  40       double precision varia(maxvar)
  41       double precision hrtime,mintime,sectime
  42       logical eof
  43       integer ilen
  44       external ilen
  45 #ifdef MPI
  46       call MPI_Init( IERROR )
  47       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
  48       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
  49       Master = 0
  50       if (ierror.gt.0) then
  51         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
  52         call mpi_finalize(ierror)
  53         stop
  54       endif
  55       if (nprocs.gt.MaxProcs+1) then
  56         write (2,*) "Error - too many processors",
  57      &   nprocs,MaxProcs+1
  58         write (2,*) "Increase MaxProcs and recompile"
  59         call MPI_Finalize(IERROR)
  60         stop
  61       endif
  62 #endif
  63
  64       call initialize
  65       call openunits
  66       call cinfo
  67       call read_control
  68       call parmread
  69       call molread
  70 c      write (iout,*) "Main: refstr ",refstr
  71       if (refstr) call read_ref_structure(*30)
  72       do i=1,nres
  73         phi(i)=0.0D0
  74         theta(i)=0.0D0
  75         alph(i)=0.0D0
  76         omeg(i)=0.0D0
  77       enddo
  78       if (nclust.gt.0) then
  79         PRINTANG(1)=.TRUE.
  80         PRINTPDB(1)=outpdb
  81         printmol2(1)=outmol2
  82         ncut=0
  83       else
  84       DO I=1,NCUT
  85         PRINTANG(I)=.FALSE.
  86         PRINTPDB(I)=0
  87         printmol2(i)=0
  88         IF (RCUTOFF(I).LT.0.0) THEN
  89           RCUTOFF(I)=ABS(RCUTOFF(I))
  90           PRINTANG(I)=.TRUE.
  91           PRINTPDB(I)=outpdb
  92           printmol2(i)=outmol2
  93         ENDIF
  94       ENDDO
  95       endif
  96       if (ncut.gt.0) then
  97       write (iout,*) 'Number of cutoffs:',NCUT
  98       write (iout,*) 'Cutoff values:'
  99       DO ICUT=1,NCUT
 100         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
 101      &    printpdb(icut),printmol2(icut)
 102       ENDDO
 103       else if (nclust.gt.0) then
 104       write (iout,'("Number of clusters requested",i5)') nclust
 105       else
 106       if (me.eq.Master)
 107      & write (iout,*) "ERROR: Either nclust or ncut must be >0"
 108       stop
 109       endif
 110       DO I=1,NRES-3
 111         MULT(I)=1
 112       ENDDO
 113       do i=1,maxconf
 114         list_conf(i)=i
 115       enddo
 116       call read_coords(ncon,*20)
 117       write (iout,*) 'from read_coords: ncon',ncon
 118
 119       write (iout,*) "nT",nT
 120       do iT=1,nT
 121       write (iout,*) "Temperature",1.0d0/(beta_h(iT)*1.987D-3)
 122 #ifdef MPI
 123       call work_partition(.true.,ncon)
 124 #endif
 125       call probabl(iT,ncon_work,ncon,*20)
 126
 127       if (ncon_work.lt.2) then
 128         write (iout,*) "Too few conformations; clustering skipped"
 129         exit
 130       endif
 131 #ifdef MPI
 132       ndis=ncon_work*(ncon_work-1)/2
 133       call work_partition(.true.,ndis)
 134 #endif
 135       DO I=1,NCON_work
 136         ICC(I)=I
 137       ENDDO
 138       WRITE (iout,'(A80)') TITEL
 139       t1=tcpu()
 140 C
 141 C CALCULATE DISTANCES
 142 C
 143       call daread_ccoords(1,ncon_work)
 144       ind1=0
 145       DO I=1,NCON_work-1
 146 c        if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 147         DO J=I+1,NCON_work
 148           IND=IOFFSET(NCON_work,I,J)
 149 #ifdef MPI
 150           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
 151 #endif
 152           ind1=ind1+1
 153           DISS(IND1)=DIFCONF(I,J)
 154 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 155 #ifdef MPI
 156           endif
 157 #endif
 158         ENDDO
 159       ENDDO
 160       t2=tcpu()
 161       WRITE (iout,'(/a,1pe14.5,a/)')
 162      & 'Time for distance calculation:',T2-T1,' sec.'
 163       t1=tcpu()
 164 c      PRINT '(a)','End of distance computation'
 165
 166       scount_buf=scount(me)
 167
 168       do ijk=1, ndis
 169       diss_buf(ijk)=diss(ijk)
 170       enddo
 171
 172
 173 #ifdef MPI
 174       WRITE (iout,*) "Wchodze do call MPI_Gatherv"
 175       call MPI_Gatherv(diss_buf(1),scount_buf,MPI_REAL,diss(1),
 176      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 177       if (me.eq.master) then
 178 #endif
 179       open(80,file=scratchdir(:ilen(scratchdir))//'/distance',
 180      &   form='unformatted')
 181       do i=1,ndis
 182         write(80) diss(i)
 183       enddo
 184       if (punch_dist) then
 185         do i=1,ncon_work-1
 186           do j=i+1,ncon_work
 187             IND=IOFFSET(NCON,I,J)
 188             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
 189      &        energy(j)-energy(i)
 190           enddo
 191         enddo
 192       endif
 193 C
 194 C Print out the RMS deviation matrix.
 195 C
 196       if (print_dist) CALL DISTOUT(NCON_work)
 197 C
 198 C  call hierarchical clustering HC from F. Murtagh
 199 C
 200       N=NCON_work
 201       LEN = (N*(N-1))/2
 202       write(iout,*) "-------------------------------------------"
 203       write(iout,*) "HIERARCHICAL CLUSTERING using"
 204       if (iopt.eq.1) then
 205         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 206       elseif (iopt.eq.2) then
 207         write(iout,*) "SINGLE LINK METHOD"
 208       elseif (iopt.eq.3) then
 209         write(iout,*) "COMPLETE LINK METHOD"
 210       elseif (iopt.eq.4) then
 211         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 212       elseif (iopt.eq.5) then
 213         write(iout,*) "MCQUITTY'S METHOD"
 214       elseif (iopt.eq.6) then
 215         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 216       elseif (iopt.eq.7) then
 217         write(iout,*) "CENTROID METHOD"
 218       else
 219         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 220         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 221         stop
 222       endif
 223       write(iout,*)
 224       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 225       write(iout,*) "February 1986"
 226       write(iout,*) "References:"
 227       write(iout,*) "1. Multidimensional clustering algorithms"
 228       write(iout,*) "   Fionn Murtagh"
 229       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 230       write(iout,*) "2. Multivariate data analysis"
 231       write(iout,*) "   Fionn Murtagh and Andre Heck"
 232       write(iout,*) "   Kluwer Academic Publishers, 1987"
 233       write(iout,*) "-------------------------------------------"
 234       write(iout,*)
 235
 236 #ifdef DEBUG
 237       write (iout,*) "The TOTFREE array"
 238       do i=1,ncon_work
 239         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
 240       enddo
 241 #endif
 242       call flush(iout)
 243       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 244       LEV = N-1
 245       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 246       if (lev.lt.2) then
 247         write (iout,*) "Too few conformations to cluster."
 248         goto 192
 249       endif
 250       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 251 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 252 c 3/3/16 AL: added explicit number of cluters
 253       if (nclust.gt.0) then
 254         is=nclust-1
 255         ie=nclust-1
 256         icut=1
 257       else
 258         is=1
 259         ie=lev-1
 260       endif
 261       do i=1,maxgr
 262         licz(i)=0
 263       enddo
 264       icut=1
 265       i=is
 266       NGR=is+1
 267       do j=1,n
 268         licz(iclass(j,i))=licz(iclass(j,i))+1
 269         nconf(iclass(j,i),licz(iclass(j,i)))=j
 270 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 271 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 272       enddo
 273 c      do i=1,lev-1
 274       do i=is,ie
 275          idum=lev-i
 276          DO L=1,LEV
 277             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 278          ENDDO
 279  190     IDUM=L
 280          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 281      &    " icut",icut," cutoff",rcutoff(icut)
 282          IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 283          if (nclust.le.0)
 284      &    WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 285           write (iout,'(a,f8.2)') 'Maximum distance found:',
 286      &              CRITVAL(IDUM)
 287           CALL SRTCLUST(ICUT,ncon_work,iT)
 288           CALL TRACK(ICUT)
 289           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 290           icut=icut+1
 291           if (icut.gt.ncut) goto 191
 292          ENDIF
 293          NGR=i+1
 294          do l=1,maxgr
 295           licz(l)=0
 296          enddo
 297          do j=1,n
 298           licz(iclass(j,i))=licz(iclass(j,i))+1
 299           nconf(iclass(j,i),licz(iclass(j,i)))=j
 300 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 301 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 302 cd          print *,j,iclass(j,i),
 303 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 304          enddo
 305       enddo
 306  191  continue
 307 C
 308       if (plot_tree) then
 309         CALL WRITRACK
 310         CALL PLOTREE
 311       endif
 312 C
 313       t2=tcpu()
 314       WRITE (iout,'(/a,1pe14.5,a/)')
 315      & 'Total time for clustering:',T2-T1,' sec.'
 316
 317 #ifdef MPI
 318       endif
 319 #endif
 320  192  continue
 321       enddo
 322 C
 323       close(icbase,status="delete")
 324 #ifdef MPI
 325       call MPI_Finalize(IERROR)
 326 #endif
 327       stop '********** Program terminated normally.'
 328    20 write (iout,*) "Error reading coordinates"
 329 #ifdef MPI
 330       call MPI_Finalize(IERROR)
 331 #endif
 332       stop
 333    30 write (iout,*) "Error reading reference structure"
 334 #ifdef MPI
 335       call MPI_Finalize(IERROR)
 336 #endif
 337       stop
 338       end
 339 c---------------------------------------------------------------------------
 340       double precision function difconf(icon,jcon)
 341       implicit none
 342       include 'DIMENSIONS'
 343       include 'sizesclu.dat'
 344       include 'COMMON.CONTROL'
 345       include 'COMMON.CLUSTER'
 346       include 'COMMON.CHAIN'
 347       include 'COMMON.INTERACT'
 348       include 'COMMON.VAR'
 349       include 'COMMON.IOUNITS'
 350       integer ipermmin
 351       double precision przes(3),obrot(3,3)
 352       double precision rmscalc
 353       integer icon,jcon,k,l
 354 c      write (iout,*) "DIFCONF: ICON",icon," JCON",jcon
 355       do k=1,2*nres
 356         do l=1,3
 357           cref(l,k)=allcart(l,k,icon)
 358           c(l,k)=allcart(l,k,jcon)
 359         enddo
 360       enddo
 361       difconf=rmscalc(c(1,1),cref(1,1),przes,obrot,ipermmin)
 362       RETURN
 363       END
 364 C------------------------------------------------------------------------------
 365       subroutine distout(ncon)
 366       implicit none
 367       include 'DIMENSIONS'
 368       include 'sizesclu.dat'
 369       integer ncol,ncon
 370       parameter (ncol=10)
 371       include 'COMMON.IOUNITS'
 372       include 'COMMON.CLUSTER'
 373       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 374       real*4 b
 375       dimension b(ncol)
 376       write (iout,'(a)') 'The distance matrix'
 377       do 1 i=1,ncon,ncol
 378       nlim=min0(i+ncol-1,ncon)
 379       write (iout,1000) (k,k=i,nlim)
 380       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 381  1000 format (/8x,10(i4,3x))
 382  1020 format (/1x,80(1h-)/)
 383       do 2 j=i,ncon
 384       jlim=min0(j,nlim)
 385       if (jlim.eq.j) then
 386         b(jlim-i+1)=0.0d0
 387         jlim1=jlim-1
 388       else
 389         jlim1=jlim
 390       endif
 391       do 3 k=i,jlim1
 392        if (j.lt.k) then
 393           IND=IOFFSET(NCON,j,k)
 394        else
 395           IND=IOFFSET(NCON,k,j)
 396        endif
 397     3  b(k-i+1)=diss(IND)
 398       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 399     2 continue
 400     1 continue
 401  1010 format (i5,3x,10(f6.2,1x))
 402       return
 403       end