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