!     Last change:  RCD   5 May 2010    2:10 pm
!-----------------------------------------------------------------------------------------------------------
	PROGRAM ASA_molecular_set_Kruskal_tree_and_ordering
!
    parameter(NatMX=151)			!Maximum atoms per Molecule
!
    Integer*4 MolMolMAX
    parameter(MOLMAX=70)  			!Maximum number of Molecules
    parameter(MolMolMAX=(MOLMAX*MOLMAX+MOLMAX)/2)
!
    CHARACTER (LEN=80):: GENERAL_TITLE
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
   	 double precision :: zero; PARAMETER(zero=0.0d0)
!
	integer*4 :: MoloM
        integer*4 :: MM,iii,II,IJ
        Double precision  SZP(Molmax)
	Double precision  Zim,CRB_INDX,Deig
	allocatable:: zim(:),CRB_INDX(:),Deig(:)
!
	Allocatable:: Izord(:,:),Kruskal(:),ISequoiaS(:,:)
!
	Logical TaK,TaL
!----------------------------------------------------------------------------------------------------------------
    Character*70 :: ID_Label_Pre,ID_Label_Suf
    Character*70 :: KRUSKAL_TREE,KT_ORDER,CRB_INDX_ORD
    Character*70 :: QSMatrix_IN,CINDX_IN   
    Character*70 :: INPUT_QSM,INPUT_CRBINDX,OUTPUT_KMST,OUTPUT_QSM,OUTPUT_CRBINDX,OUTPUT_KTO
!----------------------------------------------------------------------------------------------------------------
!		Reads a label file, where it is found 
!			input/output file names
!
		open(304,File='Labelstaff_3.txt')
!
		Read(304,*)ID_Label_Pre
        Read(304,*)ID_Label_Suf
        Read(304,*)KRUSKAL_TREE
        READ(304,*)KT_ORDER
        Read(304,*)QSMatrix_IN
        Read(304,*)CINDX_IN
        Read(304,*)CRB_INDX_ORD
        close(304)
!
!
!       Reads the optimal similarity matrix and the Carbó Index matrix
!
		INPUT_QSM=TRIM(ID_Label_Pre)//TRIM(QSMatrix_IN)//TRIM(ID_Label_Suf)//'.txt'
 		Open(111,File=INPUT_QSM)
        READ(111,'(a80)')General_title;READ(111,*)Molom,MM
!
		ALLOCATE(Zim(MM),CRB_INDX(MM),Deig(Molom))
!
   	 	Read(111,*)(iii,Zim(iii),kkk=1,MM); CLOSE(111)
!
		INPUT_CRBINDX=TRIM(ID_Label_Pre)//TRIM(CINDX_IN)//TRIM(ID_Label_Suf)//'.txt'
		OPEN(222,FILE=INPUT_CRBINDX);READ(222,'(a80)')General_title;READ(222,*)Molom,MM
    		Read(222,*)(iii,CRB_INDX(iii),kkk=1,MM);CLOSE(222)
                MilMil=MM
!
!-------------------------------------------------------------------------------------------------------------------
!	Diagonalize similarity matrices
!	write eigenvalues of each Similarity matrix
!
		OUTPUT_KMST=TRIM(ID_Label_Pre)//TRIM(KRUSKAL_TREE)//TRIM(ID_Label_Suf)//'.txt'
        OPEN(777,FILE=OUTPUT_KMST)
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-------Similarity Eigenvalues----------'
        WRITE(777,*)'---------------------------------------'
        Call Jadiag (Molom,milmil,Zim,Deig)
		write(777,*)'	Overlap Similarity Matrix'
        write(777,*)(I,Deig(I),I=1,Molom)
        Solveig=Zero; do I=1,Molom; Solveig=Solveig+Deig(I); end do
        write(777,*)'   Similarity Matrix Eigenvalues Sum:',Solveig
        WRITE(777,*)'---------------------------------------'
		Solveig=Zero;MM=0; do I=1,Molom; MM=MM+I; Solveig=Solveig+Zim(MM); end do
        write(777,*)'   Similarity Matrix  Trace:',Solveig!
        WRITE(777,*)'---------------------------------------'
        Solmiro=Zero
        do I=1,Molom; II=(I*I-I)/2
        DO J=1,I; IJ=J+II
        ZUM=Zim(IJ)
        IF(I/=J)ZUM=ZUM+ZUM
        Solmiro=Solmiro+ZUM
        END DO
        End Do
        Soloutd=Solmiro-Solveig
		write(777,*)'   Similarity Matrix Complete Sum: ',Solmiro
        WRITE(777,*)'   SM OutDiag Complete Sum:  ',Soloutd
        WRITE(777,*)'---------------------------------------'        
!
!		Save the Optimal Overlap Matrix on file Z_OVERLAP.txt
!		and the Carbó Index one on file R_CarboIndx.txt
!
		OUTPUT_QSM=TRIM(ID_Label_Pre)//'OUT_QSM'//TRIM(ID_Label_Suf)//'.txt'
        OUTPUT_CRBINDX=TRIM(ID_Label_Pre)//'OUT_CRBINDX'//TRIM(ID_Label_Suf)//'.txt'
        OPEN(34,FILE=OUTPUT_QSM)
        OPEN(35,FILE=OUTPUT_CRBINDX)
        WRITE(34,'(A80)')GENERAL_TITLE
        write(34,*)Molom,milmil
        WRITE(35,'(A80)')GENERAL_TITLE
        write(35,*)Molom,milmil
        DO I=1,Molom
        II=(I*I-I)/2 
        Write(34,*)I
        Write(35,*)I
        do J=1,I; IJ=J+II ; Deig(J)=Zim(IJ); SZP(J)=CRB_INDX(IJ);end do
        write(34,*)(J,Deig(J),J=1,I)
        write(35,*)(J,SZP(J),J=1,I)
        END DO
        Close(34)
        Close(35)
!
        deallocate(Zim,Deig)
		allocate(IZord(milmil,2),Kruskal(Molom),ISequoiaS(Molom,2))
!
!
!		Order Carbó Similarity Index Matrix 
!
		call Ordering_CRB(Molom,milmil,IZord,CRB_INDX,ID_Label_Pre,CRB_INDX_ORD,ID_Label_Suf)
!
!			IZord contains now the old off diagonal indexes of the ordered 
!		CRB_INDX elements all this information is written in the file Z_Auxili.txt
!			Then start the coherent similarity matrix optimization
!		  Stablish first a new sequence of molecules to be procesed 
!
!------------------------------------------------------------------------------------
!
!		{Order(IJ,1);Order(IJ,2} contains the (I,J) pairs of elements in 
!		a distance-like matrixordered from shortest to largest distances or  
!		from a similarity index matrix from largest to lesser values.
!		(S(P,1);S(P,2))(P=1,M-1) contains the Kruskal's tree links.
!
!				Kruskal's minimum spanning tree: 
!		Source: J. A. Hartigan "Clustering Algorithms" 
!				John Wiley & Sons (New York) (1975) pp 201.
!		Origin: J. B. Kruskal "On the shortest spaning tree 
!				of a graph, and a travelling salesman problem"
!				Proc. Amer. Math. Soc. 7 (1956) 48-50.
!		See also the algorithm described in:
!				R. Carbó, Ll. Domingo "LCAO-MO Similarity Measures and Taxonomy"
!				Intl. J. Quantum Chem. 32 (1987) 517-545.
!------------------------------------------------------------------------------------------
	   do I=1,Molom
             ISequoias(I,1)=0
              ISequoiaS(I,2)=0
               End do
!
        	IZI=1 
            ISequoiaS(1,1)=IZord(2,1)
            ISequoiaS(1,2)=IZord(2,2)
!            
 	    do I=1,milmil
            K=IZord(I,1)
            L=IZord(I,2)
!            
!			This is here because of the diagonal (I,I) pairs
!
	            if(K==L)Cycle	
!
                TaK=.False. !K is not into the tree
		do KL=1,IZI
             if(ISequoiaS(KL,1)==K.OR.ISequoiaS(KL,2)==K)TaK=.True. !K is into the tree
                end do
                TaL=.False. !L is not into the tree
        do KL=1,IZI
             if(ISequoiaS(KL,1)==L.OR.ISequoiaS(KL,2)==L)TaL=.True. !L is into the tree 
                end do
!
!				To admit the pair (K,L) as two nodes into Kruskal 
!				tree, one of the two molecules has to be allready  
!				into the tree node list and the other has to be 
!						not yet into the tree node list.
!
				if((TaK.and..not.TaL).or.(.not.TaK.and.TaL))then  				
!
              	IZI=IZI+1 
				ISequoiaS(IZI,1)=K
                ISequoias(IZI,2)=L
                If(IZI==Molom-1) Exit	!When M-1 links are reached the tree is done
!
                end if  
!
			end do             
!
!------------------------------------------------------------------------------------
!
!	 {ISequoiaS(*,1);ISequoiaS(*,2)} are the molecular index pairs in the new order
!			     Given by Kruskal's minimum spanning tree        
!
        WRITE(777,*)'---------------------------------------------'
		Write(777,*)'   Kruskal''s Minimum Spanning Tree Links'
        WRITE(777,*)'---------------------------------------------'
		do i=1,Molom-1		
		Write(777,*)(ISequoiaS(i,k),k=1,2)
        	End do
        WRITE(777,*)'---------------------------------------------'
!
!		Use Kruskal tree index ordering for coherent optimization.
!			  Construct the new molecular sequence first
!
		do i=1,Molom; Kruskal(i)=0; end do
!
		Klux=0
		do i=1,Molom-1
        	K=ISequoiaS(i,1); L=ISequoiaS(i,2)
            TaK=.False.; TaL=.False.
          		do J=1,Klux 
            	If(Kruskal(J)==K)TaK=.True.
            	If(Kruskal(J)==L)TaL=.True.
            	end do
		    if(.not.TaK)then; Klux=Klux+1; Kruskal(Klux)=K; end if
            if(.not.TaL)then; Klux=Klux+1; Kruskal(Klux)=L; end if
            if(Klux==Molom)exit        
        end do
!
        WRITE(777,*)'---------------------------------------'
		Write(777,*)'   Kruskal''s Molecular Ordering'
        WRITE(777,*)'---------------------------------------'
		Write(777,'(25(1x,I3))')(Kruskal(I),I=1,Molom)
        WRITE(777,*)'---------------------------------------'

!**************************************************************************
!
        CLOSE(777)
!
		OUTPUT_KTO=TRIM(ID_Label_Pre)//TRIM(KT_ORDER)//TRIM(ID_Label_Suf)//'.txt'
        OPEN (666,FILE=OUTPUT_KTO)
                WRITE(666,'(a80)')General_title
                WRITE(666,*)Molom
        	Write(666,'(25(1x,I3))')(Kruskal(I),I=1,Molom)
        CLOSE(666)
!
                WRITE(*,*)' Molecular QUantum Similarity Suite'
		write(*,*)' Kruskal Tree and Order Success!'
		pause        
!
!--------------------------------------------------------------------------------------------------------------------
	end                !program
!...............................................................................................................
!..................................................................................................
!--------------------------------------------------------------------------
!----------------------------------------------------------------------------------------------
      SUBROUTINE JADIAG(N,NN,O,D)
!-----------------------------------------------------------------------
!     SYMMETRIC MATRIX DIAGONALIZATION BY JACOBI'S METHOD
!-----------------------------------------------------------------------
!     N         MATRIX DIMENSION
!     O(NN)     MATRIX TO BE DIAGONALIZED
!     V(N,N)    EIGENVECTOR MATRIX
!     D(N)      EIGENVALUE MATRIX
!     NN        (N*N+N)/2
!-----------------------------------------------------------------------
	Implicit Double precision (A-H,O-Z), integer*4 (i-n)
	integer*4 p,q
        LOGICAL la,lb,lc
!
	parameter(zero=0.0d0,one=1.0d0)
	parameter(half=.50d0,potwo=.20d0)
	parameter(tenm8=1.0d-8,tenm10=1.0d-10)
!
	dimension O(nn),D(n)
	allocatable:: A(:), Z(:), V(:,:)
!
	IF (N.LE.0) STOP ' Bad dimension in JADIAG'
!
! Transfers matrix
!
! preserves matrix to be diagonalized
!
        ALLOCATE (A(NN),Z(N),V(N,N))
!
	do j=1,NN; A(j)=O(j); end do
!
	N1=N-1
	IPP=0
	do P=1,N
   	do Q=1,N
      	V(P,Q)=ZERO
   	end do
   	V(P,P)=ONE
   	IPP=IPP+P
   	D(P)=A(IPP)
   	Z(P)=ZERO
	end do
!
	if (N==1) RETURN
	iter=0
	sm=one
	do while(sm>tenm8.and.iter<100)
	   iter=iter+1
!
! CONVERGENCE CRITERION
!
	   ipq=1
   		SM=ZERO
   		DO P=2,N
      		DO Q=1,P-1
        	 ipq=ipq+1
         	SM=SM+ABS(A(IPQ))
      		end do
     		 ipq=ipq+1
   		end do
   	TRESH=ZERO
   	IF (Iter<4) TRESH=POTWO*SM/(N*N)
!
! GET OFF-DIAGONAL ELEMENTS BELOW THRESHOLD
!
 	  ipp=0
   	do P=1,N1
      iqq=ipp+p
      do Q=P+1,N
         IPQ=iqq+p
         G=100*ABS(A(IPQ))
         la=Iter>4
         lb=abs(ABS(D(P))+G)-ABS(D(P))<tenm10
         lc=abs(ABS(D(Q))+G)-ABS(D(Q))<tenm10
         IF (la.and.lb.and.lc) then
            A(IPQ)=ZERO
         else
            IF (ABS(A(IPQ))>TRESH) then
               H=D(Q)-D(P)
               IF ((ABS(H)+G)<=ABS(H)) then
!
! GENERATE THE SHIFT FOR THE JACOBI ROTATION
!
                  T=A(IPQ)/H
               else
                  THETA=HALF*H/A(IPQ)
                  T=ONE/(ABS(THETA)+SQRT(ONE+THETA*THETA))
                  IF (THETA<ZERO)T=-T
               end if
               C=ONE/SQRT(ONE+T*T)
               S=T*C
               TAU=S/(ONE+C)
               H=T*A(IPQ)
               Z(P)=Z(P)-H
               Z(Q)=Z(Q)+H
               D(P)=D(P)-H
               D(Q)=D(Q)+H
               A(IPQ)=ZERO
               jj=0
               do J=1,N
                  if (j/=p) then
                     if (j<p) then
                        jp=ipp+j
                     else
                        jp=jj+p
                     endif
                     if (j/=q) then
                        if (j<q) then
                           jq=iqq+j
                        else
                           jq=jj+q
                        endif
                        G=A(JP)
                        H=A(JQ)
                        A(JP)=G-S*(H+G*TAU)
                        A(JQ)=H+S*(G-H*TAU)
                     endif
                  endif
                  jj=jj+j
!
! STORE ORTHOGONAL TRANSFORMATION IN EIGENVECTOR MATRIX
!
                  G=V(J,P)
                  H=V(J,Q)
                  V(J,P)=G-S*(H+G*TAU)
                  V(J,Q)=H+S*(G-H*TAU)
               end do
            endif
         endif
         iqq=iqq+q
      end do
      ipp=ipp+p
 	  end do
!
! ADD INCREMENTS TO DIAGONAL ELEMENTS TO PRESERVE PRECISION
!
 	  IPP=0
   	do P=1,N
      IPP=IPP+P
      D(P)=A(IPP)+Z(P)
      A(IPP)=D(P)
      Z(P)=ZERO
   	end do
!
 	end do ! while (sm>toler.and.i<100)
!
! ORDER EIGENVALUES AND EIGENVECTORS from lesser to greater eigenvalue
!
	DO P=1,N
   	DO Q=P,N
      	if (D(Q)<D(P)) then
         G=D(P)
         D(P)=D(Q)
         D(Q)=G
         DO J=1,N
            G=V(j,p)
            V(j,p)=V(j,q)
            V(j,q)=G
         end do
         endif
   	end do
	end do
!
        DEALLOCATE(A,Z,V)
!
	if(iter>=90)then
    write(*,*)'***----Jacobi Diagonalization Warning----***'
    write(*,*)'***----Procedure Sweep Number is large---***=',iter
    end if
!

	END subroutine JADIAG
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
		subroutine Ordering_CRB(M,mm,Order,CRB,FILPRE,Filis,FILSUB)
!
		implicit double precision(a-h,o-z),integer*4(i-n)
       	Dimension CRB(mm)
        integer*4 Order(mm,2),Zorder(mm)
        CHARACTER (LEN=70) FILPRE,Filis,FILSUB,Fillas
        logical finis
!
!		Reorders the off-diagonal elements of a compact index matrix CRB(*)
!				      and returns the indices into Order(*,3)
!
		Fillas=TRIM(FILPRE)//TRIM(Filis)//TRIM(FILSUB)//'.txt'
		OPEN(25,FILE=Fillas)
!
		do II=1,mm; Order(II,1)=0; Order(II,2)=0; Zorder(II)=II; end do
!        
		write(25,*)'	Ordered Carbó Index Matrix'
        DO I=1,M
          II=(I*I-I)/2
          do J=1,I-1;  IJ=J+II 
        DO K=1,M
          KK=(K*K-K)/2
          do L=1,K-1; KL=L+KK 
            if(KL<IJ)Cycle 
!
		 if(CRB(IJ)<CRB(KL))then
			KKK=Zorder(KL)         
			Zorder(KL)=Zorder(IJ)
         	Zorder(IJ)=KKK
         	Crack=CRB(KL)
         	CRB(KL)=CRB(IJ)
         	CRB(IJ)=Crack
       	 end if
!          
         end do
         END DO     
        end do
        END DO		
!
!		 Then get the row-column indexes associated to the new order
!		     return them into the matrix Order(*,1) and Order(*,2)
!	writting the whole reordered indices and matrix CRB in Filis//.txt file
!
		do i=1,mm
        Index=Zorder(i);finis=.false.
          do k=1,M
            kk=(k*k-k)/2
            do l=1,k
              lk=l+kk; 
              if(lk==Index)then; Order(i,1)=l; Order(i,2)=k; finis=.true.; exit; end if
            end do	!l
            if(finis)exit
          end do  !k
        end do	!i
!        
	write(25,'(5(1x,i3,1x,i3,1x,i3,1x,i3,e15.8))')(IJ,Order(IJ,1),Order(IJ,2),Zorder(IJ),CRB(IJ),IJ=1,mm)
        close(25)
!
		end
!-----------------------------------------------------------------------------------------
!-!
!--------------------------------------------------------------------------------------
