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