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)
      
      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
      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
C      
C CALCULATE DISTANCES
C
   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
      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()
      DO I=1,NCON-1
        if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
        if (from_cart) 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)
      LEV = N-1
      CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)

      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)
      do i=1,nres
        phi(i)=phiall(i,jcon)
        theta(i)=thetall(i,jcon)
        alph(i)=alphall(i,jcon)
        omeg(i)=omall(i,jcon)
      enddo
      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)=cref(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(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
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,c(1,nnt),cref(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
      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