!     Last change:  RCD  18 Nov 2009    6:43 pm
!
	Program ASA Constrained Overlap triple density Similarity for atoms
	implicit double precision(a-h,o-z),INTEGER*4 (i-n)
	Parameter (PI=3.14159265358979D0)
        PARAMETER (Zero=0.0d0)
        PARAMETER(Half=0.5d0)
        PARAMETER(Toler=1.0d-7)
	parameter (NASAMAX=10)
	dimension Ap(NASAMAX),Cp(NASAMAX)
	Character*9 ASAFUN
        CHARACTER*16 Kali_Durga
	logical nollegit,LLORO,Komput
        LLORO=.true.
        Komput=.true.
!
        lumen=0
        do
!
	write(*,*)' Enter the ASA file: 1= 3-21g OR 2=6-311'
	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.
	 else
	  write(*,*)' Enter the ASA file: 1= 3-21g OR 2=6-311'
	  read(*,*)nix
	 endif
	enddo
!
!       attach ASA functions to each atom
!
        WRITE(*,*)' Read atom nuclear charge: '
	READ(*,*)nzp
!
	Open(1,file=ASAFUN)
	 do
	  read(1,*)iz,nf
	  read(1,*)(Cp(i),Ap(i),i=1,nf)
          IF(iz==nzp)THEN; Nfunp=nf; EXIT; END if
	enddo
        REWIND(1)
        CLOSE(1)
!
!       write information on ASA functions
!
        IF(LLORO)then
	  Write(*,*)nzp,Nfunp
          check=zero
          do i=1,Nfunp
          check=check+Cp(i)
	  write(*,*)i,Cp(i),Ap(i)
          end do
          WRITE(*,*)' Check coefficients: ',check
          WRITE(*,*)'-'
        endif
!
        Lumen=Lumen+1; IF(Lumen>25)exit
        call Philename(lumen,nzp,ASAFUN,Kali_Durga)
        write(*,*)' File name:'
        WRITE(*,'(1x,a15)')Kali_Durga
        OPEN(UNIT=7,File=Kali_Durga)
        WRITE(7,'(1x,a15)')Kali_Durga
!
!       Compute overlap triple density selfsimilarity between s-type functions
!
        SSOVLP=ZERO
        do k=1,nfunp
        AA=Ap(k)!/2
        CPK=Cp(k)
        do l=1,nfunp
        BB=Ap(l)!/2
        CPL=Cp(l)
        tirm=CPK*CPL
        do m=1,nfunp
        CC=Ap(m)
        term=tirm*Cp(m);
        alam=4*AA*BB*CC/(AA+BB+CC); Fecit=(SQRT(alam)/PI)**3
        termp=term*Fecit
        SSOVLP=SSOVLP+termp
        end do
        end do
        END do
!
        WRITE(*,*)' Overlap selfsimilarity : ',SSOVLP
        WRITE(*,*)'           in Z**3 units: ',SSOVLP/NZP**3
        WRITE(*,*)'*'
        WRITE(7,*)' Overlap selfsimilarity : ',SSOVLP
        WRITE(7,*)'           in Z**3 units: ',SSOVLP/NZP**3
        WRITE(7,*)'*'

!
!       constrained triple selfsimilarity for atoms
!
        IF(komput)then
!
        do KR=1,8
        R=KR*half
        WRITE(*,*)'******'
        WRITE(*,*)' Constraining Radius : ',R
!
        WRITE(7,*)'******'
        WRITE(7,*)' Constraining Radius : ',R
!
        SSKOVLP=ZERO
        do k=1,nfunp
        AA=Ap(k)!/2
        CPK=Cp(k)
        do l=1,nfunp
        BB=Ap(l)!/2
        tirm=CPK*Cp(l)
        do m=1,nfunp
        CC=Ap(m)
        term=tirm*Cp(m)
        sab=AA+BB+CC; RAB=R*SQRT(2*sab)
        alam=4*AA*BB*CC/SAB; Fecit=(SQRT(alam)/PI)**3
        termp=term*Fecit*erf(RAB)**3
        SSKOVLP=SSKOVLP+termp
        end do
        end do
        END do
!
        DSK=SSOVLP-SSKOVLP
        WRITE(*,*)' Constrained Overlap selfsimilarity  : ',SSKOVLP 
        WRITE(*,*)'           in Z**3 units: ',SSKOVLP/NZP**3
        WRITE(*,*)' Difference Unconstrained-Constrained: ',DSK
!
        WRITE(*,*)'*'
!
        WRITE(7,*)' Constrained Overlap selfsimilarity : ',SSKOVLP
        WRITE(7,*)'           in Z**3 units: ',SSKOVLP/NZP**3
        WRITE(7,*)' Difference Unconstrained-Constrained: ',DSK
!
        WRITE(7,*)'*' ; IF(DSK<toler)exit
!
        END do
        END if !Komput
        CLOSE(7)
!
        END do
!
	End !Program Overlap and Gradient Similarity for atoms
!
        subroutine Philename(Kip,Nz,Philis,Black_goddess)
!
!       Constructs a .txt file name (Black_goddess) according to
!	atom Z (Nz), run number<= (Kip) and ASA basis set Id (Philis)
!
	CHARACTER*2 Korda(25)
        DATA Korda/'1_','2_','3_','4_','5_','6_','7_','8_','9_','10',&
                   '11','12','13','14','15','16','17','18','19','20',&
                   '21','22','23','24','25'/
        CHARACTER*2 Period(54)
        data Period/'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'/
        CHARACTER*16 Black_goddess
        CHARACTER*9 Philis
!
        Black_Goddess=Period(Nz)//Korda(Kip)//'_'//Philis
!
        end subroutine

