source/cluster/wham/src-M/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,kkk
  38       double precision t1,t2,tcpu,difconf
  39
  40       double precision varia(maxvar)
  41       double precision hrtime,mintime,sectime
  42       logical eof
  43 #ifdef MPI
  44       call MPI_Init( IERROR )
  45       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
  46       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
  47       Master = 0
  48       if (ierror.gt.0) then
  49         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
  50         call mpi_finalize(ierror)
  51         stop
  52       endif
  53       if (nprocs.gt.MaxProcs+1) then
  54         write (2,*) "Error - too many processors",
  55      &   nprocs,MaxProcs+1
  56         write (2,*) "Increase MaxProcs and recompile"
  57         call MPI_Finalize(IERROR)
  58         stop
  59       endif
  60 #endif
  61
  62       call initialize
  63       call openunits
  64       call parmread
  65       call read_control
  66       call molread
  67 c      if (refstr) call read_ref_structure(*30)
  68       do i=1,nres
  69         phi(i)=0.0D0
  70         theta(i)=0.0D0
  71         alph(i)=0.0D0
  72         omeg(i)=0.0D0
  73       enddo
  74 c      write (iout,*) "Before permut"
  75 c       write (iout,*) "symetr", symetr
  76 c      call flush(iout)
  77       call permut(symetr)
  78 c      write (iout,*) "after permut"
  79 c      call flush(iout)
  80       print *,'MAIN: nnt=',nnt,' nct=',nct
  81
  82       DO I=1,NCUT
  83         PRINTANG(I)=.FALSE.
  84         PRINTPDB(I)=0
  85         printmol2(i)=0
  86         IF (RCUTOFF(I).LT.0.0) THEN
  87           RCUTOFF(I)=ABS(RCUTOFF(I))
  88           PRINTANG(I)=.TRUE.
  89           PRINTPDB(I)=outpdb
  90           printmol2(i)=outmol2
  91         ENDIF
  92       ENDDO
  93       write (iout,*) 'Number of cutoffs:',NCUT
  94       write (iout,*) 'Cutoff values:'
  95       DO ICUT=1,NCUT
  96         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
  97      &    printpdb(icut),printmol2(icut)
  98       ENDDO
  99       DO I=1,NRES-3
 100         MULT(I)=1
 101       ENDDO
 102       do i=1,maxconf
 103         list_conf(i)=i
 104       enddo
 105       call read_coords(ncon,*20)
 106       write (iout,*) 'from read_coords: ncon',ncon
 107
 108       write (iout,*) "nT",nT
 109       do iT=1,nT
 110       write (iout,*) "iT",iT
 111 #ifdef MPI
 112       call work_partition(.true.,ncon)
 113 #endif
 114
 115       call probabl(iT,ncon_work,ncon,*20)
 116
 117       if (ncon_work.lt.2) then
 118         write (iout,*) "Too few conformations; clustering skipped"
 119         exit
 120       endif
 121 #ifdef MPI
 122       ndis=ncon_work*(ncon_work-1)/2
 123       call work_partition(.true.,ndis)
 124 #endif
 125
 126       DO I=1,NCON_work
 127         ICC(I)=I
 128       ENDDO
 129       WRITE (iout,'(A80)') TITEL
 130       t1=tcpu()
 131 C
 132 C CALCULATE DISTANCES
 133 C
 134       call daread_ccoords(1,ncon_work)
 135       ind1=0
 136       DO I=1,NCON_work-1
 137         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 138         do k=1,2*nres
 139           do l=1,3
 140             c(l,k)=allcart(l,k,i)
 141           enddo
 142         enddo
 143         kkk=1
 144         do k=1,nres
 145           do l=1,3
 146             cref(l,k,kkk)=c(l,k)
 147           enddo
 148         enddo
 149         DO J=I+1,NCON_work
 150           IND=IOFFSET(NCON_work,I,J)
 151 #ifdef MPI
 152           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
 153 #endif
 154           ind1=ind1+1
 155           DISS(IND1)=DIFCONF(I,J)
 156 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 157 #ifdef MPI
 158           endif
 159 #endif
 160         ENDDO
 161       ENDDO
 162       t2=tcpu()
 163       WRITE (iout,'(/a,1pe14.5,a/)')
 164      & 'Time for distance calculation:',T2-T1,' sec.'
 165       t1=tcpu()
 166       PRINT '(a)','End of distance computation'
 167
 168 #ifdef MPI
 169       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
 170      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 171       if (me.eq.master) then
 172 #endif
 173       open(80,file='/tmp/distance',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       call flush(iout)
 236       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 237       LEV = N-1
 238       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 239       if (lev.lt.2) then
 240         write (iout,*) "Too few conformations to cluster."
 241         goto 192
 242       endif
 243       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 244 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 245
 246       do i=1,maxgr
 247         licz(i)=0
 248       enddo
 249       icut=1
 250       i=1
 251       NGR=i+1
 252       do j=1,n
 253         licz(iclass(j,i))=licz(iclass(j,i))+1
 254         nconf(iclass(j,i),licz(iclass(j,i)))=j
 255 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 256 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 257       enddo
 258       do i=1,lev-1
 259
 260          idum=lev-i
 261          DO L=1,LEV
 262             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 263          ENDDO
 264  190     IDUM=L
 265          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 266      &    " icut",icut," cutoff",rcutoff(icut)
 267          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 268           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 269           write (iout,'(a,f8.2)') 'Maximum distance found:',
 270      &              CRITVAL(IDUM)
 271           CALL SRTCLUST(ICUT,ncon_work,iT)
 272           CALL TRACK(ICUT)
 273           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 274           icut=icut+1
 275           if (icut.gt.ncut) goto 191
 276          ENDIF
 277          NGR=i+1
 278          do l=1,maxgr
 279           licz(l)=0
 280          enddo
 281          do j=1,n
 282           licz(iclass(j,i))=licz(iclass(j,i))+1
 283           nconf(iclass(j,i),licz(iclass(j,i)))=j
 284 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 285 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 286 cd          print *,j,iclass(j,i),
 287 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 288          enddo
 289       enddo
 290  191  continue
 291 C
 292       if (plot_tree) then
 293         CALL WRITRACK
 294         CALL PLOTREE
 295       endif
 296 C
 297       t2=tcpu()
 298       WRITE (iout,'(/a,1pe14.5,a/)')
 299      & 'Total time for clustering:',T2-T1,' sec.'
 300
 301 #ifdef MPI
 302       endif
 303 #endif
 304  192  continue
 305       enddo
 306 C
 307       close(icbase,status="delete")
 308 #ifdef MPI
 309       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 310 #endif
 311       stop '********** Program terminated normally.'
 312    20 write (iout,*) "Error reading coordinates"
 313 #ifdef MPI
 314       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 315 #endif
 316       stop
 317    30 write (iout,*) "Error reading reference structure"
 318 #ifdef MPI
 319       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 320 #endif
 321       stop
 322       end
 323 c---------------------------------------------------------------------------
 324       double precision function difconf(icon,jcon)
 325       implicit none
 326       include 'DIMENSIONS'
 327       include 'sizesclu.dat'
 328       include 'COMMON.CONTROL'
 329       include 'COMMON.CLUSTER'
 330       include 'COMMON.CHAIN'
 331       include 'COMMON.INTERACT'
 332       include 'COMMON.VAR'
 333       include 'COMMON.IOUNITS'
 334       logical non_conv
 335       double precision przes(3),obrot(3,3)
 336       double precision xx(3,maxres2),yy(3,maxres2)
 337       integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
 338       integer iaperm,ibezperm,run
 339       double precision rms,rmsmina
 340 c      write (iout,*) "tu dochodze"
 341       rmsmina=10d10
 342       nperm=1
 343       do i=1,symetr
 344       nperm=i*nperm
 345       enddo
 346 c      write (iout,*) "nperm",nperm
 347       call permut(symetr)
 348 c      write (iout,*) "tabperm", tabperm(1,1)
 349       do kkk=1,nperm
 350       if (lside) then
 351         ii=0
 352         chalen=int((nend-nstart+2)/symetr)
 353         do run=1,symetr
 354          do i=nstart,(nstart+chalen-1)
 355           zzz=tabperm(kkk,run)
 356 c          write (iout,*) "tutaj",zzz
 357           ii=ii+1
 358           iaperm=(zzz-1)*chalen+i
 359           ibezperm=(run-1)*chalen+i
 360           do j=1,3
 361             xx(j,ii)=allcart(j,iaperm,jcon)
 362             yy(j,ii)=cref(j,ibezperm,kkk)
 363           enddo
 364          enddo
 365         enddo
 366         do run=1,symetr
 367          do i=nstart,(nstart+chalen-1)
 368           zzz=tabperm(kkk,run)
 369           ii=ii+1
 370           iaperm=(zzz-1)*chalen+i
 371           ibezperm=(run-1)*chalen+i
 372 c          if (itype(i).ne.10) then
 373             ii=ii+1
 374             do j=1,3
 375               xx(j,ii)=allcart(j,iaperm+nres,jcon)
 376               yy(j,ii)=cref(j,ibezperm+nres,kkk)
 377             enddo
 378            enddo
 379 c          endif
 380         enddo
 381         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 382       else
 383         chalen=int((nct-nnt+2)/symetr)
 384         do run=1,symetr
 385          do i=nnt,(nnt+chalen-1)
 386           zzz=tabperm(kkk,run)
 387 c           write (iout,*) "tu szukaj", zzz,run,kkk
 388           iaperm=(zzz-1)*chalen+i
 389           ibezperm=(run-1)*chalen+i
 390 c        do i=nnt,nct
 391           do j=1,3
 392             c(j,i)=allcart(j,iaperm,jcon)
 393           enddo
 394          enddo
 395         enddo
 396         call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,
 397      &       przes,
 398      &       obrot,non_conv)
 399       endif
 400       if (rms.lt.0.0) then
 401         print *,'error, rms^2 = ',rms,icon,jcon
 402         stop
 403       endif
 404       if (non_conv) print *,non_conv,icon,jcon
 405       if (rmsmina.gt.rms) rmsmina=rms
 406       enddo
 407       difconf=dsqrt(rmsmina)
 408       RETURN
 409       END
 410 C------------------------------------------------------------------------------
 411       subroutine distout(ncon)
 412       implicit none
 413       include 'DIMENSIONS'
 414       include 'sizesclu.dat'
 415       integer ncol,ncon
 416       parameter (ncol=10)
 417       include 'COMMON.IOUNITS'
 418       include 'COMMON.CLUSTER'
 419       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 420       real*4 b
 421       dimension b(ncol)
 422       write (iout,'(a)') 'The distance matrix'
 423       do 1 i=1,ncon,ncol
 424       nlim=min0(i+ncol-1,ncon)
 425       write (iout,1000) (k,k=i,nlim)
 426       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 427  1000 format (/8x,10(i4,3x))
 428  1020 format (/1x,80(1h-)/)
 429       do 2 j=i,ncon
 430       jlim=min0(j,nlim)
 431       if (jlim.eq.j) then
 432         b(jlim-i+1)=0.0d0
 433         jlim1=jlim-1
 434       else
 435         jlim1=jlim
 436       endif
 437       do 3 k=i,jlim1
 438        if (j.lt.k) then
 439           IND=IOFFSET(NCON,j,k)
 440        else
 441           IND=IOFFSET(NCON,k,j)
 442        endif
 443     3  b(k-i+1)=diss(IND)
 444       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 445     2 continue
 446     1 continue
 447  1010 format (i5,3x,10(f6.2,1x))
 448       return
 449       end