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       call initialize
  64       call openunits
  65       call cinfo
  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 #ifdef AIX
 251       call flush_(iout)
 252 #else
 253       call flush(iout)
 254 #endif
 255       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 256       LEV = N-1
 257       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 258       if (lev.lt.2) then
 259         write (iout,*) "Too few conformations to cluster."
 260         goto 192
 261       endif
 262       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 263 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 264
 265 c 3/3/16 AL: added explicit number of cluters
 266       if (nclust.gt.0) then
 267         is=nclust-1
 268         ie=nclust-1
 269         icut=1
 270       else
 271         is=1
 272         ie=lev-1
 273       endif
 274       do i=1,maxgr
 275         licz(i)=0
 276       enddo
 277       icut=1
 278       i=is
 279       NGR=is+1
 280       do j=1,n
 281         licz(iclass(j,i))=licz(iclass(j,i))+1
 282         nconf(iclass(j,i),licz(iclass(j,i)))=j
 283 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 284 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 285       enddo
 286 c      do i=1,lev-1
 287       do i=is,ie
 288          idum=lev-i
 289          DO L=1,LEV
 290             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 291          ENDDO
 292  190     IDUM=L
 293 c         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 294 c     &    " icut",icut," cutoff",rcutoff(icut)
 295          IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 296          if (nclust.le.0)
 297      &    WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 298           write (iout,'(a,f8.2)') 'Maximum distance found:',
 299      &              CRITVAL(IDUM)
 300           CALL SRTCLUST(ICUT,ncon_work,iT)
 301           CALL TRACK(ICUT)
 302           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 303           icut=icut+1
 304           if (icut.gt.ncut) goto 191
 305          ENDIF
 306          NGR=i+1
 307          do l=1,maxgr
 308           licz(l)=0
 309          enddo
 310          ii=i-is+1
 311          do j=1,n
 312           licz(iclass(j,ii))=licz(iclass(j,ii))+1
 313           nconf(iclass(j,ii),licz(iclass(j,ii)))=j
 314 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 315 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 316 cd          print *,j,iclass(j,i),
 317 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 318          enddo
 319       enddo
 320  191  continue
 321 C
 322       if (plot_tree) then
 323         CALL WRITRACK
 324         CALL PLOTREE
 325       endif
 326 C
 327       t2=tcpu()
 328       WRITE (iout,'(/a,1pe14.5,a/)')
 329      & 'Total time for clustering:',T2-T1,' sec.'
 330
 331 #ifdef MPI
 332       endif
 333 #endif
 334  192  continue
 335       enddo
 336 C
 337       close(icbase,status="delete")
 338 #ifdef MPI
 339       call MPI_Finalize(IERROR)
 340 #endif
 341       stop '********** Program terminated normally.'
 342    20 write (iout,*) "Error reading coordinates"
 343 #ifdef MPI
 344       call MPI_Finalize(IERROR)
 345 #endif
 346       stop
 347    30 write (iout,*) "Error reading reference structure"
 348 #ifdef MPI
 349       call MPI_Finalize(IERROR)
 350 #endif
 351       stop
 352       end
 353 c---------------------------------------------------------------------------
 354       double precision function difconf(icon,jcon)
 355       implicit none
 356       include 'DIMENSIONS'
 357       include 'sizesclu.dat'
 358       include 'COMMON.CONTROL'
 359       include 'COMMON.CLUSTER'
 360       include 'COMMON.CHAIN'
 361       include 'COMMON.INTERACT'
 362       include 'COMMON.VAR'
 363       include 'COMMON.IOUNITS'
 364       logical non_conv
 365       double precision przes(3),obrot(3,3)
 366       double precision xx(3,maxres2),yy(3,maxres2)
 367       integer i,ii,j,icon,jcon
 368       double precision rms
 369       if (lside) then
 370         ii=0
 371         do i=nstart,nend
 372           ii=ii+1
 373           do j=1,3
 374             xx(j,ii)=allcart(j,i,jcon)
 375             yy(j,ii)=cref(j,i)
 376           enddo
 377         enddo
 378         do i=nstart,nend
 379 c          if (itype(i).ne.10) then
 380             ii=ii+1
 381             do j=1,3
 382               xx(j,ii)=allcart(j,i+nres,jcon)
 383               yy(j,ii)=cref(j,i+nres)
 384             enddo
 385 c          endif
 386         enddo
 387         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 388       else
 389         do i=nnt,nct
 390           do j=1,3
 391             c(j,i)=allcart(j,i,jcon)
 392           enddo
 393         enddo
 394         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
 395      &       obrot,non_conv)
 396       endif
 397       if (rms.lt.0.0) then
 398         print *,'error, rms^2 = ',rms,icon,jcon
 399         stop
 400       endif
 401       if (non_conv) print *,non_conv,icon,jcon
 402       difconf=dsqrt(rms)
 403       RETURN
 404       END
 405 C------------------------------------------------------------------------------
 406       subroutine distout(ncon)
 407       implicit none
 408       include 'DIMENSIONS'
 409       include 'sizesclu.dat'
 410       integer ncol,ncon
 411       parameter (ncol=10)
 412       include 'COMMON.IOUNITS'
 413       include 'COMMON.CLUSTER'
 414       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 415       real*4 b
 416       dimension b(ncol)
 417       write (iout,'(a)') 'The distance matrix'
 418       do 1 i=1,ncon,ncol
 419       nlim=min0(i+ncol-1,ncon)
 420       write (iout,1000) (k,k=i,nlim)
 421       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 422  1000 format (/8x,10(i4,3x))
 423  1020 format (/1x,80(1h-)/)
 424       do 2 j=i,ncon
 425       jlim=min0(j,nlim)
 426       if (jlim.eq.j) then
 427         b(jlim-i+1)=0.0d0
 428         jlim1=jlim-1
 429       else
 430         jlim1=jlim
 431       endif
 432       do 3 k=i,jlim1
 433        if (j.lt.k) then
 434           IND=IOFFSET(NCON,j,k)
 435        else
 436           IND=IOFFSET(NCON,k,j)
 437        endif
 438     3  b(k-i+1)=diss(IND)
 439       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 440     2 continue
 441     1 continue
 442  1010 format (i5,3x,10(f6.2,1x))
 443       return
 444       end