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