!     Last change:  RCD   6 May 2010   12:40 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_Circular_and_Random_Permutations_Quantum_Similarity_Matrices
	use set_molecule
!
    Integer*4 MolMolMAX
    parameter(MOLMAX=160)  			!Maximum number of Molecules
    parameter(MolMolMAX=(MOLMAX*MOLMAX+MOLMAX)/2)
!
    DOUBLE PRECISION:: zero, one 
    PARAMETER (zero=0.0d0, one=1.0d0)
    double precision::  toler
    PARAMETER (toler=1.0D-7)
!
    CHARACTER (LEN=80):: GENERAL_TITLE
    CHARACTER (LEN=12):: ASAFUN		!Name of the file for ASA function typer
    CHARACTER (LEN=7):: Kind_perm
!-------------------------------------------------------------------------------------------------
        Logical:: LLORO						! .true. write more info
        logical:: Zenith					! .true. constructs file with optimal coordinates
        LOGICAL:: Azimut                   			! .true. Makes bold optimization
        LOGICAL:: Moldone(MOLMAX)				!  To know which molecules are allready done				! .true. write more info
!-------------------------------------------------------------------------------------------------
	double precision :: Sobmax, SSmx
   	double precision :: tran(3),xyzta(NatMX,3),xyztb(NatMX,3)
   	double precision :: xyztaa(NatMX,3),xyztbb(NatMX,3)
    	double precision :: xa(3),xb(3)
	double precision :: ra(3),caa,saa,cza,sza
	double precision :: rb(3),cab,sab,czb,szb
	double precision :: TraceZ
   	double precision :: SOVLPM,Scaled,rpm,distpm
   	DOUBLE PRECISION :: Optxyz
!
	Integer :: CirPer, RanPer
!
	integer zum(MOLMAX)
	integer*4 :: MoloM, natmax, nbtmax 
   	integer*4 :: MM,iii,MilMil,KTOrder
	double precision :: ZZ, CRB_INDX
   	ALLOCATABLE  ZZ(:),CRB_INDX(:),KTOrder(:),CirPer(:,:),Optxyz(:,:,:),RanPer(:,:)
!
    integer IXYZA(0:7,3),IXYZB(0:7,3)
!    
    allocatable ira(:),done(:)
    logical :: done
!	Silverfrost fortran    
!    real :: Random@
    real :: Rn1,Delta,Rick
!    
!--------------------------------------------------------------------------------------------------	
	Type(Molecule)Mol_Set(MOLMAX)
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
    Character*70 :: ID_Label_Pre,ID_Label_Suf
    Character*70 :: PREPARED_FIL, Prepared_fil_IN
    Character*70 :: OUTPUT_FIL, FOLLOW_UP
    CHARACTER*70 :: QSMatrix_OUT, CR_QSMatrix_OUT
    Character*70 :: CRBINDX_OUT, CR_CRBINDX_OUT
    Character*70 :: OPT_XYZ_OUT, CR_OPT_XYZ_OUT
    CHARACTER*70 :: CRICRI
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!
!       read flags in Flagstaff_4b.txt
!
	LLORO=.false.
    	Zenith=.false.
    	Azimut=.false.
        OPEN(333,FILE='Flagstaff_5.txt')
        READ(333,*)iii
        IF(iii==1)LLoro=.true.
        READ(333,*)III
        IF(III==1)Zenith=.true.
        READ(333,*)iii
        IF(iii==1)Azimut=.true.
!
        read(333,'(a7)')Kind_perm
	if(Kind_perm=='CIRPERM')then
        KPindx=1
        CRICRI='CIRC_PERM_'
        else if(Kind_perm=='RANPERM')then
        KPindx=2
        CRICRI='RAND_PERM_'
        READ(333,*)Nact_perm
	if(Nact_perm<=0)Nact_perm=MOLMAX
	  else
!
      write(*,*)' ***++**Error in Permutation Kind choice!***'
      write(*,*)' Please, correctFlagstaff_5.txt contents!***'
      write(*,*)' **Only Keywords: RANPERM or CIRPERM allowed'
      pause
      STOP
!
        END if
!
        CLOSE(333)
!
!		second reads a label file, where it is found
!		information on the file system used in the
!		preceeding Molecular Quantum Similarity Suite 
!
		open(304,File='Labelstaff_5.txt')
!
		Read(304,*)ID_Label_Pre
        Read(304,*)ID_Label_Suf
        Read(304,*)PREPARED_FIL
        Read(304,*)OUTPUT_FIL
        Read(304,*)QSMatrix_OUT
        READ(304,*)CRBINDX_OUT
        READ(304,*)OPT_XYZ_OUT
!
        close(304)
!------------------------------------------------------------------------------
!	Reads a batch of ASA defined molecules
!	as given by Molecule-Type 1 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)
!
                    READ(232,'(a80)')General_Title 
                    Read(232,'(a12)')ASAFUN 
                    READ(232,*)Kind_ASA
!
!     once the molecular set dimension is known allocate working matrices
!              
                    READ(232,*)MoloM
                    MilMil=(Molom+Molom*Molom)/2
                    ALLOCATE(ZZ(MilMil),CRB_INDX(MilMil),KTOrder(MoloM))
!                                        
		    MM=0; Maxatom=0 ; TraceZ=zero
		    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
                    IF(Mol_Set(Mol)%N_atoms>Maxatom)Maxatom=Mol_Set(Mol)%N_atoms    !Maxatom is the maximal number of 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		! Overlap Quantum Self-Similarity
                                TraceZ=TraceZ+Mol_Set(Mol)%overlap_QSS
!
            		END DO
           	close(232)
!
!	construct simple XYZ rotation information
!
	do KX=0,1
         do LX=0,1
           do MX=0,1
 	    INX=KX+2*LX+4*MX
            IXYZA(INX,1)=KX
            IXYZA(INX,2)=LX
            IXYZA(INX,3)=MX
            IXYZB(INX,1)=KX
            IXYZB(INX,2)=LX
            IXYZB(INX,3)=MX
           END DO 
	 END DO
         END DO
!
!	Molecular set data retrieved succesfully and in memory now
!      	    Construct Circular or Random Permutations of Molecular indices
!              The original order in the QSM will be kept
!

      IF(KPindx==1)then
      Allocate(CirPer(MoloM,MoloM))
	Nact_perm=MoloM
      do ICP=1,MoloM; IICP=ICP-1
        do JCP=1,MoloM-IICP
          JJCp=JCP+IICP
        CirPer(ICP,JCP)=JJCP
        end do
        do JCP=1,IICP
          JJCP=MoloM-IICP+JCP
          CirPer(ICP,JJCP)=JCP
          end do
          end do
!          
        else IF(KPindx==2)then
        Allocate(RanPer(Nact_perm,MoloM))
!
!		constructs the active random permutations
!
	Rn1=MoloM+1
   	 Delta=1.0d0/Rn1
!    
	allocate(ira(Molom),done(Molom))
!
!       Silverfrost Fortran   95
!
!	call  DATE_TIME_SEED@
!
!       Lahey-Fujitsu Fortran 95
!
    	call  random_SEED
!
	do I=1,Nact_Perm
    done=.false.
    JN=0
    do
!
!       Silverfrost Fortran   95
!
!      Rick=random@()
!
!       Lahey-Fujitsu Fortran 95
!
	 call random_number(Rick)
!
      K=Rick/Delta
      if(K==0)Cycle
      if(.not.done(K))then
        done(K)=.true.
        JN=JN+1
        ira(JN)=K
        end if
        if(JN==MoloM)exit
      end do
      do J=1,MoloM; RanPer(I,J)=ira(J); end do
   	end do
!    
	deallocate(ira,done)
!
	end if  !KPindx
!
!	Allocate the optimal coordinates
!
	Allocate(optxyz(MoloM,Maxatom,3))
!    

        FOLLOW_UP = TRIM(ID_Label_Pre)//TRIM(CRICRI)//TRIM(OUTPUT_FIL)//TRIM(ID_Label_Suf)//'.txt'
              	Open(777,File= FOLLOW_UP)
                IF(KPindx==1)then
    		WRITE(777,*)' OPTIMAL CIRCULAR PERMUTATIONS OF MOLECULAR INDICES WITH FROZEN COORDINATES'
                else if(KPindx==2)then
 	        WRITE(777,*)' OPTIMAL RANDOM PERMUTATIONS OF MOLECULAR INDICES WITH FROZEN COORDINATES'
                end if
!
            Write(777,'(a80)')General_Title
            Write(777,'(a12)')ASAFUN
            Write(777,*)Kind_ASA
!
	SumZMX=Zero;MXCRIN=0
        do NICP=1,Nact_Perm
!
        do I=1,MoloM; 
          select case(KPindx)
            case(1)
          KTOrder(I)=CirPer(NICP,I)
  			case(2) 
          KTOrder(I)=RanPer(NICP,I)
          end select     
        END do
!
!       put first ordered diagonal elements
!
        do Mol=1,MoloM
             Molord=KTOrder(Mol)
             MM=(Molord*Molord+Molord)/2
             Moldone(Molord)=.false.
             ZZ(MM)=Mol_Set(Molord)%overlap_QSS			!Put into the diagonal of the similarity matrix
             CRB_INDX(MM)=one
         end do
				
!-------------------------------------------------------------------------------------------------------------------
!	Starting Sinilarity Matrix Computation with Aufbau in mind
!-------------------------------------------------------------------------------------------------------------------
!      			Overlap Matrix in this case
!-------------------------------------------------------------------------------------------------------------------
!
        WRITE(777,*)' Circular Permutation:',NICP
        WRITE(777,*)(KTOrder(I),I=1,MoloM)
!
!	put original coordinates in OPTXYZ
!
		do Mol=1,MoloM
			MolA=KTOrder(Mol)
             		do ia=1,Mol_Set(MolA)%N_atoms
               	do i=1,3; Optxyz(MolA,ia,i)=Mol_Set(MolA)%XYZ(ia,i); end do
        	END do
        ENd do
!
!       Start bimolecular quantum similarity measures first translation only
!
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'----------Bimolecular QSM--------------'
        WRITE(777,*)'---------------------------------------'
!
        do Mol=2,Molom
        MolA=KTorder(Mol)   !KTOrder contains original indices Kruskal Tree ordered or in other arbitrary order
        do Mom=1,Mol-1
! *
!       Take the ordering of Kruskal Tree or other
!
        MomB=KTOrder(Mom)     !KTOrder contains original indices Kruskal Tree ordered or in other arbitrary order
        Scaling=Mol_Set(MolA)%N_electrons*Mol_Set(MomB)%N_electrons
!
        IF(MolA<=MomB)then
        Molmom=MolA+(MomB*MomB-MomB)/2
        else
        Molmom=MomB+(MolA*MolA-MolA)/2
        END if
!
!       test the optimization done
!
        do i=1,MoloM; IF(moldone(i))THEN; zum(i)=1;ELSE;zum(i)=0; END IF; end do
        WRITE(*,'(I4,I4,40(1x,I1))')MolA,MomB,(zum(i),i=1,Molom)
!
	WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------Molecules: '
        write(777,'(1x,"-->",i4,"-*-",a80)')MolA,Mol_Set(MolA)%Identity
        write(777,'(1x,"-->",i4,"-*-",a80)')MomB,Mol_Set(MomB)%Identity
        WRITE(777,*)'---------------------------------------'
!
!       Move heavy atoms of Molecule A over the heavy atoms of molecule B
!
		Sobmax=zero; natmax=0; nbtmax=0
!
		if(LLORO)then
        WRITE(777,*)'---------------------------------------------------'
        WRITE(777,*)'--------------Atomic Translations------------------'
        WRITE(777,*)'---------------------------------------------------'
		end if	!LLORO
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!<************
	IF(.not.moldone(MolA).and..not.moldone(MomB))then; NRT_AB=0;
        else if(.not.Moldone(MolA).and.Moldone(MomB))then; NRT_AB=1;
        else if(Moldone(MolA).and..not.Moldone(MomB))then; NRT_AB=2;
        else if(Moldone(MolA).and.Moldone(MomB))then; NRT_AB=3;
        end if  
!
		select case(NRT_AB)
!
		case(0)			!both molecules have to be TR        
!
        do iata=1,Mol_Set(MolA)%N_atoms; naz=Mol_Set(MolA)%Z(iata)
!
        do ibtb=1,Mol_Set(MomB)%N_atoms; nbz=Mol_Set(MomB)%Z(ibtb)
!
!		Before anything try to find out simple optimal positioning
!
		Sobmax=zero
!
!		initialize spare coordinates
!
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 
        xyztaa(j,itran)=Mol_Set(MolA)%XYZ(j,itran)
	xyzta(j,itran)=Mol_Set(MolA)%XYZ(j,itran)
	end do
        end do
!
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 ; xyztbb(j,itran)=Mol_Set(MomB)%XYZ(j,itran)
        end do
        end do
!
	if(Azimut)then
        Do JNX=0,7
!
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 
           If (IXYZB(JNX,ITRAN)==0)then
             xyztb(j,itran)=Mol_Set(MomB)%XYZ(j,itran)
             else
              xyztb(j,itran)=-Mol_Set(MomB)%XYZ(j,itran) 
              end if
              end do
        end do
!
!	compute position overlap similarity element with simply rotated coordinates
!
		Call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!    
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM
!
!	Put spare coordinates into maximal positioning
!
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 ; xyztaa(j,itran)=xyzta(j,itran)
        end do
        end do
!
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 ; xyztbb(j,itran)=xyztb(j,itran)
        end do
        end do
!        
        END IF
!
	end do  !JNX
        END if  !Azimut
!
        IF(naz>1.and.nbz>1)then							! Just when both atoms are heavy
!
!       move molecule B over A put atom Ibtb on Iata
!
	if(LLORO)then
        WRITE(777,*)'---------------------------------------------------'
        WRITE(777,*)' translation of atom B: ',ibtb,' over atom A: ',iata
        WRITE(777,*)'---------------------------------------------------'
        end if	!LLORO
!
!		START TRANSLATING-ROTATING MOLECULAR FRAMES For Similarity Integral Optimal Search
!
!	    form the global orientation vectors on A and B
!		          use non-hydrogen atoms only
!
        do itran=1,3 ; tran(itran)=xyztaa(iata,itran);end do   
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 ; xyzta(j,itran)=xyztaa(j,itran)-tran(itran); end do
        end do
!
!		compute alignment coordinates sum vector using heavy atoms only
!
        xa=zero                    
        do ia=1,Mol_Set(MolA)%N_atoms; nza=Mol_Set(MolA)%Z(ia)
        IF(nza>1)then
        do i=1,3; xa(i)=xa(i)+xyzta(ia,i); end do
        END if
        END do
!
!		search for rotations to align sum vectors to X axis
!
		call zalign(xa,ra,caa,saa,cza,sza)
!
!		rotate ALL atomic coordinates of molecule A
!
        do ia=1,Mol_Set(MolA)%N_atoms
        	do i=1,3; xa(i)=xyzta(ia,i); end do 
				call malign(xa,ra,caa,saa,cza,sza)
        	do i=1,3; xyzta(ia,i)=ra(i); end do 
        END do
!
        do itran=1,3 ; tran(itran)=xyztbb(ibtb,itran);end do
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 ; xyztb(j,itran)=xyztbb(j,itran)-tran(itran); end do
        end do
!
        xb=zero
        do ib=1,Mol_Set(MomB)%N_atoms; nzb=Mol_Set(MomB)%Z(ib)
        IF(nzb>1)then
        do i=1,3; xb(i)=xb(i)+xyztb(ib,i); end do
        END if
        END do
!
!		search for rotations to align sum vectors to X axis
!
		call zalign(xb,rb,cab,sab,czb,szb)
!
        do ib=1,Mol_Set(MomB)%N_atoms
        	do i=1,3; xb(i)=xyztb(ib,i); end do 
 				call malign(xb,rb,cab,sab,czb,szb)
        	do i=1,3; xyztb(ib,i)=rb(i); end do 
        END do
!
		call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM; natmax=iata; nbtmax=ibtb
!
!		Store Optimal A TR Coordinates if not done yet
!
        do j=1,Mol_Set(MolA)%N_atoms; do k=1,3; Optxyz(MolA,j,k)=xyzta(j,k);END DO; END do
        do j=1,Mol_Set(MomB)%N_atoms; do k=1,3; Optxyz(MomB,j,k)=xyztb(j,k);END DO; END do
!        
        END IF !SOVLPM>Sobmax
!
        if(LLORO)then
		WRITE(777,*)'---------------------------------------------------'
		write(777,*)'----- Translated-Rotated Molecules-----------------'
       	        WRITE(777,*)' Overlap               : ',SOVLPM
                Scaled=SOVLPM/Scaling
		WRITE(777,*)'        Scaled by NA*NB: ',Scaled
                IF(SOVLPM>Sobmax)THEN;
                WRITE(777,*)'----- Succesful Rotation-Translation Maximization!-'
                write(777,*)'----- New Overlap Value Used in Similarity Matrix--'
                WRITE(777,*)'---------------------------------------------------'
		END IF
	END if !LLORO
!
	End if  !Natoms A and B

        END do  !ibtb
        END do  !iata
!
!       This is to do not optimize anymore A and B coordinates 
!			in case they are allready optimallly aligned
!
        Moldone(MolA)=.True.
        Moldone(MomB)=.True.
!
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------*OPTIMAL RESULTS*-----------'
        WRITE(777,*)'---------------------------------------'
!
		if(LLORO)then
        WRITE(777,*)'-------for Molecules: '
        write(777,'(1x,"-->",i4,"-*-",a80)')MolA,Mol_Set(MolA)%Identity
        write(777,'(1x,"-->",i4,"-*-",a80)')MomB,Mol_Set(MomB)%Identity
        WRITE(777,*)'---------------------------------------'
		END IF	!LLORO
!
       IF(natmax/=0.and.nbtmax/=0)then
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)'MAXIMUM Position Overlap Similarity (Translated-Rotated)'
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)' Atoms A - B Displacement    : ',natmax,nbtmax
        WRITE(777,*)' Overlap                     : ',Sobmax
        IF(LLORO)then
        SSmx=Sobmax/Scaling
        WRITE(777,*)' Overlap  scaled             : ',SSmx
        rpm=Sobmax/SQRT(Mol_Set(MolA)%overlap_QSS*Mol_Set(MomB)%overlap_QSS)
        WRITE(777,*)' Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(Mol_Set(MolA)%overlap_QSS+Mol_Set(MomB)%overlap_QSS-2*Sobmax)
        WRITE(777,*)' Overlap Distance Index A-B  : ',distpm
        WRITE(777,*)'---------------------------------------'
        END if !LLORO
!
!	Form the similarity matrices and preserve optimal TR coordinates 
!				for every molecular pair
!
    	ZZ(molmom)=Sobmax
    	CRB_INDX(molmom)=rpm
!
        else
        WRITE(777,*)'++++++++++> NO OPTIMAL VALUE FOUND!'
!
         end if
!
		case(1)	!A shall be TR
!*
!
		do ib=1,Mol_Set(MomB)%N_atoms
        	do i=1,3; xyztb(ib,i)=Optxyz(MomB,ib,i); end do
        	END do
!
!		Before anything try to find out simple optimal positioning
!
		Sobmax=zero
!
!		initialize spare coordinates
!
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 ; xyztaa(j,itran)=Mol_Set(MolA)%XYZ(j,itran)
        end do
        end do
!
        IF(Azimut)then
!
	Do INX=0,7
!
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 
           If (IXYZA(INX,ITRAN)==0)then
          xyzta(j,itran)=Mol_Set(MolA)%XYZ(j,itran)
          else
            xyzta(j,itran)=-Mol_Set(MolA)%XYZ(j,itran)
            end if
           end do
        end do
!
!	compute position overlap similarity element with simply rotated coordinates
!
		Call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!    
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM
!
!		Put spare coordinates into maximal positioning
!
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 ; xyztaa(j,itran)=xyzta(j,itran)
        end do
        end do
!
        END IF
!
        END DO  !INX
        END IF  !AZIMUT
!
	ibtb=0
        do iata=1,Mol_Set(MolA)%N_atoms; naz=Mol_Set(MolA)%Z(iata)
!
        IF(naz>1)then							! Just when an atom is heavy
!
!       move molecule  A put atom  Iata at origin
!
	if(LLORO)then
        WRITE(777,*)'---------------------------------------------------'
        WRITE(777,*)' translation of atom B: ',ibtb,' over atom A: ',iata
        WRITE(777,*)'---------------------------------------------------'
        end if	!LLORO
!
!		START TRANSLATING-ROTATING MOLECULAR FRAMES For Similarity Integral Optimal Search
!
!	    form the global orientation vectors on A and B
!		          use non-hydrogen atoms only
!
        do itran=1,3 ; tran(itran)=xyztaa(iata,itran);end do   
        do j=1,Mol_Set(MolA)%N_atoms
        do itran=1,3 ; xyzta(j,itran)=xyztaa(j,itran)-tran(itran); end do
        end do
!
!		compute alignment coordinates sum vector using heavy atoms only
!
        xa=zero                    
        do ia=1,Mol_Set(MolA)%N_atoms; nza=Mol_Set(MolA)%Z(ia)
        IF(nza>1)then
        do i=1,3; xa(i)=xa(i)+xyzta(ia,i); end do
        END if
        END do
!
!		search for rotations to align sum vectors to X axis
!
		call zalign(xa,ra,caa,saa,cza,sza)
!
!		rotate ALL atomic coordinates of molecule A
!
        do ia=1,Mol_Set(MolA)%N_atoms
        	do i=1,3; xa(i)=xyzta(ia,i); end do 
				call malign(xa,ra,caa,saa,cza,sza)
        	do i=1,3; xyzta(ia,i)=ra(i); end do 
        END do
!
		call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM; natmax=iata
!
!		Store Optimal A TR Coordinates if not done yet
!
        do j=1,Mol_Set(MolA)%N_atoms; do k=1,3; Optxyz(MolA,j,k)=xyzta(j,k);END DO; END do
!        
        END IF !SOVLPM>Sobmax
!
        if(LLORO)then
				WRITE(777,*)'---------------------------------------------------'
				write(777,*)'----- Translated-Rotated Molecules-----------------'
       	        WRITE(777,*)' Overlap               : ',SOVLPM
                Scaled=SOVLPM/Scaling
		WRITE(777,*)'        Scaled by NA*NB: ',Scaled
                IF(SOVLPM>Sobmax)THEN;
                WRITE(777,*)'----- Succesful Rotation-Translation Maximization!-'
                write(777,*)'----- New Overlap Value Used in Similarity Matrix--'
                WRITE(777,*)'---------------------------------------------------'
		END IF
	END if !LLORO
!
	End if  !Natoms A
!
        END do  !iata
!
!       This is to do not optimize anymore A and B coordinates 
!			in case they are allready optimallly aligned
!
        Moldone(MolA)=.True.    
!
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------*OPTIMAL RESULTS*-----------'
        WRITE(777,*)'---------------------------------------'
!
		if(LLORO)then
        WRITE(777,*)'-------for Molecules: '
        write(777,'(1x,"-->",i4,"-*-",a80)')MolA,Mol_Set(MolA)%Identity
        write(777,'(1x,"-->",i4,"-*-",a80)')MomB,Mol_Set(MomB)%Identity
        WRITE(777,*)'---------------------------------------'
		END IF	!LLORO
!
       IF(natmax/=0)then
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)'MAXIMUM Position Overlap Similarity (Translated-Rotated)'
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)' Atoms A - B Displacement    : ',natmax
        WRITE(777,*)' Overlap                     : ',Sobmax
        IF(LLORO)THEN
        SSmx=Sobmax/Scaling
        WRITE(777,*)' Overlap  scaled             : ',SSmx
        rpm=Sobmax/SQRT(Mol_Set(MolA)%overlap_QSS*Mol_Set(MomB)%overlap_QSS)
        WRITE(777,*)' Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(Mol_Set(MolA)%overlap_QSS+Mol_Set(MomB)%overlap_QSS-2*Sobmax)
        WRITE(777,*)' Overlap Distance Index A-B  : ',distpm
        WRITE(777,*)'---------------------------------------'
        END if   !LLORO
!
!	Form the similarity matrices and preserve optimal TR coordinates 
!				for every molecular pair
!
    	ZZ(molmom)=Sobmax
    	CRB_INDX(molmom)=rpm
!
        else
        WRITE(777,*)'++++++++++> NO OPTIMAL VALUE FOUND!'
!
            end if
!*
        case(2)	!B Shall be TR
!**
!		Put the allready TR optimal coordinates in as they are computed beforehand
!
		do ia=1,Mol_Set(MolA)%N_atoms
        	do i=1,3; xyzta(ia,i)=Optxyz(MolA,ia,i); end do
        	END do
!
!		Before anything try to find out simple optimal positioning
!
		Sobmax=zero
!
!		initialize spare coordinates
!
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 ; xyztbb(j,itran)=Mol_Set(MomB)%XYZ(j,itran)
        end do
        end do
!
        IF(Azimut)THEN
        Do JNX=0,7
!
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 
           If (IXYZB(JNX,ITRAN)==0)then
             xyztb(j,itran)=Mol_Set(MomB)%XYZ(j,itran)
             else
              xyztb(j,itran)=-Mol_Set(MomB)%XYZ(j,itran) 
              end if
              end do
        end do
!
!	compute position overlap similarity element with simply rotated coordinates
!
		Call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!    
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM
!
!		Put spare coordinates into maximal positioning
!
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 ; xyztbb(j,itran)=xyztb(j,itran)
        end do
        end do
!        
        END IF
!
	end do  !JNX
        END if	!Azimut
!
	iata=0
        do ibtb=1,Mol_Set(MomB)%N_atoms; nbz=Mol_Set(MomB)%Z(ibtb)
!
        IF(nbz>1)then							! Just when atoms are heavy
!
!       move molecule B over A put atom Ibtb on Iata
!
		if(LLORO)then
        WRITE(777,*)'---------------------------------------------------'
        WRITE(777,*)' translation of atom B: ',ibtb,' over atom A: ',iata
        WRITE(777,*)'---------------------------------------------------'
        end if	!LLORO
!
!		START TRANSLATING-ROTATING MOLECULAR FRAME For Similarity Integral Optimal Search
!		start TR molecule B
!
        do itran=1,3 ; tran(itran)=xyztbb(ibtb,itran);end do
        do j=1,Mol_Set(MomB)%N_atoms
        do itran=1,3 ; xyztb(j,itran)=xyztbb(j,itran)-tran(itran); end do
        end do
!
        xb=zero
        do ib=1,Mol_Set(MomB)%N_atoms; nzb=Mol_Set(MomB)%Z(ib)
        IF(nzb>1)then
        do i=1,3; xb(i)=xb(i)+xyztb(ib,i); end do
        END if
        END do
!
!		search for rotations to align sum vectors to X axis
!
		call zalign(xb,rb,cab,sab,czb,szb)
!
        do ib=1,Mol_Set(MomB)%N_atoms
        	do i=1,3; xb(i)=xyztb(ib,i); end do 
 				call malign(xb,rb,cab,sab,czb,szb)
        	do i=1,3; xyztb(ib,i)=rb(i); end do 
        END do
!
		call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM; nbtmax=ibtb
!
        IF(.not.moldone(MomB))THEN;   !Store Optimal B TR Coordinates if not done yet
        do j=1,Mol_Set(MomB)%N_atoms; do k=1,3; Optxyz(MomB,j,k)=xyztb(j,k);END DO; END do
        END if
!        
        END IF !SOVLPM>Sobmax
!
        if(LLORO)then
		WRITE(777,*)'---------------------------------------------------'
		write(777,*)'----- Translated-Rotated Molecules-----------------'
       	        WRITE(777,*)' Overlap               : ',SOVLPM
                Scaled=SOVLPM/Scaling
		WRITE(777,*)'        Scaled by NA*NB: ',Scaled
                IF(SOVLPM>Sobmax)THEN;
                WRITE(777,*)'----- Succesful Rotation-Translation Maximization!-'
                write(777,*)'----- New Overlap Value Used in Similarity Matrix--'
                WRITE(777,*)'---------------------------------------------------'
		END IF
	END if !LLORO
!
	End if  !Natoms B

        END do  !ibtb
!
!       This is to do not optimize anymore A and B coordinates 
!			in case they are allready optimallly aligned
!
        Moldone(MomB)=.True.     
!
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------*OPTIMAL RESULTS*-----------'
        WRITE(777,*)'---------------------------------------'
!
		if(LLORO)then
        WRITE(777,*)'-------for Molecules: '
        write(777,'(1x,"-->",i4,"-*-",a80)')MolA,Mol_Set(MolA)%Identity
        write(777,'(1x,"-->",i4,"-*-",a80)')MomB,Mol_Set(MomB)%Identity
        WRITE(777,*)'---------------------------------------'
		END IF	!LLORO
!
       IF(nbtmax/=0)then
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)'MAXIMUM Position Overlap Similarity (Translated-Rotated)'
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)' Atoms A - B Displacement    : ',nbtmax
        WRITE(777,*)' Overlap                     : ',Sobmax
        IF(LLORO)THEN
        SSmx=Sobmax/Scaling
        WRITE(777,*)' Overlap  scaled             : ',SSmx
        rpm=Sobmax/SQRT(Mol_Set(MolA)%overlap_QSS*Mol_Set(MomB)%overlap_QSS)
        WRITE(777,*)' Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(Mol_Set(MolA)%overlap_QSS+Mol_Set(MomB)%overlap_QSS-2*Sobmax)
        WRITE(777,*)' Overlap Distance Index A-B  : ',distpm
        WRITE(777,*)'---------------------------------------'
        END IF 	!LLORO
!
!	Form the similarity matrices and preserve optimal TR coordinates 
!				for every molecular pair
!
    	ZZ(molmom)=Sobmax
    	CRB_INDX(molmom)=rpm
!
        else
        WRITE(777,*)'++++++++++> NO OPTIMAL VALUE FOUND!'
!
        end if	!Sobmax
!**
        case(3) !None shall be TR
!***
!		Put the allready TR optimal coordinates in as they are computed beforehand
!
		do ia=1,Mol_Set(MolA)%N_atoms
        	do i=1,3; xyzta(ia,i)=Optxyz(MolA,ia,i); end do
        	END do
		do ib=1,Mol_Set(MomB)%N_atoms
        	do i=1,3; xyztb(ib,i)=Optxyz(MomB,ib,i); end do
        	END do
!
		call ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!
!		No Maximul Value Shall be found
!
        Sobmax=SOVLPM; 
!
        if(LLORO)then
		WRITE(777,*)'---------------------------------------------------'
		write(777,*)'----- Translated-Rotated Molecules-----------------'
       	        WRITE(777,*)' Overlap               : ',SOVLPM
                Scaled=SOVLPM/Scaling
		WRITE(777,*)'        Scaled by NA*NB: ',Scaled
                END if !LLORO
!
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------*RESULTS*-----------'
        WRITE(777,*)'---------------------------------------'
!
		if(LLORO)then
        WRITE(777,*)'-------for Molecules: '
        write(777,'(1x,"-->",i4,"-*-",a80)')MolA,Mol_Set(MolA)%Identity
        write(777,'(1x,"-->",i4,"-*-",a80)')MomB,Mol_Set(MomB)%Identity
        WRITE(777,*)'---------------------------------------'
		END IF	!LLORO
!
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)' Overlap Similarity '
        WRITE(777,*)'----------------------------------------'
        WRITE(777,*)' Overlap                     : ',Sobmax
        rpm=Sobmax/SQRT(Mol_Set(MolA)%overlap_QSS*Mol_Set(MomB)%overlap_QSS)
        IF(LLORO)THEN
        SSmx=Sobmax/Scaling
        WRITE(777,*)' Overlap  scaled             : ',SSmx
        WRITE(777,*)' Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(Mol_Set(MolA)%overlap_QSS+Mol_Set(MomB)%overlap_QSS-2*Sobmax)
        WRITE(777,*)' Overlap Distance Index A-B  : ',distpm
        WRITE(777,*)'---------------------------------------'
        END IF  !LLORO
!
!	Form the similarity matrices and preserve optimal TR coordinates 
!				for every molecular pair
!
    	ZZ(molmom)=Sobmax
    	CRB_INDX(molmom)=rpm
!***
        end select
!
        END do  !Mom
        end do  !Mol
!
!<*********************************************
!
!       Form sum of Outdiagonal Elements of the Similarity Matrix
!
   	MM=(Molom*Molom+Molom)/2
        SumZ=zero; Do III=1,MM; SumZ=SumZ+ZZ(iii); END do
        SumZ=2*(SumZ-TraceZ)
!
        if(SumZMX<SumZ)then
          SumZMX=SumZ
          MXCRIN=NICP
!
    CR_QSMatrix_OUT= TRIM(ID_Label_Pre)//TRIM(CRICRI)//TRIM(QSMatrix_OUT)//TRIM(ID_Label_Suf)//'.txt'
	Open(111,File=CR_QSMatrix_OUT)
!
        WRITE(111,*)'-------------------------------------'
        WRITE(111,*)'Quantum Similarity Matrix'
        WRITE(111,*)'-------------------------------------'
        WRITE(111,'(a80)')General_title
        WRITE(111,*)MoloM,MM
        WRITE(111,*)'------------------------------------'
        if(KPindx==1)then
        WRITE(111,*)' Circular Permutation:',NICP
          else if(KPindx==2) then
        WRITE(111,*)' Random Permutation:',NICP
          end if
        WRITE(111,*)' Sum of Z Triangles: ',SumZ
        WRITE(111,*)(KTOrder(I),I=1,MoloM)
        write(111,*)(iii,ZZ(iii),iii=1,MM)
        WRITE(111,*)'-------------------------------------'
        CLOSE(111)
!
        CR_CRBINDX_OUT= TRIM(ID_Label_Pre)//TRIM(CRICRI)//TRIM( CRBINDX_OUT)//TRIM(ID_Label_Suf)//'.txt'
        OPEN(222,FILE= CR_CRBINDX_OUT)
!
        WRITE(222,*)'-------------------------------------'
        WRITE(222,*)' Carbó Index Matrix'
        WRITE(222,*)'-------------------------------------'
        WRITE(222,'(a80)')General_title
        WRITE(222,*)MoloM,MM
        WRITE(222,*)'-------------------------------------'
       	if(KPindx==1)then
          WRITE(222,*)' Circular Permutation:',NICP
          else if(KPindx==2) then
          WRITE(222,*)' Random Permutation:',NICP
          end if
        SumZ=zero; Do III=1,MM; SumZ=SumZ+CRB_INDX(iii); END do
        SumZ=2*(SumZ-MoloM)
        WRITE(222,*)' Sum of Carbó Index Triangles: ',SumZ
        WRITE(222,*)(KTOrder(I),I=1,MoloM)
 	WRITE(222,*)(iii,CRB_INDX(iii),iii=1,MM)
        WRITE(222,*)'-------------------------------------'
        CLOSE(222)
!
	if(Zenith)then
!        
        CR_OPT_XYZ_OUT= TRIM(ID_Label_Pre)//TRIM(CRICRI)//TRIM(OPT_XYZ_OUT)//TRIM(ID_Label_Suf)//'.txt'
	OPEN(333,FILE= CR_OPT_XYZ_OUT)
!                               
         WRITE(333,*)'-------------------------------------'
         WRITE(333,*)' Optimal FA Coordinates'
         WRITE(333,*)'-------------------------------------'
         WRITE(333,'(a80)')General_title
         WRITE(333,*)MoloM
         WRITE(333,*)'-------------------------------------'
       	 if(KPindx==1)then
          WRITE(333,*)' Circular Permutation:',NICP
          else if(KPindx==2) then
          WRITE(333,*)' Random Permutation:',NICP
          end if !KPindex
       	WRITE(333,*)(KTOrder(I),I=1,MoloM)
        do Mol=1,MoloM
        MolA=KTOrder(Mol)
        Nmol=Mol_Set(MolA)%N_atoms
        WRITE(333,*)'---->Molecule : ',Mol
        WRITE(333,'(1x,a80)')Mol_Set(MolA)%Identity
        WRITE(333,*)Nmol
        do j=1,Nmol
        WRITE(333,*)(Mol_Set(Mola)%Z(j),Optxyz(Mola,j,k),k=1,3)
        END do
        END do
        WRITE(333,*)'-------------------------------------'
        CLOSE(333)
        end if	!Zenith
!
        END if  !SUMZ
!
        END do 	!NICP
!
	WRITE(777,*)'****************************************'
        WRITE(777,*)'---------------------------------------'
	write(777,*)'Optimal CP Superposition  : ',MXCRIN
	write(777,*)'Maximal CP Z triangle SUM : ',SumZMX
        WRITE(777,*)'---------------------------------------'
!
        CLOSE(777)
!
!--------------------------------------------------------------------------------------------------------------------
    WRITE(*,*)' Molecular Quantum Similarity Suite'
    write(*,*)'Circular or Random Permutations Success!'
    pause
    stop
!	
	end                !program
!...............................................................................................................
!..................................................................................................
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
	subroutine zalign(a,ra,ca,sa,cz,sz)
!    
!	Aligns any 3D vector with the X axis:
!	1) rotates clockwise around Z axis to put it into the ZX plane
!		using the cosine ca and the sine sa
!   2) rotates clockwise around the new Y axis to put it into the X axis
!		using the cosine cz and the sime sz.
!
	implicit double precision(a-h,o-z) 
!
	parameter(zero=0.0d0,toler=1.0d-4,one=1.0d0)
!      
	DIMENSION a(3),a2(3),ra(3)
!
    a2=a*a
!
!	compute rotation gear for vector a
!
	pam=a2(1)+a2(2); sqpam=sqrt(pam)	
	if(sqpam>=toler)then	
		ca=a(1)/sqpam;sa=a(2)/sqpam
    else
!      
!	the vector is already in the Z axis
!
        ca=one;sa=zero
	end if
	sqam=sqrt(pam+a2(3))
    if(sqam>=toler)then
    cz=sqpam/sqam;sz=a(3)/sqam
    else
!      
!	the vector is already in the Z axis
!
        cz=one;sz=zero
	end if
!
!	rotate A
!  
    ax=ca*a(1)+sa*a(2)
    ra(2)=-sa*a(1)+ca*a(2)
	ra(1)=cz*ax+sz*a(3)
    ra(3)=-sz*ax+cz*a(3)
!
	end
!-----------------------------------------------------------------------------------        
	subroutine malign(a,ra,ca,sa,cz,sz)
!
!	just rotates along Z axis with {ca,sa}  
!   and along new Y axis using {cz,sz}
!
	implicit double precision(a-h,o-z) 
!
	DIMENSION a(3),ra(3)
!
!	rotate vector a
!  
    ax=ca*a(1)+sa*a(2)
    ra(2)=-sa*a(1)+ca*a(2)
	ra(1)=cz*ax+sz*a(3)
    ra(3)=-sz*ax+cz*a(3)
!
	end
!--------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
!--------------------------------------------------------------------------
	Subroutine ASA_Overlap(Mol_Set,MolA,MomB,xyzta,xyztb,SOVLPM)
!
!	compute position overlap similarity element with translated-rotated coordinates
!
    use set_molecule
!
    parameter(MOLMAX=160)  			!Maximum number of Molecules
    DOUBLE PRECISION:: zero 
    PARAMETER (zero=0.0d0)
    double precision:: toler,PI
    PARAMETER (toler=1.0D-7)
    Parameter (PI=3.14159265358979D0)
!
    double precision :: xyzta(NatMX,3),xyztb(NatMX,3)
    double precision :: SOVLPM
    double precision :: RAB,term,AA,BB,alam,Fecit,SIMOVLP
!
	Type(Molecule)Mol_Set(MOLMAX)
!
		SOVLPM=zero
        do i=1,Mol_Set(MolA)%N_atoms
        nfi=Mol_Set(MolA)%NASAfun(i)
        do j=1,Mol_Set(MomB)%N_atoms
        nfj=Mol_Set(MomB)%NASAfun(j)
        RAB=zero
        do kik=1,3;RAB=RAB+(xyzta(i,kik)-xyztb(j,kik))**2;end do
        term=zero
        do k=1,nfi
        AA=Mol_Set(MolA)%Ep(i,k)
        do l=1,nfj
        BB=Mol_Set(MomB)%Ep(j,l)
!
!        overlap similarity between s-type functions
!
        alam=2*AA*BB/(AA+BB); Fecit=SQRT(alam/PI)**3
        IF(ABS(RAB)<=toler)then
        SIMOVLP=Fecit
        else
        SIMOVLP=Fecit*EXP(-alam*RAB)
        END if
        term=term+Mol_Set(MolA)%Op(i,k)*Mol_Set(MomB)%Op(j,l)*SIMOVLP
!
        end do
        end do
        SOVLPM=SOVLPM+term
        end do
        end do
!
		end
!--------------------------------------------------------------------------
!
!--------------------------------------------------------------------------
!-!