source/cluster/wham/src-M-SAXS-homology/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='/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 c 3/3/16 AL: added explicit number of cluters
 250       if (nclust.gt.0) then
 251         is=nclust-1
 252         ie=nclust-1
 253         icut=1
 254       else
 255         is=1
 256         ie=lev-1
 257       endif
 258       do i=1,maxgr
 259         licz(i)=0
 260       enddo
 261       icut=1
 262       i=is
 263       NGR=is+1
 264       do j=1,n
 265         licz(iclass(j,i))=licz(iclass(j,i))+1
 266         nconf(iclass(j,i),licz(iclass(j,i)))=j
 267 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 268 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 269       enddo
 270 c      do i=1,lev-1
 271       do i=is,ie
 272          idum=lev-i
 273          DO L=1,LEV
 274             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 275          ENDDO
 276  190     IDUM=L
 277          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 278      &    " icut",icut," cutoff",rcutoff(icut)
 279          IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 280          if (nclust.le.0)
 281      &    WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 282           write (iout,'(a,f8.2)') 'Maximum distance found:',
 283      &              CRITVAL(IDUM)
 284           CALL SRTCLUST(ICUT,ncon_work,iT)
 285           CALL TRACK(ICUT)
 286           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 287           icut=icut+1
 288           if (icut.gt.ncut) goto 191
 289          ENDIF
 290          NGR=i+1
 291          do l=1,maxgr
 292           licz(l)=0
 293          enddo
 294          do j=1,n
 295           licz(iclass(j,i))=licz(iclass(j,i))+1
 296           nconf(iclass(j,i),licz(iclass(j,i)))=j
 297 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 298 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 299 cd          print *,j,iclass(j,i),
 300 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 301          enddo
 302       enddo
 303  191  continue
 304 C
 305       if (plot_tree) then
 306         CALL WRITRACK
 307         CALL PLOTREE
 308       endif
 309 C
 310       t2=tcpu()
 311       WRITE (iout,'(/a,1pe14.5,a/)')
 312      & 'Total time for clustering:',T2-T1,' sec.'
 313
 314 #ifdef MPI
 315       endif
 316 #endif
 317  192  continue
 318       enddo
 319 C
 320       close(icbase,status="delete")
 321 #ifdef MPI
 322       call MPI_Finalize(IERROR)
 323 #endif
 324       stop '********** Program terminated normally.'
 325    20 write (iout,*) "Error reading coordinates"
 326 #ifdef MPI
 327       call MPI_Finalize(IERROR)
 328 #endif
 329       stop
 330    30 write (iout,*) "Error reading reference structure"
 331 #ifdef MPI
 332       call MPI_Finalize(IERROR)
 333 #endif
 334       stop
 335       end
 336 c---------------------------------------------------------------------------
 337       double precision function difconf(icon,jcon)
 338       implicit none
 339       include 'DIMENSIONS'
 340       include 'sizesclu.dat'
 341       include 'COMMON.CONTROL'
 342       include 'COMMON.CLUSTER'
 343       include 'COMMON.CHAIN'
 344       include 'COMMON.INTERACT'
 345       include 'COMMON.VAR'
 346       include 'COMMON.IOUNITS'
 347       integer ipermmin
 348       double precision przes(3),obrot(3,3)
 349       double precision rmscalc
 350       integer icon,jcon,k,l
 351 c      write (iout,*) "DIFCONF: ICON",icon," JCON",jcon
 352       do k=1,2*nres
 353         do l=1,3
 354           cref(l,k)=allcart(l,k,icon)
 355           c(l,k)=allcart(l,k,jcon)
 356         enddo
 357       enddo
 358       difconf=rmscalc(c(1,1),cref(1,1),przes,obrot,ipermmin)
 359       RETURN
 360       END
 361 C------------------------------------------------------------------------------
 362       subroutine distout(ncon)
 363       implicit none
 364       include 'DIMENSIONS'
 365       include 'sizesclu.dat'
 366       integer ncol,ncon
 367       parameter (ncol=10)
 368       include 'COMMON.IOUNITS'
 369       include 'COMMON.CLUSTER'
 370       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 371       real*4 b
 372       dimension b(ncol)
 373       write (iout,'(a)') 'The distance matrix'
 374       do 1 i=1,ncon,ncol
 375       nlim=min0(i+ncol-1,ncon)
 376       write (iout,1000) (k,k=i,nlim)
 377       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 378  1000 format (/8x,10(i4,3x))
 379  1020 format (/1x,80(1h-)/)
 380       do 2 j=i,ncon
 381       jlim=min0(j,nlim)
 382       if (jlim.eq.j) then
 383         b(jlim-i+1)=0.0d0
 384         jlim1=jlim-1
 385       else
 386         jlim1=jlim
 387       endif
 388       do 3 k=i,jlim1
 389        if (j.lt.k) then
 390           IND=IOFFSET(NCON,j,k)
 391        else
 392           IND=IOFFSET(NCON,k,j)
 393        endif
 394     3  b(k-i+1)=diss(IND)
 395       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 396     2 continue
 397     1 continue
 398  1010 format (i5,3x,10(f6.2,1x))
 399       return
 400       end