source/cluster/wham/src-new/main_clust.F.org

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