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