source/cluster/wham/src-NEWSC/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
  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
  63       call initialize
  64       call openunits
  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 #ifdef SCALREP
  78       write (iout,*) "1,4 SCSC repulsive interactions sacled down by 10"
  79 #endif
  80       DO I=1,NCUT
  81         PRINTANG(I)=.FALSE.
  82         PRINTPDB(I)=0
  83         printmol2(i)=0
  84         IF (RCUTOFF(I).LT.0.0) THEN
  85           RCUTOFF(I)=ABS(RCUTOFF(I))
  86           PRINTANG(I)=.TRUE.
  87           PRINTPDB(I)=outpdb
  88           printmol2(i)=outmol2
  89         ENDIF
  90       ENDDO
  91       write (iout,*) 'Number of cutoffs:',NCUT
  92       write (iout,*) 'Cutoff values:'
  93       DO ICUT=1,NCUT
  94         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
  95      &    printpdb(icut),printmol2(icut)
  96       ENDDO
  97       DO I=1,NRES-3
  98         MULT(I)=1
  99       ENDDO
 100       do i=1,maxconf
 101         list_conf(i)=i
 102       enddo
 103       call read_coords(ncon,*20)
 104       write (iout,*) 'from read_coords: ncon',ncon
 105
 106       write (iout,*) "nT",nT
 107       do iT=1,nT
 108       write (iout,*) "iT",iT
 109 #ifdef MPI
 110       call work_partition(.true.,ncon)
 111 #endif
 112
 113       call probabl(iT,ncon_work,ncon,*20)
 114
 115       if (ncon_work.lt.2) then
 116         write (iout,*) "Too few conformations; clustering skipped"
 117         exit
 118       endif
 119 #ifdef MPI
 120       ndis=ncon_work*(ncon_work-1)/2
 121       call work_partition(.true.,ndis)
 122 #endif
 123
 124       DO I=1,NCON_work
 125         ICC(I)=I
 126       ENDDO
 127       WRITE (iout,'(A80)') TITEL
 128       t1=tcpu()
 129 C
 130 C CALCULATE DISTANCES
 131 C
 132       call daread_ccoords(1,ncon_work)
 133       ind1=0
 134       DO I=1,NCON_work-1
 135         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 136         do k=1,2*nres
 137           do l=1,3
 138             c(l,k)=allcart(l,k,i)
 139           enddo
 140         enddo
 141         do k=1,nres
 142           do l=1,3
 143             cref(l,k)=c(l,k)
 144           enddo
 145         enddo
 146         DO J=I+1,NCON_work
 147           IND=IOFFSET(NCON_work,I,J)
 148 #ifdef MPI
 149           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
 150 #endif
 151           ind1=ind1+1
 152           DISS(IND1)=DIFCONF(I,J)
 153 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
 154 #ifdef MPI
 155           endif
 156 #endif
 157         ENDDO
 158       ENDDO
 159       t2=tcpu()
 160       WRITE (iout,'(/a,1pe14.5,a/)')
 161      & 'Time for distance calculation:',T2-T1,' sec.'
 162       t1=tcpu()
 163       PRINT '(a)','End of distance computation'
 164
 165 #ifdef MPI
 166       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
 167      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
 168       if (me.eq.master) then
 169 #endif
 170       open(80,file='/tmp/distance',form='unformatted')
 171       do i=1,ndis
 172         write(80) diss(i)
 173       enddo
 174       if (punch_dist) then
 175         do i=1,ncon_work-1
 176           do j=i+1,ncon_work
 177             IND=IOFFSET(NCON,I,J)
 178             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
 179      &        energy(j)-energy(i)
 180           enddo
 181         enddo
 182       endif
 183 C
 184 C Print out the RMS deviation matrix.
 185 C
 186       if (print_dist) CALL DISTOUT(NCON_work)
 187 C
 188 C  call hierarchical clustering HC from F. Murtagh
 189 C
 190       N=NCON_work
 191       LEN = (N*(N-1))/2
 192       write(iout,*) "-------------------------------------------"
 193       write(iout,*) "HIERARCHICAL CLUSTERING using"
 194       if (iopt.eq.1) then
 195         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 196       elseif (iopt.eq.2) then
 197         write(iout,*) "SINGLE LINK METHOD"
 198       elseif (iopt.eq.3) then
 199         write(iout,*) "COMPLETE LINK METHOD"
 200       elseif (iopt.eq.4) then
 201         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 202       elseif (iopt.eq.5) then
 203         write(iout,*) "MCQUITTY'S METHOD"
 204       elseif (iopt.eq.6) then
 205         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 206       elseif (iopt.eq.7) then
 207         write(iout,*) "CENTROID METHOD"
 208       else
 209         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 210         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 211         stop
 212       endif
 213       write(iout,*)
 214       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 215       write(iout,*) "February 1986"
 216       write(iout,*) "References:"
 217       write(iout,*) "1. Multidimensional clustering algorithms"
 218       write(iout,*) "   Fionn Murtagh"
 219       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 220       write(iout,*) "2. Multivariate data analysis"
 221       write(iout,*) "   Fionn Murtagh and Andre Heck"
 222       write(iout,*) "   Kluwer Academic Publishers, 1987"
 223       write(iout,*) "-------------------------------------------"
 224       write(iout,*)
 225
 226 #ifdef DEBUG
 227       write (iout,*) "The TOTFREE array"
 228       do i=1,ncon_work
 229         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
 230       enddo
 231 #endif
 232       call flush(iout)
 233       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 234       LEV = N-1
 235       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
 236       if (lev.lt.2) then
 237         write (iout,*) "Too few conformations to cluster."
 238         goto 192
 239       endif
 240       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 241 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 242
 243       do i=1,maxgr
 244         licz(i)=0
 245       enddo
 246       icut=1
 247       i=1
 248       NGR=i+1
 249       do j=1,n
 250         licz(iclass(j,i))=licz(iclass(j,i))+1
 251         nconf(iclass(j,i),licz(iclass(j,i)))=j
 252 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 253 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 254       enddo
 255       do i=1,lev-1
 256
 257          idum=lev-i
 258          DO L=1,LEV
 259             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 260          ENDDO
 261  190     IDUM=L
 262          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
 263      &    " icut",icut," cutoff",rcutoff(icut)
 264          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 265           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 266           write (iout,'(a,f8.2)') 'Maximum distance found:',
 267      &              CRITVAL(IDUM)
 268           CALL SRTCLUST(ICUT,ncon_work,iT)
 269           CALL TRACK(ICUT)
 270           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
 271           icut=icut+1
 272           if (icut.gt.ncut) goto 191
 273          ENDIF
 274          NGR=i+1
 275          do l=1,maxgr
 276           licz(l)=0
 277          enddo
 278          do j=1,n
 279           licz(iclass(j,i))=licz(iclass(j,i))+1
 280           nconf(iclass(j,i),licz(iclass(j,i)))=j
 281 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
 282 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
 283 cd          print *,j,iclass(j,i),
 284 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 285          enddo
 286       enddo
 287  191  continue
 288 C
 289       if (plot_tree) then
 290         CALL WRITRACK
 291         CALL PLOTREE
 292       endif
 293 C
 294       t2=tcpu()
 295       WRITE (iout,'(/a,1pe14.5,a/)')
 296      & 'Total time for clustering:',T2-T1,' sec.'
 297
 298 #ifdef MPI
 299       endif
 300 #endif
 301  192  continue
 302       enddo
 303 C
 304       close(icbase,status="delete")
 305 #ifdef MPI
 306       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 307 #endif
 308       stop '********** Program terminated normally.'
 309    20 write (iout,*) "Error reading coordinates"
 310 #ifdef MPI
 311       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 312 #endif
 313       stop
 314    30 write (iout,*) "Error reading reference structure"
 315 #ifdef MPI
 316       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 317 #endif
 318       stop
 319       end
 320 c---------------------------------------------------------------------------
 321       double precision function difconf(icon,jcon)
 322       implicit none
 323       include 'DIMENSIONS'
 324       include 'sizesclu.dat'
 325       include 'COMMON.CONTROL'
 326       include 'COMMON.CLUSTER'
 327       include 'COMMON.CHAIN'
 328       include 'COMMON.INTERACT'
 329       include 'COMMON.VAR'
 330       include 'COMMON.IOUNITS'
 331       logical non_conv
 332       double precision przes(3),obrot(3,3)
 333       double precision xx(3,maxres2),yy(3,maxres2)
 334       integer i,ii,j,icon,jcon
 335       double precision rms
 336       if (lside) then
 337         ii=0
 338         do i=nstart,nend
 339           ii=ii+1
 340           do j=1,3
 341             xx(j,ii)=allcart(j,i,jcon)
 342             yy(j,ii)=cref(j,i)
 343           enddo
 344         enddo
 345         do i=nstart,nend
 346 c          if (itype(i).ne.10) then
 347             ii=ii+1
 348             do j=1,3
 349               xx(j,ii)=allcart(j,i+nres,jcon)
 350               yy(j,ii)=cref(j,i+nres)
 351             enddo
 352 c          endif
 353         enddo
 354         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 355       else
 356         do i=nnt,nct
 357           do j=1,3
 358             c(j,i)=allcart(j,i,jcon)
 359           enddo
 360         enddo
 361         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
 362      &       obrot,non_conv)
 363       endif
 364       if (rms.lt.0.0) then
 365         print *,'error, rms^2 = ',rms,icon,jcon
 366         stop
 367       endif
 368       if (non_conv) print *,non_conv,icon,jcon
 369       difconf=dsqrt(rms)
 370       RETURN
 371       END
 372 C------------------------------------------------------------------------------
 373       subroutine distout(ncon)
 374       implicit none
 375       include 'DIMENSIONS'
 376       include 'sizesclu.dat'
 377       integer ncol,ncon
 378       parameter (ncol=10)
 379       include 'COMMON.IOUNITS'
 380       include 'COMMON.CLUSTER'
 381       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
 382       real*4 b
 383       dimension b(ncol)
 384       write (iout,'(a)') 'The distance matrix'
 385       do 1 i=1,ncon,ncol
 386       nlim=min0(i+ncol-1,ncon)
 387       write (iout,1000) (k,k=i,nlim)
 388       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 389  1000 format (/8x,10(i4,3x))
 390  1020 format (/1x,80(1h-)/)
 391       do 2 j=i,ncon
 392       jlim=min0(j,nlim)
 393       if (jlim.eq.j) then
 394         b(jlim-i+1)=0.0d0
 395         jlim1=jlim-1
 396       else
 397         jlim1=jlim
 398       endif
 399       do 3 k=i,jlim1
 400        if (j.lt.k) then
 401           IND=IOFFSET(NCON,j,k)
 402        else
 403           IND=IOFFSET(NCON,k,j)
 404        endif
 405     3  b(k-i+1)=diss(IND)
 406       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 407     2 continue
 408     1 continue
 409  1010 format (i5,3x,10(f6.2,1x))
 410       return
 411       end