!     Last change:  RCD  11 Aug 2009   11:30 pm
!
	Program ASA Gradient Similarity
	implicit double precision(a-h,o-z),INTEGER*4 (i-n)
	Parameter (PI=3.14159265358979D0)
        PARAMETER (Zero=0.0d0)
        PARAMETER (toau=0.529177249D0)
	parameter (NASAMAX=10);PARAMETER(NATOMX=200)
        Allocatable DPM(:,:)
        DIMENSION SUS(3,3)
	Dimension Nfunm(NATOMX),natomsm(NATOMX)
	dimension Am(NASAMAX,NATOMX),Cm(NASAMAX,NATOMX)
	Dimension Nfunp(NATOMX),natomsp(NATOMX)
	dimension Ap(NASAMAX,NATOMX),Cp(NASAMAX,NATOMX)
        DIMENSION Charp(NATOMX),Charm(NATOMX)
	Dimension xyzp(NATOMX,3),xyzm(NATOMX,3)
        DIMENSION CCHP(3),CCHM(3)
        DIMENSION tran(3),xyzm0(NATOMX,3)
	Character*9 ASAFUN
        CHARACTER*30 moltitp,moltitm
	logical nollegit,fine,LLORO
!
        LLORO=.false.
	write(*,*)' Enter the ASA file: 1= 3-21g OR 2=6-311 OR 3=Huzinaga H-Xe'
	read(*,*)nix
	nollegit=.true.
	do while (nollegit)
        IF(nix.eq.1)then
         ASAFUN='3-21g.txt'
	  nollegit=.false.
	 elseif(nix.eq.2)then
	  ASAFUN='6-311.txt'
	  nollegit=.false.
          ELSEIF(nix.eq.3)then
	  ASAFUN='Huzin.txt'
          nollegit=.false.
	 else
	  write(*,*)' Enter the ASA file: 1= 3-21g OR 2=6-311'
	  read(*,*)nix
	 endif
	enddo
!
!       Enter molecules and attach ASA functions to each atom
!
	Open(1,file=ASAFUN)
	open(2,file='Moleculep.txt')
        READ(2,'(A30)')Moltitp
	READ(2,*)moleculesp,ltoaup
	do ia=1,moleculesp
	 READ(2,*)natomsp(ia),(xyzp(ia,j),j=1,3)
	 nz=natomsp(ia)
!
	 do
	  read(1,*)iz,nf
          IF(nix<=2)then
	  read(1,*)(Cp(i,ia),Ap(i,ia),i=1,nf)
          else
	  read(1,*)(Ap(i,ia),Cp(i,ia),i=1,nf)
          END if
          IF(iz==nz)THEN; Nfunp(ia)=nf;
          do i=1,nf; Ap(i,ia)=Ap(i,ia)
          IF(nix<=2)Cp(i,ia)=Cp(i,ia)/natomsp(ia)
          END do
          EXIT; END if
	enddo
        REWIND(1)
!
        END do !ia
        CLOSE(1)
!
        IF(ltoaup/=0)then
        do ia=1,moleculesp
        do j=1,3
        xyzp(ia,j)=xyzp(ia,j)/toau
        end do
        end do
        END if
!
         READ(2,*)nxar
         IF(nxar/=0)THEN
          do i=1,moleculesp; READ(2,*)Charp(i); END DO
          ELSE
          do i=1,moleculesp; Charp(i)=natomsp(i); END DO
          END if
	close(2)
!
        IF(LLORO)then
        WRITE(*,'(A30)')moltitp
	write(*,*)'Molecule A: ',moleculesp
	do ia=1,moleculesp
         WRITE(*,*)'Atom: ',ia
	 WRITE(*,*)natomsp(ia),(xyzp(ia,j),j=1,3)
	 nz=natomsp(ia); nf=nfunp(ia)
	  Write(*,*)nz,nf
	  write(*,*)(Cp(i,ia),Ap(i,ia),i=1,nf)
	enddo
        endif
!
	Open(1,file=ASAFUN)
	open(2,file='Moleculem.txt')
        READ(2,'(A30)')Moltitm
	READ(2,*)moleculesm,ltoaum
	do ia=1,moleculesm
	 READ(2,*)natomsm(ia),(xyzm(ia,j),j=1,3)
	 nz=natomsm(ia)
!
	 do
	  read(1,*)iz,nf
          IF(nix<=2)then
	  read(1,*)(Cm(i,ia),Am(i,ia),i=1,nf)
          else
	  read(1,*)(Am(i,ia),Cm(i,ia),i=1,nf)
          END if
          IF(iz==nz)THEN; Nfunm(ia)=nf;
          do i=1,nf; Am(i,ia)=Am(i,ia)
          IF(nix<=2)Cm(i,ia)=Cm(i,ia)/natomsm(ia)
          END do
          EXIT; END if
   	  enddo
          REWIND(1)
!
        END do
	close(1)
!
        IF(ltoaum/=0)then
        do ia=1,moleculesm
        do j=1,3
        xyzm(ia,j)=xyzm(ia,j)/toau
        end do
        end do
        END if
!

         READ(2,*)nxar
         IF(nxar/=0)THEN
          do i=1,moleculesm; READ(2,*)Charm(i); END DO
          ELSE
          do i=1,moleculesm; Charm(i)=natomsm(i); END DO
          END if
	close(2)
!
        IF(LLORO)then
        WRITE(*,'(A30)')moltitm
	write(*,*)'Molecule B: ',moleculesm
	do ia=1,moleculesm
        WRITE(*,*)'Atom: ',ia
	 WRITE(*,*)natomsm(ia),(xyzm(ia,j),j=1,3)
	 nz=natomsm(ia);nf= Nfunm(ia)
	  Write(*,*)nz,nf
	  write(*,*)(Cm(i,ia),Am(i,ia),i=1,nf)
	  Nfunp(ia)=nf
	enddo
        END if
!
!       Test the ASA coefficients for each atom in both molecules
!
        ssip=zero;Testp=zero ;szp=zero; ccp=zero
	do j=1,moleculesp
         tistp=zero
         szp=szp+natomsp(j) ;ccp=ccp+Charp(j)
        sip=zero
	 do i=1,Nfunp(j)
         sip=sip+Cp(i,j)
         Cp(i,j)=Cp(i,j)*Charp(j)
          tistp=tistp+Cp(i,j)
	 enddo
         ssip=ssip+sip
         Testp=Testp+tistp
	enddo
        ssim=zero;Testm=zero ;szm=zero; ccm=zero
	do j=1,moleculesm
        tistm=zero
        sim=zero;szm=szm+natomsm(j);ccm=ccm+Charm(j)
	 do i=1,Nfunm(j)
         sim=sim+Cm(i,j)
         Cm(i,j)=Cm(i,j)*Charm(j)
          tistm=tistm+Cm(i,j)
          enddo
          ssim=ssim+sim
          testm=testm+tistm
	enddo
        WRITE(*,*)' Total Molecular number of atomsfrom raw coefficients:  ',ssip,ssim
        WRITE(*,*)' Total Molecular number of electrons from atomic numbers: ',szp,szm
        WRITE(*,*)' Total Molecular number of electrons from chargess: ',ccp,ccm
        WRITE(*,*)' Test of Coefficient coherence: ',Testp,Testm
!
!       compute similarity
!
!       Put both molecules with origin at center of charges
!
        WRITE(*,'(A30)')moltitp
!
        IF(moleculesp>1)then
        Ztotp=zero; do ia=1,moleculesp; Ztotp=Ztotp+Charp(ia);
        do k=1,3; CCHP(k)=CCHP(k)+Charp(ia)*xyzp(ia,k); end do
        END do
        do k=1,3
        CCHP(k)=CCHP(k)/Ztotp
        do ia=1,moleculesp
        xyzp(ia,k)=xyzp(ia,k)-CCHP(k)
        end do
        end do
        WRITE(*,*)' Center of Charges A: '
        WRITE(*,*)(CCHP(i),i=1,3)
!
        END if
!
        WRITE(*,'(A30)')moltitm
!
        IF(moleculesm>1)then
        Ztotm=zero; do ia=1,moleculesm; Ztotm=Ztotm+Charm(ia);
        do k=1,3; CCHM(k)=CCHM(k)+Charm(ia)*xyzm(ia,k); end do
        END do
        do k=1,3
        CCHM(k)=CCHM(k)/Ztotm
        do ia=1,moleculesm
        xyzm(ia,k)=xyzm(ia,k)-CCHM(k)
        end do
        end do
!
        WRITE(*,*)' Center of Charges B: '
        WRITE(*,*)(CCHM(i),i=1,3)
        END if
!
!       save original coordinates of molecule B
!
        do kiki=1,3
        do j=1,moleculesm
        xyzm0(j,kiki)=xyzm(j,kiki)
        end do
        end do
!
!       Compute overlap and gradient selfsimilarity between s-type functions
!
!       compute overlap selfsimilarities
!
        SSOVLP=zero
        do i=1,moleculesp
        nfi=Nfunp(i)
        do j=1,i
        nfj=Nfunp(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(xyzp(i,kik)-xyzp(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        do l=1,nfj
        term=term+Cp(k,i)*Cp(l,j)*SIMOVLP(RAB,Ap(k,i),Ap(l,j))
        end do
        end do
        IF(i/=j)term=term+term
        SSOVLP=SSOVLP+term
        end do
        end do
        WRITE(*,*)' Overlap selfsimilarity Molecule A: ',SSOVLP
        WRITE(*,*)'                   Scaled by Na**2: ',SSOVLP/szp**2

        SSOVLM=zero
        do i=1,moleculesm
        nfi=Nfunm(i)
        do j=1,i
        nfj=Nfunm(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(xyzm(i,kik)-xyzm(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        do l=1,nfj
        term=term+Cm(k,i)*Cm(l,j)*SIMOVLP(RAB,Am(k,i),Am(l,j))
        end do
        end do
        IF(i/=j)term=term+term
        SSOVLM=SSOVLM+term
        end do
        end do
        WRITE(*,*)' Overlap selfsimilarity Molecule B: ',SSOVLM
        WRITE(*,*)'                   Scaled by Nb**2: ',SSOVLM/szm**2
!
!       compute gradient selfsimilarities
!
        SSgradP=zero
        do i=1,moleculesp
        nfi=Nfunp(i)
        do j=1,i
        nfj=Nfunp(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(xyzp(i,kik)-xyzp(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        do l=1,nfj
        term=term+Cp(k,i)*Cp(l,j)*SIMgrad(RAB,Ap(k,i),Ap(l,j))
        end do
        end do
        IF(i/=j)term=term+term
        SSgradP=SSgradP+term
        end do
        end do
        WRITE(*,*)' Gradient selfsimilarity Molecule A: ',SSgradP
        WRITE(*,*)'                    Scaled by Na**2: ',SSgradP/szp**2
        SSgradM=zero
        do i=1,moleculesm
        nfi=Nfunm(i)
        do j=1,i
        nfj=Nfunm(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(xyzm(i,kik)-xyzm(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        do l=1,nfj
        term=term+Cm(k,i)*Cm(l,j)*SIMgrad(RAB,Am(k,i),Am(l,j))
        end do
        end do
        IF(i/=j)term=term+term
        SSgradM=SSgradM+term
        end do
        end do
        WRITE(*,*)' Gradient selfsimilarity Molecule B: ',SSgradM
        WRITE(*,*)'                    Scaled by Nb**2: ',SSgradM/szm**2
!
!*****************************************
!       Similarities
!*****************************************
!
!       compute distances between atoms of both molecules A-B
!
!       try to optimize moving the center of charge of Molecule B
!       over the Z axis only
!
!
        WRITE(*,*)'Enter : Tmin, Tmax, DeltaT: '
        READ(*,*)Tmin,Tmax,DeltaT
        akuku=Tmin
        sakamax=zero; Sobmax=zero
        gakamax=zero; Gobmax=zero
        do
!
!       move molecule B
!
        tran=0.0d0;
        WRITE(*,*)' translation: ',akuku
        tran(3)=tran(3)+akuku
!
        do j=1,moleculesm
        do itran=1,3 ; xyzm(j,itran)=xyzm0(j,itran)+tran(itran); end do
        end do
!
        ALLOCATE(DPM(moleculesp,moleculesm))
        do ia=1,moleculesp
        do ib=1,moleculesm
        DPM(ia,ib)=zero
        do k=1,3;DPM(ia,ib)=DPM(ia,ib)+(xyzp(ia,k)-xyzm(ib,k))**2;end do
        end do
        end do
        distAB=SUM(DPM)
        WRITE(*,*)'SuM of A-B Distances: ',DistAB
        DEALLOCATE(DPM)
!
!        Compute matrix of condensed coordinates
!
        Tootsy=zero
	do k=1,3
        do l=1,3
        SUS(k,l)=zero
        do ia=1,moleculesp
        do ib=1,moleculesm
        SUS(k,l)=SUS(k,l)+xyzp(ia,k)*xyzm(ib,l)
        end do
        end do
        Tootsy=Tootsy+SUS(k,l)
        END do
        END do
!
        IF(lloro)then
        WRITE(*,*)'Matrix of A-B Scalar products'
        do k=1,3
        WRITE(*,*)(SUS(l,k),l=1,3)
        end do
        END if
        WRITE(*,*)' Sum of A-B scalar products: ',Tootsy
!
!       Overlap Similarity
!
	SOVLPM=zero
        do i=1,moleculesp
        nfi=Nfunp(i)
        do j=1,moleculesm
        nfj=Nfunm(j)
        RAB=zero
        do kik=1,3;RAB=RAB+(xyzp(i,kik)-xyzm(j,kik))**2;end do
        term=zero
        do k=1,nfi
        do l=1,nfj
        term=term+Cp(k,i)*Cm(l,j)*SIMOVLP(RAB,Ap(k,i),Am(l,j))
        end do
        end do
        SOVLPM=SOVLPM+term
        end do
        end do
        WRITE(*,*)' Overlap similarity Molecules A-B: ',SOVLPM
        WRITE(*,*)'                  Scaled by NA*NB: ',SOVLPM/(szp*szm)
        IF(SOVLPM>Sobmax)THEN;Sobmax=SOVLPM;sakamax=akuku;END IF
!
        rpm=SOVLPM/SQRT(SSOVLP*SSOVLM)
        WRITE(*,*)'Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(SSOVLP+SSOVLM-2*SOVLPM)
        WRITE(*,*)'Overlap Distance Index A-B: ',distpm
!
!       Gradient Similarity
!
	SgradPM=zero
        do i=1,moleculesp
        nfi=Nfunp(i)
        do j=1,moleculesm
        nfj=Nfunm(j)
        RAB=zero
        do kik=1,3;RAB=RAB+(xyzp(i,kik)-xyzm(j,kik))**2;end do
        term=zero
        do k=1,nfi
        do l=1,nfj
        term=term+Cp(k,i)*Cm(l,j)*SIMgrad(RAB,Ap(k,i),Am(l,j))
        end do
        end do
        SgradPM=SgradPM+term
        end do
        end do
        WRITE(*,*)' Gradient similarity Molecules A-B: ',SgradPM
        WRITE(*,*)'                   Scaled by NA*NB: ',SgradPM/(szp*szm)
        IF(Gobmax<SgradPM)THEN;Gobmax=SgradPM;gakamax=akuku;END IF
!
        rgpm=SgradPM/SQRT(SSgradP*SSgradM)
        WRITE(*,*)'Gradient Similarity Index A-B: ',rgpm
!
        distgpm=SQRT(SSgradP+SSgradM-2*SgradPM)
        WRITE(*,*)'Gradient Distance Index A-B: ',distgpm
!
        akuku=akuku+DeltaT
        IF(akuku>Tmax)EXIT
        END do !akuku translation
!
        WRITE(*,*)' '
        WRITE(*,*)'****************************'
        WRITE(*,*)'MAXIMA'
        WRITE(*,*)'****************************'
        WRITE(*,*)' Displacement : ', Sakamax
        WRITE(*,*)' Overlap and overlap scaled : ',Sobmax,Sobmax/(szp*szm)
        rpm=Sobmax/SQRT(SSOVLP*SSOVLM)
        WRITE(*,*)'Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(SSOVLP+SSOVLM-2*Sobmax)
        WRITE(*,*)'Overlap Distance Index A-B: ',distpm
!
        WRITE(*,*)' Displacement : ', Gakamax
        WRITE(*,*)' Gradient and gradient scaled: ', Gobmax,Gobmax/(szp*szm)
	rgpm=Gobmax/SQRT(SSgradP*SSgradM)
        WRITE(*,*)'Gradient Similarity Index A-B: ',rgpm
        distgpm=SQRT(SSgradP+SSgradM-2*Gobmax)
        WRITE(*,*)'Gradient Distance Index A-B: ',distgpm
!
!
	End !Program Gradient Similarity
!
	DOUBLE PRECISION function SIMOVLP(RAB2,AA,BB)
        implicit double PRECISION (a-h,o-z)
        PARAMETER (toler=1.0D-7)
	Parameter (PI=3.14159265358979D0)
!
!       overlap similarity between s-type functions
!
        alam=2*AA*BB/(AA+BB); Fecit=SQRT(alam/PI)**3
        IF(ABS(RAB2)<=toler)then
        SIMOVLP=Fecit
        else
        SIMOVLP=Fecit*EXP(-alam*RAB2)
        END if
        end !SIMOVLP
!
!
	DOUBLE PRECISION function SIMgrad(RAB2,AA,BB)
        implicit double PRECISION (a-h,o-z)
        PARAMETER (toler=1.0D-7, threehalf=1.5d0)
	Parameter (PI=3.14159265358979D0,PI3=PI**3)
!
!       gradient similarity between s-type functions
!
        IF(ABS(RAB2)<=toler)then
        alam=(AA*BB/(AA+BB))**5; alf=SQRT(2*alam/pi3)
        SIMgrad=24*alf
        else
        alam=2*AA*BB/(AA+BB); alf=alam**5; bel=SQRT(alf/PI3)
        alamillo=alam*RAB2
        SIMgrad=4*bel*(threehalf-alamillo)*EXP(-alamillo)
        END if
        end !SIMgrad
!----------------------------------------------------------------------
!-------------------------------------------------------------------------------
      subroutine Jadia (n,nn,a,d,v, testis)
!
!     sYMMETRIC MATRIX DIAGONALIZATION BY JACOBI'S METHOD
!
      implicit double precision (A-H, O-Z)

      parameter( ZERO=0.0D0, TENM9=1.0D-9, ONE=1.0D0, HALF=0.5D0, POTWO=2.0D0)

      Dimension a(nn), d(n), v(n,n), z(n)

!       n              matrix dimension
!       nn             (n*n+n)/2
!       a(nn)          original matrix
!
!       v(n,n)         eigenvector matrix
!       d(n)           eigenvalue matrix
!
!       z(n)           auxiliary matrix
!
!       testis         controller of eigenvalue ordering
!			='S' : ascending
!			='L' : descending

      character*1 testis
      integer     p, q
      logical     dold, equil, equis

      if (n.le.0) return

      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.gt.1) then
        i=0
        sm=ONE
69	if ((sm.gt.TENM9) .and. (i.lt.100)) then
	  i=i+1

!	  Convergence criterion
          ipp=1
          sm=ZERO
	  do p=2, n
            ipq=ipp
	    do q=1, p-1
              ipq=ipq+1
	      sm=sm+dabs (a(ipq))
	    end do
	    ipp=ipp+p
	  end do
	  if (i.lt.4) then
            tresh=POTWO*sm/(n*n)
          else
            tresh=ZERO
	  endif

!	  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*dabs (a(ipq))
              if (i.gt.4 .and. (dabs (d(p))+g).eq.dabs (d(p))&
               .and. (dabs (d(q))+g).eq.dabs (d(q))) then
                a(ipq)=ZERO
              else
                if (dabs (a(ipq)).gt.tresh) then
                  h=d(q)-d(p)
                  if ((dabs (h)+g).le.dabs (h)) then

!		   Generate the shift for Jacobi rotation
                    t=a(ipq)/h
                  else
                    theta=HALF*h/a(ipq)
                    t=ONE/(dabs (theta)+dsqrt (ONE+theta*theta))
                    if (theta.lt.ZERO) t=-t
                  endif
                  c=ONE/dsqrt (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.ne.p) then
                      if (j.lt.p) then
                        jp=ipp+j
                      else
                        jp=jj+p
                      endif
                      if (j.ne.q) then
                        if (j.lt.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
	  go to 69
        end if

!	Order eigenvalues and eigenvectors
        equil=testis.eq.'L'
        equis=testis.eq.'S'
	do p=1, n-1
	  do q=p+1, n
            if (equil) dold=d(q).gt.d(p)
            if (equis) dold=d(q).lt.d(p)
            if (dold) 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
      endif

      end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      subroutine Q$set(MXAT,CA,C,COT)
!
!     Computes quadrupole moment and transforms coordinates
!     to principal axis components
!

      implicit double precision (A-H, O-Z)
      character*1 croc

      parameter (NEUDIM=3, CROC='S')
      parameter (NENE=(NEUDIM*NEUDIM+NEUDIM)/2)
      parameter (ZERO=0.0D0)
!
      Dimension ca(MXAT)
      dimension c(MXAT,NEUDIM)
      dimension quad(NENE), d(NEUDIM), z(neudim),v(NEUDIM,NEUDIM)
!
!	Quadrupole moment from charges
        kl=0
	do k=1, NEUDIM
	  do l=1,k
	  kl=kl+1
	  quad(kl)=ZERO
	    do i=1, MXAT
 	      quad(kl)=quad(kl)+c(i,k)*c(i,l)*ca(i)
            end do
          end do
        end do

!	Quadrupole diagonalization
	call Jadia (NEUDIM,NENE,quad,d,v,CROC)
!
!	Transformation of coordinates to
!       quadrupole principal axis
	do i=1, MXAT
	  do j=1, NEUDIM
            d(j)=c(i,j)
          end do

          do k=1, NEUDIM
	    z(k)=ZERO
	    do l=1, NEUDIM
	      z(k)=z(k)+v(l,k)*d(l)
            end do
          end do

	  do j=1, NEUDIM
            c(i,j)=z(j)
          end do
        END do
      end
!-------------------------------------------------------------------------------

