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