!     Last change:  RCD  24 Jun 2010   12:12 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*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 build_an_ASA_molecule_with_XYZ_Quadrupole_orientation
	use set_molecule
!
    parameter(MOLMAX=70)	!Maximum number of Molecules
    PARAMETER(MXPERAT=60)	!Maximum number of ASA different atoms
!
    DOUBLE PRECISION:: toau, zero
    PARAMETER (toau=0.529177249D0)
    PARAMETER (zero=0.0d0, one=1.0d0)
    double precision  toler
    PARAMETER (toler=1.0D-7)
!   
    double precision SSP
    double precision TupaK(3),TopaK(3)
!             
    CHARACTER (LEN=80):: GENERAL_TITLE
    CHARACTER (LEN=80):: Part_title(Molmax),Particulartitle
    CHARACTER (LEN=12):: ASAFUN		!Name of the file for ASA function type
!
        INTEGER:: Nuf(MXPERAT)
        DOUBLE PRECISION :: App(MXPERAT,maxfn),Cpp(MXPERAT,maxfn)
        double precision :: szp						! number of electrons
   		double precision :: ssip,ccp,Testp,tistp,ZZ	! auxiliary parameters to test ASA coefficients coherence
		double precision :: R(3),RR(3),CC,SS,CCZ,SSZ!auxiliary to rotate coordinates
        double precision :: Cent(3) , DumDum
        Logical:: LLORO								! .true. write more info
        Logical:: ThreeD_Change, Order_XYZ, Dipole_align, Positive_Octant, Quadrupole_align,Center_of_charges
!
	Type(Molecule)Mol_Set(MOLMAX)
! 
	integer MoloM, reading
!------------------------------------------------------------------------------------------------------------    
!	This is to allow Spartan output direct reading
!
    Character*2 zip
    Character*7 zup
    Parameter (MaxZan=54)
    Character*2 PTaSYm(1:MaxZan)
    Data PTaSYm/&
    'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na','Mg','Al','Si','P ','S ','Cl','Ar',&
    'K','Ca','Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr',&
    'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe'/
    INTEGER MopacInner(1:MaxZAn)
    DATA MopacInner/0,0,2,2,2,2,2,2,2,2,10,10,10,10,10,10,10,10,&
    18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,&
    36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36/
!-----------------------------------------------------------------------------------------------------------
	integer IXYZ(3),iG,iP,iQ,ixyzi
    Double precision :: Dipole(3)
	Double Precision :: QuadM(6),QPC(3,3),QEig(3),Durius
!------------------------------------------------------------------------------------------------------------
        Character*70 :: ASA_Basis_Set, Coordinate_units, Pop_Char, Pop_Char_type, Reading_Format
        Character*70 :: ID_Label_Pre,ID_Label_Suf
        Character*70 :: OUTPUT_FIL,INPUT_FIL_XYZ,INPUT_FIL_POP,PREPARED_FIL
        Character*120 :: MOLECINXYZ,MOLECINPOP,MOLECOUT,MOLECPRE,MOLECOUTXYZ
!------------------------------------------------------------------------------------------------------------
!    
!	Prepares a batch of ASA defined molecules
!
!	First reads verbose and transformation flags
!	all flags are .false. at the beginning
!
	 Lloro=.false.; ThreeD_Change=.false.
	 Order_XYZ=.false.;  Center_of_charges=.false.; Positive_Octant=.false.
     Dipole_align=.false.; Quadrupole_align=.false.
!
     open(403,File='Flagstaff_1.txt')
	 read(403,*)ixyzi;	if(ixyzi==1)LLORO=.true.
     read(403,*)ixyzi;	if(ixyzi==1)ThreeD_Change=.true.
     if(ThreeD_Change)then
     read(403,*)ixyzi;	if(ixyzi==1) Order_XYZ=.true.
     read(403,*)ixyzi;	if(ixyzi==1) Center_of_charges=.true.
     read(403,*)ixyzi;	if(ixyzi==1) Positive_Octant=.true.
     read(403,*)ixyzi;	if(ixyzi==1) Dipole_align=.true.
     read(403,*)ixyzi;	if(ixyzi==1) Quadrupole_align=.true.
     end if
	 close(403)
!
!		second reads a label file, where it is found
!		information on ASA basis set, coordinates, Populations or charges
!		Type of Populations or charges, reading format
!
		open(304,File='Labelstaff_1.txt')
        read(304,*)ASA_Basis_Set
        if(Trim(ASA_Basis_Set)=='3-21G')then
          Kind_ASA=1
        else if(Trim(ASA_Basis_Set)=='6-311G')then
          Kind_ASA=2
        else if(Trim(ASA_Basis_Set)=='Huzinaga')then  
          Kind_ASA=3
          end if
!
        ANG_AUS=0  
        read(304,*)Coordinate_units
        if(Trim(Coordinate_Units)=='Angstroms')ANG_AUS=1
!          
        read(304,*)Pop_Char
		if(TRIM(Pop_Char)=='Populations')then
		nat_char=0
 		else if(TRIM(Pop_Char)=='Charges')then
		nat_char=1
          end if
!        
        read(304,*)Pop_Char_type
        if(TRIM(Pop_Char_type)=='Promolecular')then
         NAP_TYP=0
         else if(TRIM(Pop_Char_type)=='Mulliken')then 
	    NAP_TYP=1
            else if(TRIM(Pop_Char_type)=='NO')then
              NAP_TYP=2
              else IF(TRIM(Pop_Char_type)=='ASABAD')then
                NAP_TYP=3
              end if
!
        read(304,*)Reading_Format
		if(TRIM(Reading_Format)=='Normal')then
		reading=0
        else if(TRIM(Reading_Format)=='Spartan')then
		reading=1
        else IF(TRIM(Reading_Format)=='Mopac')then
                reading=2
        else
                WRITE(*,*)' Error in XYZ file: Format of XYZ not supported STOP'
                pause
                Stop
          end if
!
	Read(304,'(a70)')ID_Label_Pre
        Read(304,'(a70)')ID_Label_Suf
        Read(304,'(a70)')INPUT_FIL_XYZ
        Read(304,'(a70)')INPUT_FIL_POP
        Read(304,'(a70)')OUTPUT_FIL
        Read(304,'(a70)')PREPARED_FIL
!
        close(304)
!
	MOLECINXYZ=TRIM(ID_Label_Pre)//TRIM(INPUT_FIL_XYZ)//TRIM(ID_Label_Suf)//'.txt'
	open(666,file=MOLECINXYZ)
!
        MOLECOUT=TRIM(ID_Label_Pre)//TRIM(OUTPUT_FIL)//TRIM(ID_Label_Suf)//'.txt'
        open(777,file=MOLECOUT)
        MOLECPRE=TRIM(ID_Label_Pre)//TRIM(PREPARED_FIL)//TRIM(ID_Label_Suf)//'.txt'
!
!		prepare the molecular set for ASA calculations
!
!       read a generic title
!
        WRITE(777,*)'   '
        WRITE(777,*)'*********************************************************'
        WRITE(777,*)'   '
!
        READ(666,fmt='(A80)')GENERAL_TITLE
        READ(666,*)Molom
        WRITE(777,fmt='(A80)')GENERAL_TITLE
!
!       read number of molecules
!
        IF(Molom>Molmax)then
        WRITE(777,*)' * ERROR 0: Number of Molecules exceeds: ',Molmax,' *'
        WRITE(*,*)' ERROR 0: Number of Molecules exceeds: ',Molmax
        stop
        END if
        WRITE(777,*)' Number of molecules in this run: ',MoloM
!
!       enter the kind of ASA functions you will use
!		KIND_ASA: 1=> 3-21g  2=> 6-31g  3=> Huzinaga
!
        SELECT CASE(KIND_ASA)
        CASE(1)
        WRITE(777,*)' ASA basis set: 3-21G'
        ASAFUN='3-21g___.txt'
        CASE(2)
        WRITE(777,*)' ASA basis set: 6-311G'
        ASAFUN='6-311g__.txt'
        CASE(3)
        WRITE(777,*)' ASA basis set: Huzinaga'
        ASAFUN='Huzinaga.txt'
        case default
        WRITE(777,*)'*  ERROR 1: Specified ASA basis set file undefined *'
        WRITE(*,*)' ERROR 1: Specified ASA basis set file undefined'
        STOP
        END SELECT
!
!       store in CPU memory contents of chosen ASA basis set file
!	these files shall be present in the running program folder
!
        WRITE(777,'(" ASA basis set file used: ",A12)')ASAFUN
        OPEN(111,FILE=ASAFUN)
!
       		 do
!
	   	 read(111,*,END=7)iz,nf		!atomic number, number of Asa functions
          	Nuf(iz)=nf			!Nuf contains the number of ASA functions
          	IF(KIND_ASA<=2)then
	  	read(111,*)(Cpp(iz,i),App(iz,i),i=1,nf)
          	else
	  	read(111,*)(App(iz,i),Cpp(iz,i),i=1,nf)
         	 END if
!
	 	enddo
!
7         CLOSE(111)
!
!       Enter molecules and attach ASA functions to each atom
!
		milmil=(Molom*Molom+Molom)/2
        WRITE(777,*)'***** Number of compact QSM elements :  ',milmil
!
	if(NAP_TYP==1.or.NAP_TYP==2)then
   	MOLECINPOP=TRIM(ID_Label_Pre)//TRIM(INPUT_FIL_POP)//TRIM(ID_Label_Suf)//'.txt'
        open(636,File=MOLECINPOP)
        READ(636,fmt='(A80)')GENERAL_TITLE
        READ(636,*)Nadiuska
        	If(Nadiuska/=MoloM)then
       		 write(*,*)' ERROR: total Number of Molecules not coincident in XYZ and POP input files!'
        	pause
       		 Stop
        	end if
        END if ! NAP_TYP=1 or 2
!
    	do Mol=1,MoloM
        write(*,*)Mol
        WRITE(777,*)'   '
        WRITE(777,*)'*********************************************************'
        WRITE(777,*)'   '
       	READ(666,fmt='(A80)')Part_title(Mol)
        Mol_Set(Mol)%Identity=Part_title(Mol)
        write(777,fmt='(A80)')Part_title(Mol)
        READ(666,*)Mol_Set(Mol)%N_Atoms         !,Mol_Set(Mol)%ANG_AU
        Mol_Set(Mol)%ANG_AU=ANG_AUS
        write(777,*)'Number of atoms: ',Mol_Set(Mol)%N_Atoms
        IF(Mol_Set(Mol)%ANG_AU/=0)WRITE(777,*)' Coordinates transformed to AU'
!
	 	do ia=1,Mol_Set(Mol)%N_Atoms
!
!               normal reading
!
	 	if(reading==0)then
		READ(666,*)Mol_Set(Mol)%Z(ia),(Mol_Set(Mol)%XYZ(ia,j),j=1,3)
!
        else if(reading==1) then
!
!		direct Spartan reading
!
		READ(666,*)iiaa,zip,zup,(Mol_Set(Mol)%XYZ(ia,j),j=1,3)
        	Do mist=1,MaxZan
          	if(zip==PTaSYm(mist))exit
            	end do
            if(mist>=1.and.mist<=MaxZan)then
         	Mol_Set(Mol)%Z(ia)=mist
      		else
            Write(*,*)' Error in Atom assignation:'
            Write(*,*)' Atomic Number not valid or out of allowed range!',mist
            Write(*,*)' Maximum Atomic Number is: ',MaxZan
            pause
            Stop
          	end if
!
        else if(reading==2) then
!
!		direct Mopac reading
!
		READ(666,*)iiaa,zip,(Mol_Set(Mol)%XYZ(ia,j),j=1,3)
        	Do mist=1,MaxZan
          	if(zip==PTaSYm(mist))exit
            	end do
            if(mist>=1.and.mist<=MaxZan)then
         	Mol_Set(Mol)%Z(ia)=mist
!
         	else
            Write(*,*)' Error in Atom assignation:'
            Write(*,*)' Atomic Number not valid or out of allowed range!',mist
            Write(*,*)' Maximum Atomic Number is: ',MaxZan
            pause
            Stop
         	 end if
!
          end if
!
        IF(Mol_Set(Mol)%ANG_AU/=0)THEN
           do j=1,3
              Mol_Set(Mol)%xyz(ia,j)=Mol_Set(Mol)%xyz(ia,j)/toau
           end DO
        END if
!
        END do    !ia : # atoms in molecule Mol
!
!       read according NAP_TYPE the indicative of Atomic Population type 
!		and nat_char which if zero corresponds to populations
!		but if one corresponds to charges
!		the populations or charges one by one
!
!		if populations or charges are read then open the corresponding file and check
!
	if(NAP_TYP==1.or.NAP_TYP==2)then
        READ(636,fmt='(A80)')Particulartitle
        read(636,*)Nadiuskilla
        if(Nadiuskilla/=Mol_Set(Mol)%N_Atoms)then
          write(*,*)' Error!'
          WRITE(*,*)' For molecule . ',Mol
          write(*,*)' Number of atoms in XYZ and POP files not coincident!'
          pause
          stop
          end if
          end if !NAP_TYP = 1 or 2
!
	Mol_Set(Mol)%NAP_TYPE=NAP_TYP
         select case (Mol_Set(Mol)%NAP_TYPE)
          case (0)			! Atomic Numbers
          do i=1,Mol_Set(Mol)%N_Atoms
             IF(reading<=1)Mol_Set(Mol)%AP(i)=Mol_Set(Mol)%Z(i)
!
!       if Promolecular and Mopac Labels on use only valence electrons
!
             IF(reading==2)Mol_Set(Mol)%AP(i)=Mol_Set(Mol)%Z(i)-MopacInner(Mol_Set(Mol)%Z(i))
!
              END DO
          WRITE(777,*)' Promolecular DF'
          case (1)			! Mulliken Atomic Gross Populations
          do i=1,Mol_Set(Mol)%N_Atoms
             if(reading==0)THEN; READ(636,*)Mol_Set(Mol)%AP(i)
               else if(reading==1)THEN; READ(636,*)Zup,Mol_Set(Mol)%AP(i)
                   else IF(reading==2)THEN;
                             IF(nat_char==0)READ(636,*)Zup,Zup,DumDum,Mol_Set(Mol)%AP(i)
              		     IF(nat_char==1)READ(636,*)Zup,Zup,Mol_Set(Mol)%AP(i)
                 end if
              END DO
          WRITE(777,*)' Promolecular Mulliken Populations DF'
!
          case (2)			! Natural Orbitals Atomic Populations
          do i=1,Mol_Set(Mol)%N_Atoms
             if(reading==0)then; READ(636,*)Mol_Set(Mol)%AP(i)
               else if(reading==1)then; READ(636,*)Zup,Mol_Set(Mol)%AP(i)
               else IF(reading==2)THEN;
               IF(nat_char==0)READ(636,*)Zup,Zup,DumDum,Mol_Set(Mol)%AP(i)
               IF(nat_char==1)READ(636,*)Zup,Zup,Mol_Set(Mol)%AP(i)
               end if
              END DO
          WRITE(777,*)' Promolecular NO Populations DF'
!
          CASE(3)
!
!         ASA Bulk Atomic Density: All ASA coefficients are chosen unit
!
          do i=1,Mol_set(Mol)%N_Atoms
          Mol_Set(Mol)%AP(i)=one
          end do
          WRITE(777,*)' Promolecular ASA Bulk Atomic Density (BAD)'
          case default
          WRITE(777,*)' ERROR 2: Specified Atomic Population type Inexistent *'
          WRITE(*,*)' ERROR 2: Specified Atomic Population type Inexistent *'
          STOP
          end select
!
	if(Mol_Set(Mol)%NAP_TYPE/=0.and.nat_char==1)then
        write(777,*)'  Atomic Charges transformed to Populations'
        do i=1,Mol_Set(Mol)%N_Atoms
        IF(reading==1)Mol_Set(Mol)%AP(i)=Mol_Set(Mol)%Z(i)-Mol_Set(Mol)%AP(i)
        IF(reading==2)Mol_Set(Mol)%AP(i)=Mol_Set(Mol)%Z(i)-(Mol_Set(Mol)%AP(i)+MopacInner(Mol_Set(Mol)%Z(i)))
        end do
        end if
!
!       Build up and store ASA exponents and coefficients for each atom in molecule Mol
!       This shall be done every time molecule Mol ASA functions shall be used anywhere
!       Note that exponents A and Coeficients C are indexed as (atom,function)
!
        do i=1,Mol_Set(Mol)%N_Atoms
        natz=Mol_Set(Mol)%Z(i)
        nf=Nuf(natz)
        Mol_Set(Mol)%NASAfun(i)=nf
        do j=1,nf
        Mol_Set(Mol)%Ep(i,j)=App(natz,j); Mol_Set(Mol)%Op(i,j)=Cpp(natz,j)
        IF(KIND_ASA<=2)Mol_Set(Mol)%Op(i,j)=Mol_Set(Mol)%Op(i,j)/natz
        end do
        end do
!
!       Test the ASA coefficients for each atom in both molecules
!
        ssip=zero; Testp=zero; szp=zero; ccp=zero
		do j=1,Mol_Set(Mol)%N_Atoms
         tistp=zero
         szp=szp+Mol_Set(Mol)%Z(j)
         ZZ=Mol_Set(Mol)%AP(j)
         ccp=ccp+ZZ
         sip=zero
       	  	do i=1,Mol_Set(Mol)%NASAfun(j)
          sip=sip+Mol_Set(Mol)%Op(j,i)
          Mol_Set(Mol)%Op(j,i)=Mol_Set(Mol)%Op(j,i)*ZZ  !ASA coefficients become multiplyed by the atomic population
          tistp=tistp+Mol_Set(Mol)%Op(j,i)
      	   	enddo
         ssip=ssip+sip
         Testp=Testp+tistp
	   	enddo
		Mol_Set(Mol)%N_electrons=szp
!
        if(LLORO)then
		write(777,*)Mol,Part_title(Mol)			!Molecule number and identity
        write(777,*)Mol_Set(Mol)%N_atoms		!Number of atoms
        write(777,*)Mol_Set(Mol)%N_electrons		!Number of electrons
    	do natl=1,Mol_Set(Mol)%N_atoms
        write(777,*)natl				!atom number
   		write(777,*)Mol_Set(Mol)%Z(natl)		!Nuclear charge
    	write(777,*)Mol_Set(Mol)%AP(natl)		!Gross atomic population
		write(777,*)(Mol_Set(Mol)%XYZ(natl,k),k=1,3)	!Nuclear Cartesian coordinates (in au)
    	na_boo=Mol_Set(Mol)%NASAfun(natl)
        write(777,*)na_boo				!Number of ASA basis functions
    	write(777,*)(Mol_Set(Mol)%Ep(natl,mfn),mfn=1,na_boo)  		!Exponents
    	write(777,*)(Mol_Set(Mol)%Op(natl,mfn),mfn=1,na_boo)		!Coefficients
		end do
		WRITE(777,*)'.......................................................................'
        WRITE(777,*)' Testing data coherence:'
        WRITE(777,*)' Total Molecular number of atoms from raw coefficients:   ',ssip
        WRITE(777,*)' Total Molecular number of electrons from atomic numbers: ',szp
        WRITE(777,*)' Total Molecular number of electrons from charges:        ',ccp
        WRITE(777,*)' Test of Coefficient coherence:                           ',Testp
        WRITE(777,*)'........................................................................'
	end if
!
	end do !MOL
!
	close(666)
    	close(777)
!
    	if(NAP_TYP==1.or.NAP_TYP==2)Close(636)
!
!		write initial coordinates
!
		MOLECINXYZ=TRIM(ID_Label_Pre)//'Original_IN_XYZ'//TRIM(ID_Label_Suf)//'.txt'        
		open(333,File=MOLECINXYZ)
!
		WRITE(333,*)'---------------------------------------------'
       	write(333,*)'-----Original Initial Coordinates------------'
		WRITE(333,*)'---------------------------------------------'
        WRITE(333,fmt='(A80)')GENERAL_TITLE
        Write(333,*)MoloM
!
        do MolA=1,MoloM
        Nmol=Mol_Set(MolA)%N_atoms
        WRITE(333,*)'---->Molecule : ',MolA
        WRITE(333,'(a80)')Part_title(MolA)
        WRITE(333,*)Nmol
        do j=1,Nmol
        WRITE(333,*)Mol_Set(MolA)%Z(j),(Mol_Set(MolA)%XYZ(j,k),k=1,3)
        END do
        END do
        WRITE(333,*)'-------------------------------------'
!        
		close(333) 
!
!		write prepared data
!
		open(232,File=MOLECPRE)
!
        WRITE(232,*)GENERAL_TITLE
        SELECT CASE(KIND_ASA)
        CASE(1)
        WRITE(232,*)' ASA basis set: 3-21G'
        CASE(2)
        WRITE(232,*)' ASA basis set: 6-311G'
        CASE(3)
        WRITE(232,*)' ASA basis set: Huzinaga'
        END SELECT
        Write(232,*)Kind_ASA
        write(232,*)MoloM					!Number of molecules
!
	    do Mol=1,MoloM
! 
!	test the coordinates to unify coordinate directions to lie 
!	linear molecules forced in the X direction
!	planar molecules forced in the XY plane
!
!	compute number of non hydrogen atoms        
!
		nonH=0
		do natl=1,Mol_Set(Mol)%N_atoms
        if(Mol_Set(Mol)%Z(natl)>1)nonH=nonH+1
        end do
!        
		do ixyzi=1,3
        Tupak(ixyzi)=zero
		do natl=1,Mol_Set(Mol)%N_atoms
!        
!	use only non hydrogen atoms to test where the molecule backbone lies          
!
        if(Mol_Set(Mol)%Z(natl)>1)Tupak(ixyzi)=Tupak(ixyzi)+Dabs(Mol_Set(Mol)%XYZ(natl,ixyzi))
        end do
        end do
!
!	order the coordinates according the greater absolute sum value
!
 		if(Order_XYZ)then
		Do iP=1,3
        IXYZ(iP)=iP
        end do
!        
		DO iP=1,3
   		DO iQ=iP,3
      	if (Tupak(iQ)>Tupak(iP)) then
         iG=IXYZ(iP)
         IXYZ(iP)=IXYZ(iQ)
         IXYZ(iQ)=iG
         Durius=Tupak(iP)
         Tupak(iP)=Tupak(iQ)
         Tupak(iQ)=Durius
        endif
   		end do
		end do
!
        do natl=1,Mol_Set(Mol)%N_atoms
        DO iP=1,3
        iG=IXYZ(iP)  
        Cent(iG)=Mol_Set(Mol)%XYZ(natl,iP)
        end do
        do iP=1,3
        Mol_Set(Mol)%XYZ(natl,iP)=Cent(iP)        
		end do
        end do
!
		end if !Order_XYZ        
!
        Nemini=0
        do ixyzi=1,3
        if(Tupak(ixyzi)<=toler)then
          Nemini=Nemini+1
          end if
        end do
!
!		Check coordinate problems and proceed to rearrange if fine 
!
        IF(Nemini<3)then
!
		if(ThreeD_Change)then
!
!		put the molecule into the positive octant as much as possible
!
		if(Positive_Octant)then
		do ixyzi=1,3
        TopaK(ixyzi)=zero
		do natl=1,Mol_Set(Mol)%N_atoms
        Topak(ixyzi)=Topak(ixyzi)+Mol_Set(Mol)%XYZ(natl,ixyzi)
        end do
        end do
!		 
		do J=1,3
        If(Topak(J)<zero)then
          	do natl=1,Mol_Set(Mol)%N_atoms
              Mol_Set(Mol)%XYZ(natl,J)=-Mol_Set(Mol)%XYZ(natl,J)
              end do
          end if
          end do
!
		end if !Positive_Octant          
!
!		compute center of charges and make it the molecular center
!
		if(Center_of_Charges)then
		Cent=zero
		do natl=1,Mol_Set(Mol)%N_atoms
        do J=1,3
        Cent(J)=Cent(J)+Mol_Set(Mol)%XYZ(natl,J)*Mol_Set(Mol)%Z(natl)
        end do
        end do
        do J=1,3
        Cent(j)=Cent(j)/Mol_Set(Mol)%N_electrons
        end do
        do natl=1,Mol_Set(Mol)%N_atoms
        do J=1,3
		Mol_Set(Mol)%XYZ(natl,J)=Mol_Set(Mol)%XYZ(natl,J)-Cent(J)
		end do
        end do
		end if	!Center_of_Charges
!
!		orient coordinates according to  nuclear dipole moment
!
		if(Dipole_align)then
		Dipole=zero
        do natl=1,Mol_Set(Mol)%N_atoms
        Do J=1,3
          Dipole(J)=Dipole(J)+Mol_Set(Mol)%XYZ(natl,J)*Mol_Set(Mol)%Z(natl)
		end do
        end do
!
!		align dipole with X axis
!      
        call zalign(Dipole,RR,CC,SS,CCZ,SSZ)
        Dipole(1)=RR(1)
!
!		rotate all coordinates according to the dipole moment alignment 
!
        do natl=1,Mol_Set(Mol)%N_atoms
        	do J=1,3; R(J)=Mol_Set(Mol)%XYZ(natl,J); end do 
			call malign(R,RR,CC,SS,CCZ,SSZ)
        	do J=1,3; Mol_Set(Mol)%XYZ(natl,J)= RR(J); end do 
        END do	
!
        end if		!Dipole_align
!
!       Orient molecular coordinates to the Nuclear
!		Quadrupole Moment principal components
!
!		Compute Nuclear Quadrupole Moment
!
		if(Quadrupole_align)then
		QuadM=zero
!
		do K=1,3
        KK=(K*K-K)/2
        do J=1,K
          JK=J+KK
		do natl=1,Mol_Set(Mol)%N_atoms
        QuadM(JK)=QuadM(JK)+Mol_Set(Mol)%XYZ(natl,J)*Mol_Set(Mol)%XYZ(natl,K)*Mol_Set(Mol)%Z(natl)
        end do
        end do
        end do       
!										
		call JADIAG(3,6,QuadM,QEig,QPC)		
!
!		Transform the coordinates to the principal components of the Quadrupole
!
		do natl=1,Mol_Set(Mol)%N_atoms
        R=zero
        do J=1,3
          do K=1,3
            R(J)=R(J)+QPC(K,J)*Mol_Set(Mol)%XYZ(Natl,K)
            end do
            end do
        	do L=1,3
          	Mol_Set(Mol)%XYZ(natl,L)=R(L)
          	end do
          end do
		end if		!Quadrupole_align
!
		else
        continue
		end if	!ThreeD_Change               
!
        else
        if(nonH/=1)then
        write(*,*)' Error in backbone coordinates: They are zero!'
        write(232,*)' Error in backbone coordinates: They are zero!'
        write(*,*)' In molecule : ',Mol
        write(232,*)' In molecule : ',Mol
       	write(232,*)' Discontinued Computation!'
        Pause
        Stop
        end if
!
        end if   	!Nemini
!        
        Mol_Set(Mol)%Identity=Part_title(Mol)
		write(232,*)Mol,Mol_Set(Mol)%Identity						!Molecule number and identity
  		write(232,*)Mol_Set(Mol)%N_atoms					    	!Number of atoms
        write(232,*)Mol_Set(Mol)%N_electrons					    !Number of electrons
    	do natl=1,Mol_Set(Mol)%N_atoms	     		 
    	write(232,*)natl											!atom number
        write(232,*)Mol_Set(Mol)%Z(natl)							!Nuclear charge
    	write(232,*)Mol_Set(Mol)%AP(natl)							!Gross atomic population
		write(232,*)(Mol_Set(Mol)%XYZ(natl,k),k=1,3)				!Nuclear Cartesian coordinates (in au)
    	na_boo=Mol_Set(Mol)%NASAfun(natl) 		 
    	write(232,*)na_boo											!Number of ASA basis functions
    	write(232,*)(Mol_Set(Mol)%Ep(natl,mfn),mfn=1,na_boo)  		!Exponents
    	write(232,*)(Mol_Set(Mol)%Op(natl,mfn),mfn=1,na_boo)		!Coefficients
		end do
!    
!		Compute overlap Quantum Selfsimilarity
!
		call ASA_Overlap_Self(Mol_Set,Mol,SSP)
!
		Mol_Set(Mol)%overlap_QSS=SSP
		write(232,*)Mol_Set(Mol)%overlap_QSS						!Position Overlap Quantum Self-Similarity
!
        end do !Mol		
		close(232)
!
		if(ThreeD_Change)then
!
!		write final TR coordinates
!
		MOLECOUTXYZ=TRIM(ID_Label_Pre)//'Original_OUT_XYZ'//TRIM(ID_Label_Suf)//'.txt'         
		open(333,File=MOLECOUTXYZ)
!
		WRITE(333,*)'-----------------------------------------------------'
       	write(333,*)'-----Original OUTPUT Coordinates Modified------------'
		WRITE(333,*)'-----------------------------------------------------'
        WRITE(333,fmt='(A80)')GENERAL_TITLE
        write(333,*)MoloM
!
        do MolA=1,MoloM
        Nmol=Mol_Set(MolA)%N_atoms
        WRITE(333,*)'---->Molecule : ',MolA
        WRITE(333,'(a80)')Mol_Set(MolA)%Identity
        WRITE(333,*)Nmol
        do j=1,Nmol
        WRITE(333,*)Mol_Set(MolA)%Z(j),(Mol_Set(MolA)%XYZ(j,k),k=1,3)
        END do
        END do
        WRITE(333,*)'----------------------------------------------------'
!        
		close(333)
        end if		!ThreeD_Change
!
        WRITE(*,*)' Molecular Quantum Similarity Suite'
        write(*,*)' Input and Prepare Data Success!'
        Pause
        Stop
!

	end                !program
!.........................................................................
!.........................................................................
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
	Subroutine ASA_Overlap_Self(Mol_Set,Mol,SSP)
!
!	compute position overlap similarity element with translated-rotated coordinates
!
    use set_molecule
!
    parameter(MOLMAX=60)  			!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 :: SSP
    double precision :: RAB,term,AA,BB,alam,Fecit,SIMOVLP
!
	Type(Molecule)Mol_Set(MOLMAX)
!
!       Position overlap similarity between s-type functions
!
        SSP=zero
        NUMA=Mol_Set(Mol)%N_atoms
        do i=1,NUMA
        nfi=Mol_Set(Mol)%NASAfun(i)        
        do j=1,i
        nfj=Mol_Set(Mol)%NASAfun(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(Mol_Set(Mol)%XYZ(i,kik)-Mol_Set(Mol)%XYZ(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        AA=Mol_Set(Mol)%Ep(i,k)
        do l=1,nfj
        BB=Mol_Set(Mol)%Ep(j,l)
!
        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(Mol)%Op(i,k)*Mol_Set(Mol)%Op(j,l)*SIMOVLP
!
        end do
        end do
        IF(i/=j)term=term+term
        SSP=SSP+term
        end do
        end do
!
		end        ! ASA overlap Selfsimilarity
!------------------------------------------------------------------------------ 
!--------------------------------------------------------------------------
	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 JADIAG(N,NN,O,D,V)
!-----------------------------------------------------------------------
!     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),V(n,n)
	allocatable:: A(:), Z(:)
!
	IF (N.LE.0) STOP ' Bad dimension in JADIAG'
!
! Transfers matrix
!
! preserves matrix to be diagonalized
!
        ALLOCATE (A(NN),Z(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  greater to lesser 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)
!
	if(iter>=90)then
    write(*,*)'***----Jacobi Diagonalization Warning----***'
    write(*,*)'***----Procedure Sweep Number is large---***=',iter
    end if
!
	END subroutine JADIAG
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
