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