+C
+C Program to cluster united-residue MCM results.
+C
+ include 'DIMENSIONS'
+ include 'sizesclu.dat'
+ include 'COMMON.TIME1'
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.GEO'
+ include 'COMMON.HEADER'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ include 'COMMON.CLUSTER'
+ include 'COMMON.IOUNITS'
+ logical printang(max_cut)
+ integer printpdb(max_cut)
+ integer printmol2(max_cut)
+ character*240 lineh
+ REAL CRIT(maxconf),MEMBR(maxconf)
+ REAL CRITVAL(maxconf-1)
+ REAL DISS(maxdist)
+ INTEGER IA(maxconf),IB(maxconf)
+ INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
+ INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
+ DIMENSION NN(maxconf),DISNN(maxconf)
+ LOGICAL FLAG(maxconf)
+ integer len
+
+ double precision varia(maxvar)
+ double precision hrtime,mintime,sectime
+ logical eof
+ double precision eee(n_ene,maxconf)
+
+ call initialize
+ call openunits
+ call read_control
+ call molread
+ do i=1,nres
+ phi(i)=0.0D0
+ theta(i)=0.0D0
+ alph(i)=0.0D0
+ omeg(i)=0.0D0
+ enddo
+
+ print *,'MAIN: nnt=',nnt,' nct=',nct
+
+ DO I=1,NCUT
+ PRINTANG(I)=.FALSE.
+ PRINTPDB(I)=0
+ printmol2(i)=0
+ IF (RCUTOFF(I).LT.0.0) THEN
+ RCUTOFF(I)=ABS(RCUTOFF(I))
+ PRINTANG(I)=.TRUE.
+ PRINTPDB(I)=outpdb
+ printmol2(i)=outmol2
+ ENDIF
+ ENDDO
+ PRINT *,'Number of cutoffs:',NCUT
+ PRINT *,'Cutoff values:'
+ DO ICUT=1,NCUT
+ PRINT '(8HRCUTOFF(,I2,2H)=,F8.1,2i2)',ICUT,RCUTOFF(ICUT),
+ & printpdb(icut),printmol2(icut)
+ ENDDO
+ DO I=1,NRES-3
+ MULT(I)=1
+ ENDDO
+ if (from_cx) then
+
+ call cxread(ncon,*101)
+ goto 111
+ 101 write (iout,*) "Error in CX file"
+ stop
+ 111 continue
+
+ else
+
+ ICON=1
+ 123 continue
+ if (from_cart) then
+ if (efree) then
+ read (intin,*,end=13,err=11) energy(icon),totfree(icon),
+ & rmstab(icon),
+ & nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
+ & i=1,nss_all(icon)),iscore(icon)
+ else
+ read (intin,*,end=13,err=11) energy(icon),rmstab(icon),
+ & nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
+ & i=1,nss_all(icon)),iscore(icon)
+ endif
+ read (intin,'(8f10.5)',end=13,err=10)
+ & ((allcart(j,i,icon),j=1,3),i=1,nres),
+ & ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
+ print *,icon,energy(icon),nss_all(icon),rmstab(icon)
+ else
+ read(intin,'(a80)',end=13,err=12) lineh
+ read(lineh(:5),*,err=8) ic
+ if (efree) then
+ read(lineh(6:),*,err=8) energy(icon)
+ else
+ read(lineh(6:),*,err=8) energy(icon)
+ endif
+ goto 9
+ 8 ic=1
+ print *,'error, assuming e=1d10',lineh
+ energy(icon)=1d10
+ nss=0
+ 9 continue
+cold read(lineh(18:),*,end=13,err=11) nss_all(icon)
+ ii = index(lineh(15:)," ")+15
+ read(lineh(ii:),*,end=13,err=11) nss_all(icon)
+ IF (NSS_all(icon).LT.9) THEN
+ read (lineh(20:),*,end=102)
+ & (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),
+ & iscore(icon)
+ ELSE
+ read (lineh(20:),*,end=102)
+ & (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
+ read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),
+ & I=9,NSS_all(icon)),iscore(icon)
+ ENDIF
+
+ 102 continue
+
+ PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
+ totfree(icon)=energy(icon)
+ call read_angles(intin,*13)
+ call chainbuild
+ do ii=1,2*nres
+ do jj=1,3
+ allcart(jj,ii,icon)=c(jj,ii)
+ enddo
+ enddo
+ do i=1,nres
+ phiall(i,icon)=phi(i)
+ thetall(i,icon)=theta(i)
+ alphall(i,icon)=alph(i)
+ omall(i,icon)=omeg(i)
+ enddo
+ endif
+ ICON=ICON+1
+ GOTO 123
+ 10 print *,'something wrong with angles'
+ goto 13
+ 11 print *,'something wrong with NSS',nss
+ goto 13
+ 12 print *,'something wrong with header'
+
+ 13 NCON=ICON-1
+
+ endif
+
+#ifdef DEBUG
+ do icon=1,ncon
+ write (iout,*) "Conformation",icon
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,
+ & (allcart(j,i,icon),j=1,3),(allcart(j,i+nres,icon),j=1,3)
+ enddo
+ enddo
+#endif
+ if (lgrp) then
+ i=0
+ do while (i.lt.ncon)
+c read (jstatin,*,err=1111,end=1111) idum,(eee(j,i),j=1,n_ene)
+ read (jstatin,'(a)') lineh
+ if (index(lineh,'#').gt.0) goto 1123
+ i=i+1
+ read (lineh,*) idum,(eee(j,i),j=1,n_ene)
+ print *,i,idum,(eee(j,i),j=1,n_ene)
+ energy(i)=eee(15,i)
+ 1123 continue
+ enddo
+ endif
+ goto 1112
+ 1111 print *,'Error in the statin file; statout will not be produced.'
+ write (iout,*)
+ & 'Error in the statin file; statout will not be produced.'
+ lgrp=.false.
+ 1112 continue
+ print *,'ncon',ncon
+ DO I=1,NCON
+ ICC(I)=I
+ ENDDO
+ WRITE (iout,'(A80)') TITEL
+ t1=tcpu()
+C
+C CALCULATE DISTANCES
+C
+ DO I=1,NCON-1
+ if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
+ if (from_cart .or. from_cx) then
+ do ii=1,2*nres
+ do jj=1,3
+ c(jj,ii)=allcart(jj,ii,i)
+ enddo
+ enddo
+ call int_from_cart(.true.,.false.)
+ do ii=1,nres
+ phiall(ii,i)=phi(ii)
+ thetall(ii,i)=theta(ii)
+ alphall(ii,i)=alph(ii)
+ omall(ii,i)=omeg(ii)
+ enddo
+ else
+ do ii=1,nres
+ phi(ii)=phiall(ii,i)
+ theta(ii)=thetall(ii,i)
+ alph(ii)=alphall(ii,i)
+ omeg(ii)=omall(ii,i)
+ enddo
+ call chainbuild
+ endif
+ do ii=1,nres
+ do jj=1,3
+ cref(jj,ii)=c(jj,ii)
+ enddo
+ enddo
+ DO J=I+1,NCON
+ IND=IOFFSET(NCON,I,J)
+ ATTALUMS(IND)=DIFCONF(I,J)
+c write (iout,'(2i4,f10.5)') i,j,ATTALUMS(IND)
+ DISS(IND)=ATTALUMS(IND)
+ ENDDO
+ ENDDO
+ t2=tcpu()
+ WRITE (iout,'(/a,1pe14.5,a/)')
+ & 'Time for distance calculation:',T2-T1,' sec.'
+ t1=tcpu()
+ PRINT '(a)','End of distance computation'
+
+ if (punch_dist) then
+ do i=1,ncon-1
+ do j=i+1,ncon
+ IND=IOFFSET(NCON,I,J)
+ write (jrms,'(2i5,2f10.5)') i,j,attalums(IND),
+ & energy(j)-energy(i)
+ enddo
+ enddo
+ endif
+C
+C Print out the RMS deviation matrix.
+C
+ if (print_dist) CALL DISTOUT(NCON)
+C
+C call hierarchical clustering HC from F. Murtagh
+C
+ N=NCON
+ LEN = (N*(N-1))/2
+ write(iout,*) "-------------------------------------------"
+ write(iout,*) "HIERARCHICAL CLUSTERING using"
+ if (iopt.eq.1) then
+ write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
+ elseif (iopt.eq.2) then
+ write(iout,*) "SINGLE LINK METHOD"
+ elseif (iopt.eq.3) then
+ write(iout,*) "COMPLETE LINK METHOD"
+ elseif (iopt.eq.4) then
+ write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
+ elseif (iopt.eq.5) then
+ write(iout,*) "MCQUITTY'S METHOD"
+ elseif (iopt.eq.6) then
+ write(iout,*) "MEDIAN (GOWER'S) METHOD"
+ elseif (iopt.eq.7) then
+ write(iout,*) "CENTROID METHOD"
+ else
+ write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
+ write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
+ stop
+ endif
+ write(iout,*)
+ write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
+ write(iout,*) "February 1986"
+ write(iout,*) "References:"
+ write(iout,*) "1. Multidimensional clustering algorithms"
+ write(iout,*) " Fionn Murtagh"
+ write(iout,*) " Vienna : Physica-Verlag, 1985."
+ write(iout,*) "2. Multivariate data analysis"
+ write(iout,*) " Fionn Murtagh and Andre Heck"
+ write(iout,*) " Kluwer Academic Publishers, 1987"
+ write(iout,*) "-------------------------------------------"
+ write(iout,*)
+ CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
+c print *,"HC"
+ LEV = N-1
+ CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
+c print *,"HCASS"
+c print *,"lev",lev
+c print *,"call HCDEN"
+ CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
+c print *,"HCDEN"
+
+ icut=1
+ i=1
+ NGR=i+1
+ do j=1,n
+ licz(iclass(j,i))=licz(iclass(j,i))+1
+ nconf(iclass(j,i),licz(iclass(j,i)))=j
+ enddo
+ do i=1,lev-1
+
+ idum=lev-i
+ DO L=1,LEV
+ IF (HEIGHT(L).EQ.IDUM) GOTO 190
+ ENDDO
+ 190 IDUM=L
+cd print *,i+1,CRITVAL(IDUM)
+ IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
+ WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
+ write (iout,'(a,f8.2)') 'Maximum distance found:',
+ & CRITVAL(IDUM)
+ CALL SRTCLUST(ICUT,ncon)
+ CALL TRACK(ICUT)
+ CALL WRTCLUST(ncon,ICUT,PRINTANG,PRINTPDB,PRINTMOL2)
+ icut=icut+1
+ if (icut.gt.ncut) goto 191
+ ENDIF
+ NGR=i+1
+ do l=1,maxgr
+ licz(l)=0
+ enddo
+ do j=1,n
+ licz(iclass(j,i))=licz(iclass(j,i))+1
+ nconf(iclass(j,i),licz(iclass(j,i)))=j
+cd print *,j,iclass(j,i),
+cd & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
+ enddo
+ enddo
+ 191 continue
+C
+ if (plot_tree) then
+ CALL WRITRACK
+ CALL PLOTREE
+ endif
+C
+ t2=tcpu()
+ WRITE (iout,'(/a,1pe14.5,a/)')
+ & 'Total time for clustering:',T2-T1,' sec.'
+ if (lgrp) then
+ write (jstatout,'(a5,18a12,10f4.1)')"# ",
+ & "EVDW SC-SC","EVDW2 SC-p","EES p-p",
+ & "EBE bending","ESC SCloc","ETORS ",
+ & "ECORR4 ","ECORR5 ","ECORR6 ",
+ & "EELLO ","ETURN3 ","ETURN4 ","ETURN6 "
+ & ,"ETOT total","RMSD","nat.contact","nnt.contact",
+ & "cont.order",(rcutoff(i),i=1,ncut)
+ do i=1,ncon
+ write(jstatout,'(i5,18f12.4,10i6)') i,(eee(j,i),j=1,n_ene),
+ & (iass_tot(i,icut),icut=1,ncut)
+ enddo
+ endif
+C
+ stop '********** Program terminated normally.'
+ end
+c---------------------------------------------------------------------------
+ double precision function difconf(icon,jcon)
+ include 'DIMENSIONS'
+ include 'sizesclu.dat'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CLUSTER'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ logical non_conv
+ double precision przes(3),obrot(3,3)
+ double precision xx(3,maxres2),yy(3,maxres2)
+c do i=1,nres
+c phi(i)=phiall(i,jcon)
+c theta(i)=thetall(i,jcon)
+c alph(i)=alphall(i,jcon)
+c omeg(i)=omall(i,jcon)
+c enddo
+c call chainbuild
+c do i=1,nres
+c print '(i4,2(3f10.5,5x))',i,(cref(j,i),j=1,3),(c(j,i),j=1,3)
+c enddo
+ if (lside) then
+ ii=0
+ do i=nnt,nct
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i,jcon)
+ yy(j,ii)=allcart(j,i,icon)
+ enddo
+ enddo
+ do i=nnt,nct
+c if (itype(i).ne.10) then
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i+nres,jcon)
+ yy(j,ii)=allcart(j,i+nres,icon)
+ enddo
+c endif
+ enddo
+ call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
+ else
+ do i=nnt,nct
+c do j=1,3
+c c(j,i)=allcart(j,i,jcon)
+c enddo
+c write (2,'(i5,3f10.5,5x,3f10.5)') i,
+c & (c(j,i),j=1,3),(cref(j,i),j=1,3)
+ enddo
+ call fitsq(rms,allcart(1,nnt,jcon),allcart(1,nnt,icon),
+ & nct-nnt+1,przes,obrot,non_conv)
+ endif
+ if (rms.lt.0.0) then
+ print *,error,'rms^2 = ',rms,icon,jcon
+ stop
+ endif
+ if (non_conv) print *,non_conv,icon,jcon
+ difconf=dsqrt(rms)
+ RETURN
+ END
+C------------------------------------------------------------------------------
+ double precision function rmsnat(jcon)
+ include 'DIMENSIONS'
+ include 'sizesclu.dat'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CLUSTER'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ logical non_conv
+ double precision przes(3),obrot(3,3)
+ double precision xx(3,maxres2),yy(3,maxres2)
+ if (lside) then
+ ii=0
+ do i=nnt,nct
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i,jcon)
+ yy(j,ii)=cref_pdb(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+c if (itype(i).ne.10) then
+ ii=ii+1
+ do j=1,3
+ xx(j,ii)=allcart(j,i+nres,jcon)
+ yy(j,ii)=cref_pdb(j,i+nres)
+ enddo
+c endif
+ enddo
+ call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
+ else
+ do i=nnt,nct
+ do j=1,3
+ c(j,i)=allcart(j,i,jcon)
+ enddo
+ enddo
+ call fitsq(rms,c(1,nnt),cref_pdb(1,nnt),nct-nnt+1,przes,obrot
+ & ,non_conv)
+ endif
+ if (rms.lt.0.0) then
+ print *,error,'rms^2 = ',rms,icon,jcon
+ stop
+ endif
+ if (non_conv) print *,non_conv,icon,jcon
+ rmsnat=dsqrt(rms)
+ RETURN
+ END
+C------------------------------------------------------------------------------
+ subroutine distout(ncon)
+ include 'DIMENSIONS'
+ include 'sizesclu.dat'
+ parameter (ncol=10)
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CLUSTER'
+ dimension b(ncol)
+ write (iout,'(a)') 'The distance matrix'
+ do 1 i=1,ncon,ncol
+ nlim=min0(i+ncol-1,ncon)
+ write (iout,1000) (k,k=i,nlim)
+ write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
+ 1000 format (/8x,10(i4,3x))
+ 1020 format (/1x,80(1h-)/)
+ do 2 j=i,ncon
+ jlim=min0(j,nlim)
+ if (jlim.eq.j) then
+ b(jlim-i+1)=0.0d0
+ jlim1=jlim-1
+ else
+ jlim1=jlim
+ endif
+ do 3 k=i,jlim1
+ if (j.lt.k) then
+ IND=IOFFSET(NCON,j,k)
+ else
+ IND=IOFFSET(NCON,k,j)
+ endif
+ 3 b(k-i+1)=attalums(IND)
+ write (iout,1010) j,(b(k),k=1,jlim-i+1)
+ 2 continue
+ 1 continue
+ 1010 format (i5,3x,10(f6.2,1x))
+ return
+ end
+C-----------------------------------------------------------------------
+ block data
+ include 'DIMENSIONS'
+ include 'COMMON.LOCAL'
+ data dsc / 1.237, ! CYS (type 1)
+ & 2.142, ! MET (type 2)
+ & 2.299, ! PHE (type 3)
+ & 1.776, ! ILE (type 4)
+ & 1.939, ! LEU (type 5)
+ & 1.410, ! VAL (type 6)
+ & 2.605, ! TRP (type 7)
+ & 2.484, ! TYR (type 8)
+ & 0.743, ! ALA (type 9)
+ & 0.000, ! GLY (type 10)
+ & 1.393, ! THR (type 11)
+ & 1.150, ! SER (type 12)
+ & 2.240, ! GLN (type 13)
+ & 1.684, ! ASN (type 14)
+ & 2.254, ! GLU (type 15)
+ & 1.709, ! ASP (type 16)
+ & 2.113, ! HIS (type 17)
+ & 3.020, ! ARG (type 18)
+ & 2.541, ! LYS (type 19)
+ & 1.345, ! PRO (type 20)
+ & 0.000 /! D (type 21)
+ end