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