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
  40       double precision varia(maxvar)
  41       double precision hrtime,mintime,sectime
  42       logical eof
  43       external ilen
  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       call initialize
  63       call openunits
  64       call cinfo
  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       scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
 168       open(80,file=scrachdir2d,form='unformatted')
 169       do i=1,ndis
 170         write(80) diss(i)
 171       enddo
 172       if (punch_dist) then
 173         do i=1,ncon_work-1
 174           do j=i+1,ncon_work
 175             IND=IOFFSET(NCON,I,J)
 176             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
 177      &        energy(j)-energy(i)
 178           enddo
 179         enddo
 180       endif
 181 C
 182 C Print out the RMS deviation matrix.
 183 C
 184       if (print_dist) CALL DISTOUT(NCON_work)
 185 C
 186 C  call hierarchical clustering HC from F. Murtagh
 187 C
 188       N=NCON_work
 189       LEN = (N*(N-1))/2
 190       write(iout,*) "-------------------------------------------"
 191       write(iout,*) "HIERARCHICAL CLUSTERING using"
 192       if (iopt.eq.1) then
 193         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 194       elseif (iopt.eq.2) then
 195         write(iout,*) "SINGLE LINK METHOD"
 196       elseif (iopt.eq.3) then
 197         write(iout,*) "COMPLETE LINK METHOD"
 198       elseif (iopt.eq.4) then
 199         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 200       elseif (iopt.eq.5) then
 201         write(iout,*) "MCQUITTY'S METHOD"
 202       elseif (iopt.eq.6) then
 203         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 204       elseif (iopt.eq.7) then
 205         write(iout,*) "CENTROID METHOD"
 206       else
 207         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 208         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 209         stop
 210       endif
 211       write(iout,*)
 212       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 213       write(iout,*) "February 1986"
 214       write(iout,*) "References:"
 215       write(iout,*) "1. Multidimensional clustering algorithms"
 216       write(iout,*) "   Fionn Murtagh"
 217       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 218       write(iout,*) "2. Multivariate data analysis"
 219       write(iout,*) "   Fionn Murtagh and Andre Heck"
 220       write(iout,*) "   Kluwer Academic Publishers, 1987"
 221       write(iout,*) "-------------------------------------------"
 222       write(iout,*)
 223
 224 #ifdef DEBUG
 225       write (iout,*) "The TOTFREE array"
 226       do i=1,ncon_work
 227         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
 228       enddo
 229 #endif
 230       call flush(iout)
 231       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 232       LEV = N-1
 233       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 234       if (lev.lt.2) then
 235         write (iout,*) "Too few conformations to cluster."
 236         goto 192
 237       endif
 238       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 239 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 240
 241       do i=1,maxgr
 242         licz(i)=0
 243       enddo
 244       icut=1
 245       i=1
 246       NGR=i+1
 247       do j=1,n
 248         licz(iclass(j,i))=licz(iclass(j,i))+1
 249         nconf(iclass(j,i),licz(iclass(j,i)))=j
 250 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 251 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 252       enddo
 253       do i=1,lev-1
 254
 255          idum=lev-i
 256          DO L=1,LEV
 257             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 258          ENDDO
 259  190     IDUM=L
 260          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 261      &    " icut",icut," cutoff",rcutoff(icut)
 262          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 263           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 264           write (iout,'(a,f8.2)') 'Maximum distance found:',
 265      &              CRITVAL(IDUM)
 266           CALL SRTCLUST(ICUT,ncon_work,iT)
 267           CALL TRACK(ICUT)
 268           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 269           icut=icut+1
 270           if (icut.gt.ncut) goto 191
 271          ENDIF
 272          NGR=i+1
 273          do l=1,maxgr
 274           licz(l)=0
 275          enddo
 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 cd          print *,j,iclass(j,i),
 282 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 283          enddo
 284       enddo
 285  191  continue
 286 C
 287       if (plot_tree) then
 288         CALL WRITRACK
 289         CALL PLOTREE
 290       endif
 291 C
 292       t2=tcpu()
 293       WRITE (iout,'(/a,1pe14.5,a/)')
 294      & 'Total time for clustering:',T2-T1,' sec.'
 295
 296 #ifdef MPI
 297       endif
 298 #endif
 299  192  continue
 300       enddo
 301 C
 302       close(icbase,status="delete")
 303 #ifdef MPI
 304       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 305 #endif
 306       stop '********** Program terminated normally.'
 307    20 write (iout,*) "Error reading coordinates"
 308 #ifdef MPI
 309       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 310 #endif
 311       stop
 312    30 write (iout,*) "Error reading reference structure"
 313 #ifdef MPI
 314       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 315 #endif
 316       stop
 317       end
 318 c---------------------------------------------------------------------------
 319       double precision function difconf(icon,jcon)
 320       implicit none
 321       include 'DIMENSIONS'
 322       include 'sizesclu.dat'
 323       include 'COMMON.CONTROL'
 324       include 'COMMON.CLUSTER'
 325       include 'COMMON.CHAIN'
 326       include 'COMMON.INTERACT'
 327       include 'COMMON.VAR'
 328       include 'COMMON.IOUNITS'
 329       logical non_conv
 330       double precision przes(3),obrot(3,3)
 331       double precision xx(3,maxres2),yy(3,maxres2)
 332       integer i,ii,j,icon,jcon
 333       double precision rms
 334       if (lside) then
 335         ii=0
 336         do i=nstart,nend
 337           ii=ii+1
 338           do j=1,3
 339             xx(j,ii)=allcart(j,i,jcon)
 340             yy(j,ii)=cref(j,i)
 341           enddo
 342         enddo
 343         do i=nstart,nend
 344 c          if (itype(i).ne.10) then
 345             ii=ii+1
 346             do j=1,3
 347               xx(j,ii)=allcart(j,i+nres,jcon)
 348               yy(j,ii)=cref(j,i+nres)
 349             enddo
 350 c          endif
 351         enddo
 352         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 353       else
 354         do i=nnt,nct
 355           do j=1,3
 356             c(j,i)=allcart(j,i,jcon)
 357           enddo
 358         enddo
 359         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
 360      &       obrot,non_conv)
 361       endif
 362       if (rms.lt.0.0) then
 363         print *,'error, rms^2 = ',rms,icon,jcon
 364         stop
 365       endif
 366       if (non_conv) print *,non_conv,icon,jcon
 367       difconf=dsqrt(rms)
 368       RETURN
 369       END
 370 C------------------------------------------------------------------------------
 371       subroutine distout(ncon)
 372       implicit none
 373       include 'DIMENSIONS'
 374       include 'sizesclu.dat'
 375       integer ncol,ncon
 376       parameter (ncol=10)
 377       include 'COMMON.IOUNITS'
 378       include 'COMMON.CLUSTER'
 379       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 380       real*4 b
 381       dimension b(ncol)
 382       write (iout,'(a)') 'The distance matrix'
 383       do 1 i=1,ncon,ncol
 384       nlim=min0(i+ncol-1,ncon)
 385       write (iout,1000) (k,k=i,nlim)
 386       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 387  1000 format (/8x,10(i4,3x))
 388  1020 format (/1x,80(1h-)/)
 389       do 2 j=i,ncon
 390       jlim=min0(j,nlim)
 391       if (jlim.eq.j) then
 392         b(jlim-i+1)=0.0d0
 393         jlim1=jlim-1
 394       else
 395         jlim1=jlim
 396       endif
 397       do 3 k=i,jlim1
 398        if (j.lt.k) then
 399           IND=IOFFSET(NCON,j,k)
 400        else
 401           IND=IOFFSET(NCON,k,j)
 402        endif
 403     3  b(k-i+1)=diss(IND)
 404       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 405     2 continue
 406     1 continue
 407  1010 format (i5,3x,10(f6.2,1x))
 408       return
 409       end