!     Last change:  RCD  13 Nov 2009   10:42 pm
!
	Program ASA Bounded Molecular Similarity
	implicit double precision(a-h,o-z),INTEGER*4 (i-n)
	Parameter (PI=3.14159265358979D0)
        PARAMETER (Zero=0.0d0, half=0.5d0, one=1.0d0, two=2.0d0)
        parameter (Toau=0.529177249D0)
	parameter (NASAMAX=10);PARAMETER(NATOMX=200, NATNAT=(NATOMX*NATOMX+NATOMX)/2)
	Dimension Nfunp(NATOMX),natomsp(NATOMX)
	dimension Ap(NASAMAX,NATOMX),Cp(NASAMAX,NATOMX)
        DIMENSION Cp0(NASAMAX,NATOMX)
        DIMENSION Charp(NATOMX)
	Dimension xyzp(NATOMX,3)
        DIMENSION ZS(NATOMX,NATOMX)
        DIMENSION PS(NATOMX)
        DIMENSION QS(NATOMX)
        DIMENSION ZSZS(NATNAT),D(NATOMX),V(NATOMX,NATOMX),QQ(NATOMX)
	Character*9 ASAFUN
        CHARACTER*20 typet
        CHARACTER*30 moltitp
        INTEGER Parrot
	logical LLORO
!
        LLORO=.false.
!
!       Enter molecule and attach ASA functions to each atom
!       the file is organized as follows:
!       1) Title (30 characteres)
!       2) nix: (ASA file) 1 (3-21G)   2 (6-311G)    3 (Huzinaga) and Parrot: /=0 verbose output  =0 normal output
!       3) Moleculesp (Number of atoms) and ltoau (transform to au (1) or not (0))
!	4) For each atom: Atomic Number, Coordinates X,Y,Z
!	5) Nxar: /=0 reads atom populations; =0 uses atomic number
!	6) when Nxar/=0: atom populations in a column (same order than coordinates) otherwise nothing.
!
	open(2,file='Molecule.txt')
        READ(2,'(A30)')Moltitp
	read(2,*)nix,Parrot
        IF(Parrot/=0)LLORO=.true.
!
	 if(nix.eq.1)then
	  ASAFUN='3-21g.txt'
	 elseif(nix.eq.2)then
	  ASAFUN='6-311.txt'
	 elseif(nix.eq.3)then
	  ASAFUN='Huzin.txt'
	 else
          WRITE(*,*)' Error in ASA fila specification!'
	  write(*,*)' Please, enter the ASA file: 1= 3-21g OR 2=6-311 OR 3=Huzinaga H-Xe'
          Stop
	 endif
	Open(1,file=ASAFUN)

!
	READ(2,*)moleculesp, ltoau
	do ia=1,moleculesp
	READ(2,*)natomsp(ia),(xyzp(ia,j),j=1,3)
	nz=natomsp(ia)
!
        IF(ltoau/=0)then
        do i=1,Moleculesp
        do j=1,3
        xyzp(i,j)=xyzp(i,j)/toau
        end do
        end do
        END if
!
          do  ! fetch ASA coefficients and exponents for each atom Z(=nz)
!
	  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)
!
        WRITE(*,'(A30)')moltitp
         WRITE(*,'(a40)')' ASA basis set: '//ASAFUN
         READ(2,'(I1,a20)')nxar,typet
         IF(nxar/=0)THEN
         WRITE(*,'(a40)')' Populations type: '//typet
          do i=1,moleculesp; READ(2,*)Charp(i); END DO
          ELSE
         WRITE(*,*)' Populations type: Atomic Number'
          do i=1,moleculesp; Charp(i)=natomsp(i); END DO
          END if
!
	close(2)
!
        IF(LLORO)then
        WRITE(*,'(A30)')moltitp
	write(*,*)'Molecule : ',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
!
     	close(1)
!
!       Test the ASA coefficients for each atom in molecule
!
        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
!
        WRITE(*,*)' Total Molecular number of atomsfrom raw coefficients:  ',ssip
        WRITE(*,*)' Total Molecular number of electrons from atomic numbers: ',szp
        WRITE(*,*)' Total Molecular number of electrons from chargess: ',ccp
        WRITE(*,*)' Test of Coefficient coherence: ',Testp
!
!       compute constrained
!       similarity matrix for the associated molecule
!       overlap selfsimilarity between s-type functions
!
        call SELFOVERLAP(Moleculesp,Nfunp,xyzp,Ap,Cp,ZS,SSOVLP,unck)
!
        WRITE(*,*)' Unconstrained Calculations'
        WRITE(*,*)' Overlap selfsimilarity:    ',SSOVLP
        WRITE(*,*)'            Scaled by N**2: ',SSOVLP/szp**2
        TATA=SQRT(SSOVLP)
        WRITE(*,*)' SQRT Overlap selfsimilarity: ',TATA
        TITUS=SQRT(TATA)
        WRITE(*,*)' SQRT(SQRT) Overlap selfsimilarity: ',TITUS
!
!       Constrained calculations
!
        WRITE(*,*)' Constrained Calculations'
!
        RORO=half
        do kap=1,8
        call CSELFOVERLAP(Moleculesp,Nfunp,xyzp,Ap,Cp,ZS,SSOVLP,RORO)
!
        WRITE(*,*)' Cicumscribed sphere radius= ',RORO
        WRITE(*,*)'Constrained Overlap selfsimilarity:    ',SSOVLP
        WRITE(*,*)'            Scaled by N**2: ',SSOVLP/szp**2
!
        RORO=RORO+half
        END do !kap
!
      	End ! Program Overlap Similarity Populations
! --------------------------------------------------------------------------------------------------------
	DOUBLE PRECISION function SIMOVLP(RAB2,AA,BB)
        implicit double PRECISION (a-h,o-z)
        PARAMETER (toler=1.0D-7)
	Parameter (PI=3.14159265358979D0)
!
!       overlap unconstrained similarity integrals 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 CSIMOVLP(RAB2,AA,BB,ROCK)
        implicit double PRECISION (a-h,o-z)
        PARAMETER (toler=1.0D-7)
	Parameter (PI=3.14159265358979D0)
!
!       overlap constrained similarity integrals
!	between s-type functions
!       in fact it is the unconstained overlap scaled by
!	a cube of the error function
!
        summum=AA+BB
        alam=2*AA*BB/summum
        Fecit=SQRT(alam/PI)**3
        IF(ABS(RAB2)<=toler)then
        SIMOVLP=Fecit
        else
        SIMOVLP=Fecit*EXP(-alam*RAB2)
        END if
        ROCKY=ROCK*SQRT(summum+summum)
        BOXY=erf(ROCKY)
        CSIMOVLP=SIMOVLP*BOXY**3
        end !CSIMOVLP
!------------------------------------------------------------------------------------------
        SUBROUTINE CSELFOVERLAP(Molecules,Nfun,xyz,A,C,ZS,SSOVL,ROCK)
!
!       computes unconstrained/constrained overlap selfsimilarity of a given molecule
!
        implicit double PRECISION (a-h,o-z)
        PARAMETER (zero=0.0d0)
       	parameter (NASAMAX=10);PARAMETER(NATOMX=200)
	Dimension Nfun(NATOMX)
	dimension A(NASAMAX,NATOMX),C(NASAMAX,NATOMX)
	Dimension xyz(NATOMX,3)
        DIMENSION ZS(NATOMX,NATOMX)
!
        SSOVL=zero
        do i=1,molecules
        nfi=Nfun(i)
        do j=1,i
        nfj=Nfun(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(xyz(i,kik)-xyz(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        do l=1,nfj
!
        	QMSSM=CSIMOVLP(RAB,A(k,i),A(l,j),ROCK)
!
        term=term+C(k,i)*C(l,j)*QMSSM
        end do
        end do
!
        ZS(i,j)=SQRT(term);ZS(j,i)=ZS(i,j)
        IF(i/=j)term=term+term
        SSOVL=SSOVL+term
        end do
        end do
!
	end   !Cselfoverlap
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
        SUBROUTINE SELFOVERLAP(Molecules,Nfun,xyz,A,C,ZS,SSOVL)
!
!       computes unconstrained/constrained overlap selfsimilarity of a given molecule
!
        implicit double PRECISION (a-h,o-z)
        PARAMETER (zero=0.0d0)
       	parameter (NASAMAX=10);PARAMETER(NATOMX=200)
	Dimension Nfun(NATOMX)
	dimension A(NASAMAX,NATOMX),C(NASAMAX,NATOMX)
	Dimension xyz(NATOMX,3)
        DIMENSION ZS(NATOMX,NATOMX)
!
        SSOVL=zero
        do i=1,molecules
        nfi=Nfun(i)
        do j=1,i
        nfj=Nfun(j)
        RAB=zero
        IF(i/=j)THEN; do kik=1,3;RAB=RAB+(xyz(i,kik)-xyz(j,kik))**2;end DO; END if
        term=zero
        do k=1,nfi
        do l=1,nfj
!
        	QMSSM=SIMOVLP(RAB,A(k,i),A(l,j))
!
        term=term+C(k,i)*C(l,j)*QMSSM
        end do
        end do
!
        ZS(i,j)=SQRT(term);ZS(j,i)=ZS(i,j)
        IF(i/=j)term=term+term
        SSOVL=SSOVL+term
        end do
        end do
!
	end   !selfoverlap

!-----------------------------------------------------------------------------------------------------------
