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
  64       call initialize
  65       call openunits
  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       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
 338       double precision rms
 339       if (lside) then
 340         ii=0
 341         do i=nstart,nend
 342           ii=ii+1
 343           do j=1,3
 344             xx(j,ii)=allcart(j,i,jcon)
 345             yy(j,ii)=cref(j,i)
 346           enddo
 347         enddo
 348         do i=nstart,nend
 349 c          if (itype(i).ne.10) then
 350             ii=ii+1
 351             do j=1,3
 352               xx(j,ii)=allcart(j,i+nres,jcon)
 353               yy(j,ii)=cref(j,i+nres)
 354             enddo
 355 c          endif
 356         enddo
 357         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 358       else
 359         do i=nnt,nct
 360           do j=1,3
 361             c(j,i)=allcart(j,i,jcon)
 362           enddo
 363         enddo
 364         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
 365      &       obrot,non_conv)
 366       endif
 367       if (rms.lt.0.0) then
 368         print *,'error, rms^2 = ',rms,icon,jcon
 369         stop
 370       endif
 371       if (non_conv) print *,non_conv,icon,jcon
 372       difconf=dsqrt(rms)
 373       RETURN
 374       END
 375 C------------------------------------------------------------------------------
 376       subroutine distout(ncon)
 377       implicit none
 378       include 'DIMENSIONS'
 379       include 'sizesclu.dat'
 380       integer ncol,ncon
 381       parameter (ncol=10)
 382       include 'COMMON.IOUNITS'
 383       include 'COMMON.CLUSTER'
 384       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 385       real*4 b
 386       dimension b(ncol)
 387       write (iout,'(a)') 'The distance matrix'
 388       do 1 i=1,ncon,ncol
 389       nlim=min0(i+ncol-1,ncon)
 390       write (iout,1000) (k,k=i,nlim)
 391       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 392  1000 format (/8x,10(i4,3x))
 393  1020 format (/1x,80(1h-)/)
 394       do 2 j=i,ncon
 395       jlim=min0(j,nlim)
 396       if (jlim.eq.j) then
 397         b(jlim-i+1)=0.0d0
 398         jlim1=jlim-1
 399       else
 400         jlim1=jlim
 401       endif
 402       do 3 k=i,jlim1
 403        if (j.lt.k) then
 404           IND=IOFFSET(NCON,j,k)
 405        else
 406           IND=IOFFSET(NCON,k,j)
 407        endif
 408     3  b(k-i+1)=diss(IND)
 409       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 410     2 continue
 411     1 continue
 412  1010 format (i5,3x,10(f6.2,1x))
 413       return
 414       end