!     Last change:  RCD   7 May 2010    7:18 am
!-----------------------------------------------------------------------------------------------------------
	PROGRAM ASA_molecular_set_Jacobi_Refinement
!
	integer*4 MMAX,MMMAX
	Parameter (MMax=151)
    Parameter (MMMax=(MMax*MMax+MMax)/2)
    CHARACTER (LEN=80):: GENERAL_TITLE
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
   	 double precision :: zero; PARAMETER(zero=0.0d0)
     double precision :: one; PARAMETER(one=1.0d0)
!
	integer*4 :: MoloM
    integer*4 :: MM,iii,II,JJ,IJ,kkk
	Double precision ::  Zim(MMMax),Zejr(MMMax),Deig(MMax),CRB(MMMax)
    Double precision :: SumZ,SumZejr,ZetI,ZetJ,ZetIJ,delta,ZetM,Prod,Testis,t
!	
!----------------------------------------------------------------------------------------------------------------
    Character*70 :: ID_Label_Pre,ID_Label_Suf
    Character*70 :: CRBINDX_IN,JR_CINDX_OUT
    CHARACTER*70 :: QSMatrix_IN,QSMatrix_ININ
    CHARACTER*70 :: JR_QSMatrix,JR_QSMatrix_OUT
    Character*70 :: JR_FOLLOW_UP,JR_FOLLOW_UP_OUT
!----------------------------------------------------------------------------------------------------------------
!		Reads a label file
!
	open(304,File='Labelstaff_6.txt')
!
	Read(304,*)ID_Label_Pre
        Read(304,*)ID_Label_Suf
        Read(304,*)QSMatrix_IN
        READ(304,*)JR_QSMatrix
        Read(304,*)CRBINDX_IN
!
        close(304)
!---------------------------------------------------------------------------------------
!       Reads the optimal similarity matrix 
!
        QSMatrix_ININ=TRIM(ID_Label_Pre)//TRIM(QSMatrix_IN)//TRIM(ID_Label_Suf)//'.txt'
 	Open(111,File=QSMatrix_ININ)
        READ(111,'(a80)')General_title
        READ(111,*)Molom,MM
        write(*,*)General_title
        Write(*,*)Molom,MM
!
   	Read(111,*)(iii,Zim(iii),kkk=1,MM); CLOSE(111)
	write(*,*)' QSM Z succesfully read!'
!
!-------------------------------------------------------------------------------------------------------------------
!
      JR_QSMatrix_OUT=TRIM(ID_Label_Pre)//TRIM(JR_QSMatrix)//TRIM(ID_Label_Suf)//'.txt'
      OPEN(777,FILE=JR_QSMatrix_OUT)
!
	SumZ=zero; SumZejr=zero
!
	do I=2,Molom
    II=(I*I-I)/2
    ZetI=Zim(II+I)
    Zejr(II+I)=ZetI
!
    do J=1,I-1
      JJ=(J*J-J)/2
      ZetJ=Zim(JJ+J)
      Zejr(JJ+J)=ZetJ
 !
      delta=abs(ZetI-ZetJ)/2
      ZetM=(ZetI+ZetJ)/2
      ProD=ZetM*delta
      IJ=J+II        
	ZetIJ=Zim(IJ)
        SumZ=SumZ+ZetIJ 
        Testis=(ProD-ZetIJ**2)/delta
! 
       if(Testis>zero)then
          t=ZetIJ/delta
          Zvar=2*t*(ZetM-ZetIJ*t)/(one+t*t)
          Zejr(IJ)=ZetIJ+Zvar
          else
            Zejr(IJ)=ZetIJ
            end if
!
            SumZejr=SumZejr+Zejr(IJ)
!
            end do !J
            end do !I
        SumZ=SumZ*2
        SumZejr=SumZejr*2
! 
	write(*,*)' Original Tri sum : ', SumZ
    	write(*,*)' EJR Tri sum      : ', SumZejr
!
	write(777,'(a80)')General_title
   	write(777,*)Molom,MM
	write(777,*)(kkk,Zejr(kkk),kkk=1,MM)
!
	write(*,*)' Succesful EJR treatment!' 
!
        WRITE(777,*)'---------------------------------------------------'
        WRITE(777,*)'-------EJR refined Similarity Eigenvalues----------'
        WRITE(777,*)'---------------------------------------------------'
        Call Jadiag (Molom,MM,Zejr,Deig)
        write(777,*)(I,Deig(I),I=1,Molom)
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)' Up+Down Triangle Sum : ',SumZejr
        WRITE(777,*)'---------------------------------------'
        CLOSE(777)
!
!       compute and save Carbó index matrix
!
        CRBIND_SUMEJR=zero
!
	II=0
        do I=1,Molom
        II=II+I
        CRB(II)=one
        end do
!
	do I=2,Molom
   	 II=(I*I-I)/2
    	ZetI=Zejr(II+I)
!
    	do J=1,I-1
      	JJ=(J*J-J)/2
      	ZetJ=Zejr(JJ+J)
!
      	IJ=J+II
	CRB(IJ)=Zejr(IJ)/SQRT(ZetI*ZetJ)
        CRBIND_SUMEJR=CRBIND_SUMEJR+CRB(IJ)
!
            end do !J
            end do !I
!
		close(777)
!

       JR_CINDX_OUT=TRIM(ID_Label_Pre)//TRIM(CRBINDX_IN)//TRIM(ID_Label_Suf)//'.txt'

        OPEN(333,FILE=JR_CINDX_OUT)
!
	write(333,'(a80)')General_title
   	write(333,*)Molom,MM
!
        CRBIND_SUMEJR=2*(CRBIND_SUMEJR-Molom)
	write(333,*)(kkk,CRB(kkk),kkk=1,MM)
        WRITE(333,*)'---------------------------------------'
        WRITE(333,*)' Up+Down Triangle Sum : ',CRBIND_SUMEJR
        WRITE(333,*)'---------------------------------------'
        CLOSE(333)
!
        WRITE(*,*)' Molecular Quantum Similarity Suite'
	write(*,*)' EJR Quantum Similarity Matrix Refinement ----> Success!'
	Pause
        Stop
!
 	end
!--------------------------------------------------------------------------------------------------------------------------------------------
!----------------------------------------------------------------------------------------------
      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
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
