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
      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,scount_buf
      real*4 DISNN, diss_buf(maxdist)
      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,kkk, ijk, is,ie
      double precision t1,t2,tcpu,difconf
      
      double precision varia(maxvar)
      double precision hrtime,mintime,sectime
      logical eof
#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 cinfo
      call read_control
      call parmread
      call molread
c      write (iout,*) "Main: refstr ",refstr
      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
      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,*) "Temperature",1.0d0/(beta_h(iT)*1.987D-3)
#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
c        if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
        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
          DISS(IND1)=DIFCONF(I,J)
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()
c      PRINT '(a)','End of distance computation'

      scount_buf=scount(me)

      do ijk=1, ndis
      diss_buf(ijk)=diss(ijk)
      enddo


#ifdef MPI
      WRITE (iout,*) "Wchodze do call MPI_Gatherv"
      call MPI_Gatherv(diss_buf(1),scount_buf,MPI_REAL,diss(1),
     &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
      if (me.eq.master) then
#endif
      open(80,file=scratchdir(:ilen(scratchdir))//'/distance',
     &   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
         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
     &    " 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
         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)))
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'
      integer ipermmin
      double precision przes(3),obrot(3,3)
      double precision rmscalc
      integer icon,jcon,k,l
c      write (iout,*) "DIFCONF: ICON",icon," JCON",jcon
      do k=1,2*nres
        do l=1,3
          cref(l,k)=allcart(l,k,icon)
          c(l,k)=allcart(l,k,jcon)
        enddo 
      enddo
      difconf=rmscalc(c(1,1),cref(1,1),przes,obrot,ipermmin)
      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