source/cluster/wham/src-M/bakup/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 #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        write (iout,*) "symetr", symetr
  75       call permut(symetr)
  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
 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           DISS(IND1)=DIFCONF(I,J)
 151 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 152 #ifdef MPI
 153           endif
 154 #endif
 155         ENDDO
 156       ENDDO
 157       t2=tcpu()
 158       WRITE (iout,'(/a,1pe14.5,a/)')
 159      & 'Time for distance calculation:',T2-T1,' sec.'
 160       t1=tcpu()
 161       PRINT '(a)','End of distance computation'
 162
 163 #ifdef MPI
 164       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
 165      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 166       if (me.eq.master) then
 167 #endif
 168       open(80,file='/tmp/distance',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,kkk,nperm,chalen,zzz
 333       integer iaperm,ibezperm,run
 334       double precision rms,rmsmina
 335 c      write (iout,*) "tu dochodze"
 336       rmsmina=10d10
 337       nperm=1
 338       do i=1,symetr
 339       nperm=i*nperm
 340       enddo
 341 c      write (iout,*) "nperm",nperm
 342       call permut(symetr)
 343 c      write (iout,*) "tabperm", tabperm(1,1)
 344       do kkk=1,nperm
 345       if (lside) then
 346         ii=0
 347         chalen=int((nend-nstart+2)/symetr)
 348         do run=1,symetr
 349          do i=nstart,(nstart+chalen-1)
 350           zzz=tabperm(kkk,run)
 351 c          write (iout,*) "tutaj",zzz
 352           ii=ii+1
 353           iaperm=(zzz-1)*chalen+i
 354           ibezperm=(run-1)*chalen+i
 355           do j=1,3
 356             xx(j,ii)=allcart(j,iaperm,jcon)
 357             yy(j,ii)=cref(j,ibezperm)
 358           enddo
 359          enddo
 360         enddo
 361         do run=1,symetr
 362          do i=nstart,(nstart+chalen-1)
 363           zzz=tabperm(kkk,run)
 364           ii=ii+1
 365           iaperm=(zzz-1)*chalen+i
 366           ibezperm=(run-1)*chalen+i
 367 c          if (itype(i).ne.10) then
 368             ii=ii+1
 369             do j=1,3
 370               xx(j,ii)=allcart(j,iaperm+nres,jcon)
 371               yy(j,ii)=cref(j,ibezperm+nres)
 372             enddo
 373            enddo
 374 c          endif
 375         enddo
 376         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 377       else
 378         chalen=int((nct-nnt+2)/symetr)
 379         do run=1,symetr
 380          do i=nnt,(nnt+chalen-1)
 381           zzz=tabperm(kkk,run)
 382 c           write (iout,*) "tu szukaj", zzz,run,kkk
 383           iaperm=(zzz-1)*chalen+i
 384           ibezperm=(run-1)*chalen+i
 385 c        do i=nnt,nct
 386           do j=1,3
 387             c(j,i)=allcart(j,iaperm,jcon)
 388           enddo
 389          enddo
 390         enddo
 391         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
 392      &       obrot,non_conv)
 393       endif
 394       if (rms.lt.0.0) then
 395         print *,'error, rms^2 = ',rms,icon,jcon
 396         stop
 397       endif
 398       if (non_conv) print *,non_conv,icon,jcon
 399       if (rmsmina.gt.rms) rmsmina=rms
 400       enddo
 401       difconf=dsqrt(rmsmina)
 402       RETURN
 403       END
 404 C------------------------------------------------------------------------------
 405       subroutine distout(ncon)
 406       implicit none
 407       include 'DIMENSIONS'
 408       include 'sizesclu.dat'
 409       integer ncol,ncon
 410       parameter (ncol=10)
 411       include 'COMMON.IOUNITS'
 412       include 'COMMON.CLUSTER'
 413       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 414       real*4 b
 415       dimension b(ncol)
 416       write (iout,'(a)') 'The distance matrix'
 417       do 1 i=1,ncon,ncol
 418       nlim=min0(i+ncol-1,ncon)
 419       write (iout,1000) (k,k=i,nlim)
 420       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 421  1000 format (/8x,10(i4,3x))
 422  1020 format (/1x,80(1h-)/)
 423       do 2 j=i,ncon
 424       jlim=min0(j,nlim)
 425       if (jlim.eq.j) then
 426         b(jlim-i+1)=0.0d0
 427         jlim1=jlim-1
 428       else
 429         jlim1=jlim
 430       endif
 431       do 3 k=i,jlim1
 432        if (j.lt.k) then
 433           IND=IOFFSET(NCON,j,k)
 434        else
 435           IND=IOFFSET(NCON,k,j)
 436        endif
 437     3  b(k-i+1)=diss(IND)
 438       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 439     2 continue
 440     1 continue
 441  1010 format (i5,3x,10(f6.2,1x))
 442       return
 443       end