source/cluster/unres/src/main_clust.f

   1 C
   2 C Program to cluster united-residue MCM results.
   3 C
   4       include 'DIMENSIONS'
   5       include 'sizesclu.dat'
   6       include 'COMMON.TIME1'
   7       include 'COMMON.INTERACT'
   8       include 'COMMON.NAMES'
   9       include 'COMMON.GEO'
  10       include 'COMMON.HEADER'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CONTACTS'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.VAR'
  15       include 'COMMON.CLUSTER'
  16       include 'COMMON.IOUNITS'
  17       logical printang(max_cut)
  18       integer printpdb(max_cut)
  19       integer printmol2(max_cut)
  20       character*240 lineh
  21       REAL CRIT(maxconf),MEMBR(maxconf)
  22       REAL CRITVAL(maxconf-1)
  23       REAL DISS(maxdist)
  24       INTEGER IA(maxconf),IB(maxconf)
  25       INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
  26       INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
  27       DIMENSION NN(maxconf),DISNN(maxconf)
  28       LOGICAL FLAG(maxconf)
  29
  30       double precision varia(maxvar)
  31       double precision hrtime,mintime,sectime
  32       logical eof
  33       double precision eee(n_ene,maxconf)
  34
  35       call initialize
  36       call openunits
  37       call read_control
  38       call molread
  39       do i=1,nres
  40         phi(i)=0.0D0
  41         theta(i)=0.0D0
  42         alph(i)=0.0D0
  43         omeg(i)=0.0D0
  44       enddo
  45
  46       print *,'MAIN: nnt=',nnt,' nct=',nct
  47
  48       DO I=1,NCUT
  49         PRINTANG(I)=.FALSE.
  50         PRINTPDB(I)=0
  51         printmol2(i)=0
  52         IF (RCUTOFF(I).LT.0.0) THEN
  53           RCUTOFF(I)=ABS(RCUTOFF(I))
  54           PRINTANG(I)=.TRUE.
  55           PRINTPDB(I)=outpdb
  56           printmol2(i)=outmol2
  57         ENDIF
  58       ENDDO
  59       PRINT *,'Number of cutoffs:',NCUT
  60       PRINT *,'Cutoff values:'
  61       DO ICUT=1,NCUT
  62         PRINT '(8HRCUTOFF(,I2,2H)=,F8.1,2i2)',ICUT,RCUTOFF(ICUT),
  63      &    printpdb(icut),printmol2(icut)
  64       ENDDO
  65       DO I=1,NRES-3
  66         MULT(I)=1
  67       ENDDO
  68       ICON=1
  69   123 continue
  70       if (from_cart) then
  71         if (efree) then
  72         read (intin,*,end=13,err=11) energy(icon),totfree(icon),
  73      &    rmstab(icon),
  74      &    nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
  75      &    i=1,nss_all(icon)),iscore(icon)
  76         else
  77         read (intin,*,end=13,err=11) energy(icon),rmstab(icon),
  78      &    nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
  79      &    i=1,nss_all(icon)),iscore(icon)
  80         endif
  81         read (intin,'(8f10.5)',end=13,err=10)
  82      &    ((allcart(j,i,icon),j=1,3),i=1,nres),
  83      &    ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
  84         print *,icon,energy(icon),nss_all(icon),rmstab(icon)
  85       else
  86         read(intin,'(a80)',end=13,err=12) lineh
  87         read(lineh(:5),*,err=8) ic
  88         if (efree) then
  89         read(lineh(6:),*,err=8) energy(icon)
  90         else
  91         read(lineh(6:),*,err=8) energy(icon)
  92         endif
  93         goto 9
  94     8   ic=1
  95         print *,'error, assuming e=1d10',lineh
  96         energy(icon)=1d10
  97         nss=0
  98     9   continue
  99 cold        read(lineh(18:),*,end=13,err=11) nss_all(icon)
 100         ii = index(lineh(15:)," ")+15
 101         read(lineh(ii:),*,end=13,err=11) nss_all(icon)
 102         IF (NSS_all(icon).LT.9) THEN
 103           read (lineh(20:),*,end=102)
 104      &    (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),
 105      &    iscore(icon)
 106         ELSE
 107           read (lineh(20:),*,end=102)
 108      &           (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
 109           read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),
 110      &      I=9,NSS_all(icon)),iscore(icon)
 111         ENDIF
 112
 113   102   continue
 114
 115         PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
 116         totfree(icon)=energy(icon)
 117         call read_angles(intin,*13)
 118         call chainbuild
 119         do ii=1,2*nres
 120           do jj=1,3
 121             allcart(jj,ii,icon)=c(jj,ii)
 122           enddo
 123         enddo
 124         do i=1,nres
 125           phiall(i,icon)=phi(i)
 126           thetall(i,icon)=theta(i)
 127           alphall(i,icon)=alph(i)
 128           omall(i,icon)=omeg(i)
 129         enddo
 130       endif
 131       ICON=ICON+1
 132       GOTO 123
 133 C
 134 C CALCULATE DISTANCES
 135 C
 136    10 print *,'something wrong with angles'
 137       goto 13
 138    11 print *,'something wrong with NSS',nss
 139       goto 13
 140    12 print *,'something wrong with header'
 141
 142    13 NCON=ICON-1
 143       if (lgrp) then
 144         i=0
 145         do while (i.lt.ncon)
 146 c          read (jstatin,*,err=1111,end=1111) idum,(eee(j,i),j=1,n_ene)
 147           read (jstatin,'(a)') lineh
 148           if (index(lineh,'#').gt.0) goto 1123
 149           i=i+1
 150           read (lineh,*) idum,(eee(j,i),j=1,n_ene)
 151           print *,i,idum,(eee(j,i),j=1,n_ene)
 152           energy(i)=eee(15,i)
 153  1123     continue
 154         enddo
 155       endif
 156       goto 1112
 157  1111 print *,'Error in the statin file; statout will not be produced.'
 158       write (iout,*)
 159      &  'Error in the statin file; statout will not be produced.'
 160       lgrp=.false.
 161  1112 continue
 162       print *,'ncon',ncon
 163       DO I=1,NCON
 164         ICC(I)=I
 165       ENDDO
 166       WRITE (iout,'(A80)') TITEL
 167       t1=tcpu()
 168       DO I=1,NCON-1
 169         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
 170         if (from_cart) then
 171           do ii=1,2*nres
 172             do jj=1,3
 173               c(jj,ii)=allcart(jj,ii,i)
 174             enddo
 175           enddo
 176           call int_from_cart(.true.,.false.)
 177           do ii=1,nres
 178             phiall(ii,i)=phi(ii)
 179             thetall(ii,i)=theta(ii)
 180             alphall(ii,i)=alph(ii)
 181             omall(ii,i)=omeg(ii)
 182           enddo
 183         else
 184           do ii=1,nres
 185             phi(ii)=phiall(ii,i)
 186             theta(ii)=thetall(ii,i)
 187             alph(ii)=alphall(ii,i)
 188             omeg(ii)=omall(ii,i)
 189           enddo
 190           call chainbuild
 191         endif
 192         do ii=1,nres
 193           do jj=1,3
 194             cref(jj,ii)=c(jj,ii)
 195           enddo
 196         enddo
 197         DO J=I+1,NCON
 198           IND=IOFFSET(NCON,I,J)
 199           ATTALUMS(IND)=DIFCONF(I,J)
 200 c          write (iout,'(2i4,f10.5)') i,j,ATTALUMS(IND)
 201           DISS(IND)=ATTALUMS(IND)
 202         ENDDO
 203       ENDDO
 204       t2=tcpu()
 205       WRITE (iout,'(/a,1pe14.5,a/)')
 206      & 'Time for distance calculation:',T2-T1,' sec.'
 207       t1=tcpu()
 208       PRINT '(a)','End of distance computation'
 209
 210       if (punch_dist) then
 211         do i=1,ncon-1
 212           do j=i+1,ncon
 213             IND=IOFFSET(NCON,I,J)
 214             write (jrms,'(2i5,2f10.5)') i,j,attalums(IND),
 215      &        energy(j)-energy(i)
 216           enddo
 217         enddo
 218       endif
 219 C
 220 C Print out the RMS deviation matrix.
 221 C
 222       if (print_dist) CALL DISTOUT(NCON)
 223 C
 224 C  call hierarchical clustering HC from F. Murtagh
 225 C
 226       N=NCON
 227       LEN = (N*(N-1))/2
 228       write(iout,*) "-------------------------------------------"
 229       write(iout,*) "HIERARCHICAL CLUSTERING using"
 230       if (iopt.eq.1) then
 231         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
 232       elseif (iopt.eq.2) then
 233         write(iout,*) "SINGLE LINK METHOD"
 234       elseif (iopt.eq.3) then
 235         write(iout,*) "COMPLETE LINK METHOD"
 236       elseif (iopt.eq.4) then
 237         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
 238       elseif (iopt.eq.5) then
 239         write(iout,*) "MCQUITTY'S METHOD"
 240       elseif (iopt.eq.6) then
 241         write(iout,*) "MEDIAN (GOWER'S) METHOD"
 242       elseif (iopt.eq.7) then
 243         write(iout,*) "CENTROID METHOD"
 244       else
 245         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
 246         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
 247         stop
 248       endif
 249       write(iout,*)
 250       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
 251       write(iout,*) "February 1986"
 252       write(iout,*) "References:"
 253       write(iout,*) "1. Multidimensional clustering algorithms"
 254       write(iout,*) "   Fionn Murtagh"
 255       write(iout,*) "   Vienna : Physica-Verlag, 1985."
 256       write(iout,*) "2. Multivariate data analysis"
 257       write(iout,*) "   Fionn Murtagh and Andre Heck"
 258       write(iout,*) "   Kluwer Academic Publishers, 1987"
 259       write(iout,*) "-------------------------------------------"
 260       write(iout,*)
 261
 262       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
 263       LEV = N-1
 264       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 265       CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 266
 267       icut=1
 268       i=1
 269       NGR=i+1
 270       do j=1,n
 271         licz(iclass(j,i))=licz(iclass(j,i))+1
 272         nconf(iclass(j,i),licz(iclass(j,i)))=j
 273       enddo
 274       do i=1,lev-1
 275
 276          idum=lev-i
 277          DO L=1,LEV
 278             IF (HEIGHT(L).EQ.IDUM) GOTO 190
 279          ENDDO
 280  190     IDUM=L
 281 cd         print *,i+1,CRITVAL(IDUM)
 282          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
 283           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
 284           write (iout,'(a,f8.2)') 'Maximum distance found:',
 285      &              CRITVAL(IDUM)
 286           CALL SRTCLUST(ICUT,ncon)
 287           CALL TRACK(ICUT)
 288           CALL WRTCLUST(ncon,ICUT,PRINTANG,PRINTPDB,PRINTMOL2)
 289           icut=icut+1
 290           if (icut.gt.ncut) goto 191
 291          ENDIF
 292          NGR=i+1
 293          do l=1,maxgr
 294           licz(l)=0
 295          enddo
 296          do j=1,n
 297           licz(iclass(j,i))=licz(iclass(j,i))+1
 298           nconf(iclass(j,i),licz(iclass(j,i)))=j
 299 cd          print *,j,iclass(j,i),
 300 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
 301          enddo
 302       enddo
 303  191  continue
 304 C
 305       if (plot_tree) then
 306         CALL WRITRACK
 307         CALL PLOTREE
 308       endif
 309 C
 310       t2=tcpu()
 311       WRITE (iout,'(/a,1pe14.5,a/)')
 312      & 'Total time for clustering:',T2-T1,' sec.'
 313       if (lgrp) then
 314         write (jstatout,'(a5,18a12,10f4.1)')"#    ",
 315      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p",
 316      &   "EBE bending","ESC SCloc","ETORS ",
 317      &   "ECORR4 ","ECORR5 ","ECORR6 ",
 318      &   "EELLO ","ETURN3 ","ETURN4 ","ETURN6 "
 319      &   ,"ETOT total","RMSD","nat.contact","nnt.contact",
 320      &   "cont.order",(rcutoff(i),i=1,ncut)
 321         do i=1,ncon
 322           write(jstatout,'(i5,18f12.4,10i6)') i,(eee(j,i),j=1,n_ene),
 323      &     (iass_tot(i,icut),icut=1,ncut)
 324         enddo
 325       endif
 326 C
 327       stop '********** Program terminated normally.'
 328       end
 329 c---------------------------------------------------------------------------
 330       double precision function difconf(icon,jcon)
 331       include 'DIMENSIONS'
 332       include 'sizesclu.dat'
 333       include 'COMMON.CONTROL'
 334       include 'COMMON.CLUSTER'
 335       include 'COMMON.CHAIN'
 336       include 'COMMON.INTERACT'
 337       include 'COMMON.VAR'
 338       logical non_conv
 339       double precision przes(3),obrot(3,3)
 340       double precision xx(3,maxres2),yy(3,maxres2)
 341       do i=1,nres
 342         phi(i)=phiall(i,jcon)
 343         theta(i)=thetall(i,jcon)
 344         alph(i)=alphall(i,jcon)
 345         omeg(i)=omall(i,jcon)
 346       enddo
 347       call chainbuild
 348 c     do i=1,nres
 349 c       print '(i4,2(3f10.5,5x))',i,(cref(j,i),j=1,3),(c(j,i),j=1,3)
 350 c     enddo
 351       if (lside) then
 352         ii=0
 353         do i=nnt,nct
 354           ii=ii+1
 355           do j=1,3
 356             xx(j,ii)=allcart(j,i,jcon)
 357             yy(j,ii)=cref(j,i)
 358           enddo
 359         enddo
 360         do i=nnt,nct
 361 c          if (itype(i).ne.10) then
 362             ii=ii+1
 363             do j=1,3
 364               xx(j,ii)=allcart(j,i+nres,jcon)
 365               yy(j,ii)=cref(j,i+nres)
 366             enddo
 367 c          endif
 368         enddo
 369         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 370       else
 371         do i=nnt,nct
 372           do j=1,3
 373             c(j,i)=allcart(j,i,jcon)
 374           enddo
 375 c          write (2,'(i5,3f10.5,5x,3f10.5)') i,
 376 c     &     (c(j,i),j=1,3),(cref(j,i),j=1,3)
 377         enddo
 378         call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obrot
 379      &       ,non_conv)
 380       endif
 381       if (rms.lt.0.0) then
 382         print *,error,'rms^2 = ',rms,icon,jcon
 383         stop
 384       endif
 385       if (non_conv) print *,non_conv,icon,jcon
 386       difconf=dsqrt(rms)
 387       RETURN
 388       END
 389 C------------------------------------------------------------------------------
 390       double precision function rmsnat(jcon)
 391       include 'DIMENSIONS'
 392       include 'sizesclu.dat'
 393       include 'COMMON.CONTROL'
 394       include 'COMMON.CLUSTER'
 395       include 'COMMON.CHAIN'
 396       include 'COMMON.INTERACT'
 397       include 'COMMON.VAR'
 398       logical non_conv
 399       double precision przes(3),obrot(3,3)
 400       double precision xx(3,maxres2),yy(3,maxres2)
 401       if (lside) then
 402         ii=0
 403         do i=nnt,nct
 404           ii=ii+1
 405           do j=1,3
 406             xx(j,ii)=allcart(j,i,jcon)
 407             yy(j,ii)=cref_pdb(j,i)
 408           enddo
 409         enddo
 410         do i=nnt,nct
 411 c          if (itype(i).ne.10) then
 412             ii=ii+1
 413             do j=1,3
 414               xx(j,ii)=allcart(j,i+nres,jcon)
 415               yy(j,ii)=cref_pdb(j,i+nres)
 416             enddo
 417 c          endif
 418         enddo
 419         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
 420       else
 421         do i=nnt,nct
 422           do j=1,3
 423             c(j,i)=allcart(j,i,jcon)
 424           enddo
 425         enddo
 426         call fitsq(rms,c(1,nnt),cref_pdb(1,nnt),nct-nnt+1,przes,obrot
 427      &       ,non_conv)
 428       endif
 429       if (rms.lt.0.0) then
 430         print *,error,'rms^2 = ',rms,icon,jcon
 431         stop
 432       endif
 433       if (non_conv) print *,non_conv,icon,jcon
 434       rmsnat=dsqrt(rms)
 435       RETURN
 436       END
 437 C------------------------------------------------------------------------------
 438       subroutine distout(ncon)
 439       include 'DIMENSIONS'
 440       include 'sizesclu.dat'
 441       parameter (ncol=10)
 442       include 'COMMON.IOUNITS'
 443       include 'COMMON.CLUSTER'
 444       dimension b(ncol)
 445       write (iout,'(a)') 'The distance matrix'
 446       do 1 i=1,ncon,ncol
 447       nlim=min0(i+ncol-1,ncon)
 448       write (iout,1000) (k,k=i,nlim)
 449       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 450  1000 format (/8x,10(i4,3x))
 451  1020 format (/1x,80(1h-)/)
 452       do 2 j=i,ncon
 453       jlim=min0(j,nlim)
 454       if (jlim.eq.j) then
 455         b(jlim-i+1)=0.0d0
 456         jlim1=jlim-1
 457       else
 458         jlim1=jlim
 459       endif
 460       do 3 k=i,jlim1
 461        if (j.lt.k) then
 462           IND=IOFFSET(NCON,j,k)
 463        else
 464           IND=IOFFSET(NCON,k,j)
 465        endif
 466     3  b(k-i+1)=attalums(IND)
 467       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 468     2 continue
 469     1 continue
 470  1010 format (i5,3x,10(f6.2,1x))
 471       return
 472       end
 473 C-----------------------------------------------------------------------
 474       block data
 475       include 'DIMENSIONS'
 476       include 'COMMON.LOCAL'
 477       data dsc /  1.237,  ! CYS (type  1)
 478      &            2.142,  ! MET (type  2)
 479      &            2.299,  ! PHE (type  3)
 480      &            1.776,  ! ILE (type  4)
 481      &            1.939,  ! LEU (type  5)
 482      &            1.410,  ! VAL (type  6)
 483      &            2.605,  ! TRP (type  7)
 484      &            2.484,  ! TYR (type  8)
 485      &            0.743,  ! ALA (type  9)
 486      &            0.000,  ! GLY (type 10)
 487      &            1.393,  ! THR (type 11)
 488      &            1.150,  ! SER (type 12)
 489      &            2.240,  ! GLN (type 13)
 490      &            1.684,  ! ASN (type 14)
 491      &            2.254,  ! GLU (type 15)
 492      &            1.709,  ! ASP (type 16)
 493      &            2.113,  ! HIS (type 17)
 494      &            3.020,  ! ARG (type 18)
 495      &            2.541,  ! LYS (type 19)
 496      &            1.345,  ! PRO (type 20)
 497      &            0.000  /! D   (type 21)
 498       end