!     Last change:  RCD   5 May 2010    2:18 pm
!-----------------------------------------------------------------------------------------------------------
	module Set_Molecule
!
    parameter(NatMX=151)    !Maximum atoms per Molecule
    parameter(maxfn=10)	    !Maximum ASA functions per atom

    type Molecule	    		 ! ASA prepared information for a molecule
    CHARACTER (LEN=80):: Identity	 ! Label to identify the molecular structure
    Integer:: N_atoms			 ! Number of atoms
    Integer:: Z(natmx)			 ! Atomic Numbers
    integer:: N_electrons                ! Number of electrons in the molecule
    integer:: NAP_Type			 ! (1: Mulliken; 2: NO)
    double precision :: AP(natmx)	 ! Gross Atomic Populations
    Integer:: ANG_AU			 ! (0: AU; 1: Angstroms)
    double precision :: XYZ(natMX,3)     ! Cartesian Coordinates
    integer:: NASAfun(natmx)		 ! Number of ASA terms in each atom
    double precision :: Ep(natmx,maxfn)  ! ASA Exponents    (atom,function)
    double precision :: Op(natmx,maxfn)  ! ASA Coefficients (atom,function)
    double precision :: overlap_QSS	 ! Overlap Quantum Self-Similarity
    end type Molecule
!
    end module
!-----------------------------------------------------------------------------------------------------------   
	PROGRAM ASA_molecular_set_Inward_Function_Kruskal_tree_and_ordering
!-----------------------------------------------------------------------------------------------------------
	use Set_Molecule
!    
    Integer*4 MolMolMAX
    parameter(MOLMAX=70)  			!Maximum number of Molecules
    parameter(MolMolMAX=(MOLMAX*MOLMAX+MOLMAX)/2)
!
    CHARACTER (LEN=80):: GENERAL_TITLE
    CHARACTER (LEN=7):: INW_FUNCTION
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
   	 double precision :: zero, one; PARAMETER(zero=0.0d0, one=1.0d0)
!
	integer*4 :: MoloM, Nelectrons
   	integer*4 :: MM,iii,II,JJ,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 :: PREPARED_FIL, Prepared_fil_IN
    CHARACTER*70 :: QSMatrix_IN, INW_QSMatrix_IN,INW_QSMatrix_OUT
    Character*70 :: CRB_INDX_OUT, INW_CRB_INDX_OUT
    Character*70 :: KT_OUT,INW_KT_OUT
    Character*70 :: KTO_OUT,INW_KTO_OUT
    character*4 :: INW_F
!--------------------------------------------------------------------------------------------------
!----------------------------------------------------------------------------------------------------------------
        integer*4 Zorder
        Allocatable Zorder(:)
        Double precision Crack
        logical finis    
!----------------------------------------------------------------------------------------------------------------
	Type(Molecule)Mol_Set(MOLMAX)
!--------------------------------------------------------------------------------------------------	
!
!		first reads a file with the function to be computed
!
	open(313,File='Flagstaff_7.txt')
        read(313,'(a7)')INW_Function
        close(313)
!
	write(*,*)' QSM Inward Function Chosen Is:'
		if(INW_Function=='XAPE   ')then
        IFUNNY=0; write(*,'(1x,a7)')INW_Function
		else if(INW_Function=='SQRT   ')then
        IFUNNY=1; write(*,'(1x,a7)')INW_Function
        else if(INW_Function=='LOG    ')then
        IFUNNY=2; write(*,'(1x,a7)')INW_Function  
        else if(INW_Function=='ATAN   ')then
        IFUNNY=3; write(*,'(1x,a7)')INW_Function     
		else
        write(*,*)' Unknown Function! Revise file:Inward_Function'
        pause
        stop
        end if
!--------------------------------------------------------------------------------------------------
!
!		second reads a label file, where it is found
!		information on the input and output file system
!		as used in the other programs
!
		open(304,File='Labelstaff_7.txt')
!
		Read(304,*)ID_Label_Pre
        Read(304,*)ID_Label_Suf
        Read(304,*)PREPARED_FIL
        Read(304,*)QSMatrix_IN
        read(304,*)CRB_INDX_OUT
        read(304,*)KT_OUT
        read(304,*)KTO_OUT
!
        close(304)
!------------------------------------------------------------------------------
!	If the funtion is XAPE
!	Reads a batch of ASA defined molecules
!	as given by Molecule-Type Input-Output
!
!	Reads molecular data in file  Prepared_fil_IN
!
        Prepared_fil_IN = TRIM(ID_Label_Pre)//TRIM(PREPARED_FIL)//TRIM(ID_Label_Suf)//'.txt'
       	open(232,File= Prepared_fil_IN)
!
		if(IFUNNY==0)then
		open(232,File=Prepared_fil_IN)
!
                    READ(232,'(a80)')General_Title 
                    Read(232,'(a80)')ASAFUN 
                    READ(232,*)Kind_ASA
!                    
                    READ(232,*)MoloM
!                   
		    MM=0
		    do Mol=1,MoloM
                    MM=MM+Mol
		    read(232,'(i3,a80)')Molim,Mol_Set(Mol)%Identity		!Molim is the molecule order here is dummy
  		    read(232,*)Mol_Set(Mol)%N_atoms
                    read(232,*)Mol_Set(Mol)%N_electrons				 
    				do natl=1,Mol_Set(Mol)%N_atoms	     		 
    				read(232,*)nat								!It is the atom number by itself here is dummy
    				read(232,*)Mol_Set(Mol)%Z(natl)
    				read(232,*)Mol_Set(Mol)%AP(natl)
				read(232,*)(Mol_Set(Mol)%XYZ(natl,k),k=1,3)
    				read(232,*)Mol_Set(Mol)%NASAfun(natl)			 	  
    				na_boo=Mol_Set(Mol)%NASAfun(natl) 		 
    				read(232,*)(Mol_Set(Mol)%Ep(natl,mfn),mfn=1,na_boo)  		 
    				read(232,*)(Mol_Set(Mol)%Op(natl,mfn),mfn=1,na_boo)
                    end do
				read(232,*)Mol_Set(Mol)%overlap_QSS				!Position Overlap Quantum Self-Similarity
!					
            		END DO
           			close(232)
                    end if !IFUNNY==0
!------------------------------------------------------------------------------------------------------------------  
!
!       Reads the optimal similarity matrix and the Carbó Index matrix
!
    	INW_QSMatrix_IN= TRIM(ID_Label_Pre)//TRIM(QSMatrix_IN)//TRIM(ID_Label_Suf)//'.txt'
		Open(111,File=INW_QSMatrix_IN)
!
        READ(111,'(a80)')General_title
        READ(111,*)Molom,MM
!
		ALLOCATE(Zim(MM),CRB_INDX(MM),Deig(Molom))
!
!       compute inward function of the elements of the original Z QSM matrix and store
!
   	 	Read(111,*)(iii,Zim(iii),kkk=1,MM)
         CLOSE(111)
!        
             if(IFUNNY==0)then   
               do I=1,MoloM
                  II=(I*I-I)/2
                  do J=1,I
                    Nelectrons=Mol_Set(I)%N_electrons*Mol_Set(J)%N_electrons
                    III=J+II
                Zim(III)=Zim(III)/Nelectrons
                end do
                end do
                else if(IFUNNY==1)then; Zim(III)=SQRT(Zim(III))
                else if(IFUNNY==2)then; Zim(III)=LOG(Zim(III))
                else if(IFUNNY==3)then; Zim(III)=ATAN(Zim(III))
                  end if
!
!       compute Carbó Index matrix for the Inward function(Z) and save
!
                do I=1,Molom
                II=(I*I-I)/2
                III=I+II
                ZII=SQRT(Zim(III))
                do J=1,I
                JJ=(J*J+J)/2
                ZJJ=SQRT(Zim(JJ))
                IJ=J+II
                IF(I/=J)then
                   CRB_INDX(IJ)=Zim(IJ)/(ZII*ZJJ)
                   else
                   CRB_INDX(IJ)=one
                   END if
                END do
                end do
!
                MilMil=MM
!
!-------------------------------------------------------------------------------------------------------------------
!		Diagonalize similarity matrices
!		write eigenvalues of each Similarity matrix
!
		INW_F=TRIM(INW_FUNCTION)
        INW_KT_OUT= TRIM(ID_Label_Pre)//TRIM(KT_OUT)//TRIM(ID_Label_Suf)//'.txt'
        OPEN(777,FILE=INW_F//'_'//INW_KT_OUT)
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-'//INW_Function//'---Similarity Eigenvalues----------'
        WRITE(777,*)'---------------------------------------'
        Call Jadiag (Molom,milmil,Zim,Deig)
		write(777,*)INW_Function //'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,*)INW_Function//' 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,*)INW_Function//'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,*)INW_F//'-Similarity Matrix Complete Sum: ',Solmiro
        WRITE(777,*)INW_F//'-SM OutDiag Complete Sum:  ',Soloutd
        WRITE(777,*)'---------------------------------------'        
!
!	output files of QS and Carbó index Matrices are the same as input 
!				with an extra prefix label
!		which is coincident with the label flag
!	associated to the type of inward function employed
!
!
!		usual format first
!
    	INW_QSMatrix_OUT= TRIM(ID_Label_Pre)//TRIM(QSMatrix_IN)//TRIM(ID_Label_Suf)//'.txt'
		Open(111,File=INW_F//'_'//INW_QSMatrix_OUT)
        write(111,'(a80)')General_title
        write(111,*)Molom,MM
   	 	write(111,*)(iii,Zim(iii),iii=1,MM)
        CLOSE(111)
! 
      	INW_CRB_INDX_OUT= TRIM(ID_Label_Pre)//TRIM( CRB_INDX_OUT)//TRIM(ID_Label_Suf)//'.txt'
        OPEN(111,FILE=INW_F//'_'//INW_CRB_INDX_OUT)
        write(111,'(a80)')General_title
        write(111,*)Molom,MM
   	 	write(111,*)(iii,CRB_INDX(iii),iii=1,MM)
        CLOSE(111)
!
        deallocate(Zim,Deig)
	allocate(IZord(milmil,2),Kruskal(Molom),ISequoiaS(Molom,2),Zorder(milmil))
!
!--------------------------------------------------------------------------------------
!
!		Reorders the off-diagonal elements of a compact index matrix CRB(*)
!				      and returns the indices into Zorder(*)
!
		OPEN(25,FILE=INW_F//'_Ordered_'//INW_CRB_INDX_OUT)
!
		do II=1,MilMil; Izord(II,1)=0; Izord(II,2)=0; Zorder(II)=II; end do
!        
		write(25,*)INW_F//'-Ordered Carbó Index Matrix'
        DO I=1,MoloM
          II=(I*I-I)/2
          do J=1,I-1;  IJ=J+II 
        DO K=1,MoloM
          KK=(K*K-K)/2
          do L=1,K-1; KL=L+KK 
            if(KL<IJ)Cycle 
!
		 if(CRB_INDX(IJ)<CRB_INDX(KL))then
			KKK=Zorder(KL)         
			Zorder(KL)=Zorder(IJ)
         	Zorder(IJ)=KKK
         	Crack=CRB_INDX(KL)
         	CRB_INDX(KL)=CRB_INDX(IJ)
         	CRB_INDX(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 IZord(*,1) and IZord(*,2)
!	writting the whole reordered indices and matrix CRB in Filis//.txt file
!
		do i=1,MilMil
        Index=Zorder(i);finis=.false.
          do k=1,MoloM
            kk=(k*k-k)/2
            do l=1,k
              lk=l+kk; 
              if(lk==Index)then; IZord(i,1)=l; IZord(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,IZord(IJ,1),IZord(IJ,2),Zorder(IJ),CRB_INDX(IJ),IJ=1,MilMil)
    	close(25)
	deallocate(Zorder)
!
!-----------------------------------------------------------------------------------------
!
!			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,*)INW_F//'-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,*)INW_F//'-Kruskal''s Molecular Ordering'
        	WRITE(777,*)'--------------------------------------------'
			Write(777,'(25(1x,I3))')(Kruskal(I),I=1,Molom)
        	WRITE(777,*)'--------------------------------------------'

!**************************************************************************
!
        	CLOSE(777)
!
       INW_KTO_OUT= TRIM(ID_Label_Pre)//TRIM(KTO_OUT)//TRIM(ID_Label_Suf)//'.txt'
        OPEN (666,FILE=INW_F//'_'//INW_KTO_OUT)
                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(*,*)' Inward Function Transformation of QSM Success!'
    	pause
        Stop
!
!--------------------------------------------------------------------------------------------------------------------
	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
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
