C
C Program to cluster united-residue MCM results.
C
      implicit none
      include 'DIMENSIONS'
      include 'sizesclu.dat'
#ifdef MPI
      include "mpif.h"
      integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
      include "COMMON.MPI"
#endif
      include 'COMMON.TIME1'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.GEO'
      include 'COMMON.HEADER'
      include 'COMMON.CONTROL'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.CLUSTER'
      include 'COMMON.IOUNITS'
      include 'COMMON.FREE'
      logical printang(max_cut)
      integer printpdb(max_cut)
      integer printmol2(max_cut)
      character*240 lineh,scrachdir2d
      REAL CRIT(maxconf),MEMBR(maxconf)
      REAL CRITVAL(maxconf-1)
      INTEGER IA(maxconf),IB(maxconf)
      INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
      INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
      integer nn,ndis
      real*4 DISNN
      DIMENSION NN(maxconf),DISNN(maxconf)
      LOGICAL FLAG(maxconf)
      integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon,
     & it,ncon_work,ind1,ilen,is,ie
      double precision t1,t2,tcpu,difconf
      real diss_(maxdist)
      
      double precision varia(maxvar)
      double precision hrtime,mintime,sectime
      logical eof
      external ilen
#ifdef MPI
      call MPI_Init( IERROR )
      call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
      call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
      Master = 0
      if (ierror.gt.0) then
        write(iout,*) "SEVERE ERROR - Can't initialize MPI."
        call mpi_finalize(ierror)
        stop
      endif
      if (nprocs.gt.MaxProcs+1) then
        write (2,*) "Error - too many processors",
     &   nprocs,MaxProcs+1
        write (2,*) "Increase MaxProcs and recompile"
        call MPI_Finalize(IERROR)
        stop
      endif
#endif

      call initialize
      call openunits
      call parmread
      call read_control
      call molread
c      if (refstr) call read_ref_structure(*30)
      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

      if (nclust.gt.0) then
        PRINTANG(1)=.TRUE.
        PRINTPDB(1)=outpdb
        printmol2(1)=outmol2
        ncut=0
      else
      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
      endif
      if (ncut.gt.0) then
      write (iout,*) 'Number of cutoffs:',NCUT
      write (iout,*) 'Cutoff values:'
      DO ICUT=1,NCUT
        WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
     &    printpdb(icut),printmol2(icut)
      ENDDO
      else if (nclust.gt.0) then
      write (iout,'("Number of clusters requested",i5)') nclust
      else
      if (me.eq.Master)
     & write (iout,*) "ERROR: Either nclust or ncut must be >0"
      stop
      endif
      DO I=1,NRES-3  
        MULT(I)=1
      ENDDO
      do i=1,maxconf
        list_conf(i)=i
      enddo
      call read_coords(ncon,*20)
      write (iout,*) 'from read_coords: ncon',ncon
      
      write (iout,*) "nT",nT
      do iT=1,nT
      write (iout,*) "iT",iT
#ifdef MPI
      call work_partition(.true.,ncon)
#endif
      call probabl(iT,ncon_work,ncon,*20)

      if (ncon_work.lt.2) then
        write (iout,*) "Too few conformations; clustering skipped"
        exit
      endif
#ifdef MPI
      ndis=ncon_work*(ncon_work-1)/2
      call work_partition(.true.,ndis)
#endif

      DO I=1,NCON_work
        ICC(I)=I
      ENDDO
      WRITE (iout,'(A80)') TITEL
      t1=tcpu()
C
C CALCULATE DISTANCES
C
      call daread_ccoords(1,ncon_work)
      ind1=0
      DO I=1,NCON_work-1
        if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
        do k=1,2*nres
          do l=1,3
            c(l,k)=allcart(l,k,i)
          enddo 
        enddo
        do k=1,nres
          do l=1,3
            cref(l,k)=c(l,k)
          enddo
        enddo
        DO J=I+1,NCON_work
          IND=IOFFSET(NCON_work,I,J)
#ifdef MPI
          if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
#endif
          ind1=ind1+1
#ifdef MPI
          DISS_(IND1)=DIFCONF(I,J)
#else
          DISS(IND1)=DIFCONF(I,J)
#endif
c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
#ifdef MPI
          endif
#endif
        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'

#ifdef MPI
      call MPI_Gatherv(diss_(1),scount(me),MPI_REAL,diss(1),
     &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
      if (me.eq.master) then
#endif
      scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
      open(80,file=scrachdir2d,form='unformatted')
      do i=1,ndis
        write(80) diss(i)
      enddo
      if (punch_dist) then
        do i=1,ncon_work-1
          do j=i+1,ncon_work
            IND=IOFFSET(NCON,I,J)
            write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
     &        energy(j)-energy(i)
          enddo
        enddo
      endif
C
C Print out the RMS deviation matrix.
C
      if (print_dist) CALL DISTOUT(NCON_work)
C
C  call hierarchical clustering HC from F. Murtagh
C
      N=NCON_work
      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,*)

#ifdef DEBUG
      write (iout,*) "The TOTFREE array"
      do i=1,ncon_work
        write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
      enddo
#endif
      call flush(iout)
      CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
      LEV = N-1
      write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
      if (lev.lt.2) then
        write (iout,*) "Too few conformations to cluster."
        goto 192
      endif
      CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)

c 3/3/16 AL: added explicit number of cluters
      if (nclust.gt.0) then
        is=nclust-1
        ie=nclust-1
        icut=1
      else
        is=1
        ie=lev-1
      endif
      do i=1,maxgr
        licz(i)=0
      enddo
      icut=1
      i=is
      NGR=is+1
      do j=1,n
        licz(iclass(j,i))=licz(iclass(j,i))+1
        nconf(iclass(j,i),licz(iclass(j,i)))=j
c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
c     &    nconf(iclass(j,i),licz(iclass(j,i)))
      enddo        
c      do i=1,lev-1
      do i=is,ie
         idum=lev-i
         DO L=1,LEV
            IF (HEIGHT(L).EQ.IDUM) GOTO 190
         ENDDO
 190     IDUM=L
c         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
c     &    " icut",icut," cutoff",rcutoff(icut)
         IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
         if (nclust.le.0) 
     &    WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
          write (iout,'(a,f8.2)') 'Maximum distance found:',
     &              CRITVAL(IDUM)
          CALL SRTCLUST(ICUT,ncon_work,iT)
          CALL TRACK(ICUT)
          CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
          icut=icut+1
          if (icut.gt.ncut) goto 191
         ENDIF
         NGR=i+1
         do l=1,maxgr
          licz(l)=0
         enddo
         ii=i-is+1
         do j=1,n
          licz(iclass(j,ii))=licz(iclass(j,ii))+1
          nconf(iclass(j,ii),licz(iclass(j,ii)))=j
c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
c     &    nconf(iclass(j,i),licz(iclass(j,i)))
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.'

#ifdef MPI
      endif
#endif
 192  continue
      enddo
C
      close(icbase,status="delete")
#ifdef MPI
      call MPI_Finalize(IERROR)
#endif
      stop '********** Program terminated normally.'
   20 write (iout,*) "Error reading coordinates"
#ifdef MPI
      call MPI_Finalize(IERROR)
#endif
      stop
   30 write (iout,*) "Error reading reference structure"
#ifdef MPI
      call MPI_Finalize(IERROR)
#endif
      stop
      end
c---------------------------------------------------------------------------
      double precision function difconf(icon,jcon)
      implicit none
      include 'DIMENSIONS'
      include 'sizesclu.dat'
      include 'COMMON.CONTROL'
      include 'COMMON.CLUSTER'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      logical non_conv
      double precision przes(3),obrot(3,3)
      double precision xx(3,maxres2),yy(3,maxres2)
      integer i,ii,j,icon,jcon
      double precision rms
      if (lside) then
        ii=0
        do i=nstart,nend
          ii=ii+1
          do j=1,3
            xx(j,ii)=allcart(j,i,jcon)
            yy(j,ii)=cref(j,i)
          enddo
        enddo
        do i=nstart,nend
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
        enddo
        call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+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------------------------------------------------------------------------------
      subroutine distout(ncon)
      implicit none
      include 'DIMENSIONS'
      include 'sizesclu.dat'
      integer ncol,ncon
      parameter (ncol=10)
      include 'COMMON.IOUNITS'
      include 'COMMON.CLUSTER'
      integer i,j,k,jlim,jlim1,nlim,ind,ioffset
      real*4 b
      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)=diss(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