!     Last change:  RCD   5 May 2010    2:09 pm
!     
!-----------------------------------------------------------------------------------------------------------
	module Set_Molecule
!
    parameter(NatMX=151)    !Maximum atoms per Molecule
    parameter(maxfn=7)	    !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_similarity_Matrix
	use set_molecule
!
    Integer*4 MolMolMAX
	parameter(MOLMAX=70)  			!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=80):: ASAFUN		!Name of the file for ASA function type
!-------------------------------------------------------------------------------------------------
        Logical:: LLORO						! .true. write more info  
        logical:: Zenith					! .true. constructs file with optimal pair coordinates
        LOGICAL:: azimut                    ! .true. Makes bold optimization
        Logical:: Found_Optimal
!-------------------------------------------------------------------------------------------------
	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 :: SOVLPM,Scaled,rpm,distpm
!
	integer*4 :: MoloM, natmax, nbtmax, Molchk
    integer*4 :: MM, iii, MolMom , K  
	double precision :: ZZ, CRB
    allocatable :: ZZ(:),CRB(:)
    integer IXYZB(3,3)
    Character*70 :: ID_Label_Pre,ID_Label_Suf
    Character*70 :: OUTPUT_FIL,PREPARED_FIL,QSM_OUT,CINDX_OUT,CRBINDX_OUT
    Character*70 :: MOLECOUT,MOLECPRE,QSMatrix,OPT_XYZ,OPT_XYZ_OUT,AUXILY
!--------------------------------------------------------------------------------------------------	
	Type(Molecule)Mol_Set(MOLMAX)
!--------------------------------------------------------------------------------------------------
!
!       read flags in Flagstaff_2.txt
!
		LLORO=.false.
    	Zenith=.false.
    	Azimut=.false.
        OPEN(333,FILE='Flagstaff_2.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.
        CLOSE(333)
!
!		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_2.txt')
!
		Read(304,*)ID_Label_Pre
        Read(304,*)ID_Label_Suf
        Read(304,*)OUTPUT_FIL
        Read(304,*)PREPARED_FIL
        Read(304,*)QSMatrix
        Read(304,*)CINDX_OUT
        READ(304,*)OPT_XYZ
!
        close(304)
!
!	Reads a batch of ASA defined molecules
!	as given by Molecule-Type Input-Output
!
!	Reads molecular data in file MOLECPRE
!
		MOLECPRE=TRIM(ID_LABEL_Pre)//TRIM(PREPARED_FIL)//TRIM(ID_Label_Suf)//'.txt'
		open(232,File=MOLECPRE)
!
                    READ(232,'(a80)')General_Title 
                    Read(232,'(a80)')ASAFUN 
                    READ(232,*)Kind_ASA
!                    
                    READ(232,*)MoloM
!                   
					MM=0					
					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
                    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						!Position Overlap Quantum Self-Similarity
!					
            		END DO
           			close(232)
!------------------------------------------------------------------------------------------------------------------        
!	data retrieved succesfully and in memory now
!-------------------------------------------------------------------------------------------------------------------
!	Starting Sinilarity Matrix Computation
!-------------------------------------------------------------------------------------------------------------------
!	Overlap Matrix in this case
!-------------------------------------------------------------------------------------------------------------------
!    
    	MOLECOUT=TRIM(ID_Label_Pre)//TRIM(OUTPUT_FIL)//TRIM(ID_Label_Suf)//'.txt'
        OPT_XYZ_OUT=TRIM(ID_Label_Pre)//TRIM(OPT_XYZ)//TRIM(ID_Label_Suf)//'.txt'
!        
		Open(777,File=MOLECOUT)
        open(313,File=OPT_XYZ_OUT)
!
!       Start bimolecular quantum similarity measures first translation only
!
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'----------Bimolecular QSM--------------'
        WRITE(777,*)'---------------------------------------'
        Write(777,'(a80)')General_Title 
        Write(777,'(a80)')ASAFUN 
        Write(777,*)' ASA Kind:',Kind_ASA
        Write(777,*)' Number of Molecules : ',MoloM
!
		if(Zenith)then
        WRITE(313,*)'-------------------------------------------------------'
        WRITE(313,*)'----------Bimolecular Optimal Coordinates--------------'
        WRITE(313,*)'-------------------------------------------------------'
        Write(313,'(a80)')General_Title 
        Write(313,'(a80)')ASAFUN 
        Write(313,*)' ASA Kind:',Kind_ASA
        Write(313,*)' Number of Molecules : ',MoloM
        end if   ! Zenith
!
!		construct simple XYZ rotation information
!
        IF(Azimut)then
          IXYZB(:,:)=1

            IXYZB(1,1)=0
            IXYZB(2,2)=0
            IXYZB(3,3)=0
!
         END if  !Azimut
!
!		create auxiliary and input-output filenemes
! 
        QSM_OUT=TRIM(ID_Label_Pre)//TRIM(QSMatrix)//TRIM(ID_Label_Suf)//'.txt'	
        CRBINDX_OUT=TRIM(ID_Label_Pre)//TRIM(CINDX_OUT)//TRIM(ID_Label_Suf)//'.txt'	
        AUXILY=TRIM(ID_Label_Pre)//'AUX_QSM_CRBI'//TRIM(ID_Label_Suf)//'.txt'
!
		open(979,file=AUXILY)
!
		Molchk=0
        do Mol=1,Molom
        do Mom=1,Mol
!
        molmom=Mom+(Mol*Mol-Mol)/2
        if(Mol==Mom)then
          write(979,*)molmom,mol,mom,Mol_Set(Mol)%overlap_QSS,one
          Molchk=Molchk+1
          cycle
          end if	!Mol=Mom
!*
        Scaling=Mol_Set(Mol)%N_electrons*Mol_Set(Mom)%N_electrons	
!       
		WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------Molecules: '
        write(777,'(1x,"-->",i4,"-*-",3x,a80)')Mol,Mol_Set(Mol)%Identity
        write(777,'(1x,"-->",i4,"-*-",3x,a80)')Mom,Mol_Set(Mom)%Identity
        WRITE(777,*)'---------------------------------------'
!
!		Before anything try to find out simple optimal positioning
!
		Sobmax=zero
!
!		initialize spare coordinates
!
        do j=1,Mol_Set(Mol)%N_atoms
        do itran=1,3 ; xyzta(j,itran)=Mol_Set(Mol)%XYZ(j,itran);end do
        do itrana=1,3 ; xyztaa(j,itrana)=Mol_Set(Mol)%XYZ(j,itrana); end do
        end do
!
        do j=1,Mol_Set(Mom)%N_atoms
        do itran=1,3 ; xyztbb(j,itran)=Mol_Set(Mom)%XYZ(j,itran)
        end do
        end do
!
		Call ASA_Overlap(Mol_Set,Mol,Mom,xyzta,xyztbb,SOVLPM)
		sobmax=SOVLPM
!
        IF(Azimut)then
!
        Do JNX=1,3
!
        do j=1,Mol_Set(Mom)%N_atoms
        do itran=1,3 
           If (IXYZB(JNX,ITRAN)==0)then
             xyztb(j,itran)=Mol_Set(Mom)%XYZ(j,itran)
             else
              xyztb(j,itran)=-Mol_Set(Mom)%XYZ(j,itran) 
              end if
              end do
        end do
!
!	compute position overlap similarity element with simply rotated coordinates
!
	Call ASA_Overlap(Mol_Set,Mol,Mom,xyzta,xyztb,SOVLPM)
!    
        IF(SOVLPM>Sobmax)THEN
       Sobmax=SOVLPM
!
!		Put spare coordinates into maximal positioning
!
        do j=1,Mol_Set(Mom)%N_atoms
        do itran=1,3 ; xyztbb(j,itran)=xyztb(j,itran)
       end do
        end do
!        
       END IF  !Max sov
!
		end do  !JNX
!
        END if  !Azimut
!
!       Move heavy atoms of Molecule A over the heavy atoms of molecule B
!
		natmax=0; nbtmax=0
!
		if(LLORO)then
        WRITE(777,*)'---------------------------------------------------'
        WRITE(777,*)'--------------Atomic Translations------------------'
        WRITE(777,*)'---------------------------------------------------'
		end if	!LLORO
!
        do iata=1,Mol_Set(Mol)%N_atoms; naz=Mol_Set(Mol)%Z(iata)
        do ibtb=1,Mol_Set(Mom)%N_atoms; nbz=Mol_Set(Mom)%Z(ibtb)
        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
!
        do itran=1,3 ; tran(itran)=xyztaa(iata,itran);end do   
        do j=1,Mol_Set(Mol)%N_atoms
        do itran=1,3 ; xyzta(j,itran)=xyztaa(j,itran)-tran(itran); end do
        end do
!
        do itran=1,3 ; tran(itran)=xyztbb(ibtb,itran);end do
        do j=1,Mol_Set(Mom)%N_atoms
        do itran=1,3 ; xyztb(j,itran)=xyztbb(j,itran)-tran(itran); end do
        end do
!
!		START ROTATING MOLECULAR FRAMES For Overlap Optimal Search
!
!	    form the global orientation vectors on A and B
!		          use non-hydrogen atoms only
!
        xa=zero
        do ia=1,Mol_Set(Mol)%N_atoms; nza=Mol_Set(Mol)%Z(ia)
        IF(nza>1)then
        do i=1,3; xa(i)=xa(i)+xyzta(ia,i); end do
        END if
        END do
!
        xb=zero
        do ib=1,Mol_Set(Mom)%N_atoms; nzb=Mol_Set(Mom)%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(xa,ra,caa,saa,cza,sza)
		call zalign(xb,rb,cab,sab,czb,szb)
!
!	rotate all coordinates of A and B
!
        do ia=1,Mol_Set(Mol)%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 ib=1,Mol_Set(Mom)%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
!
!	compute position overlap similarity element with translated-rotated coordinates
!
		Call ASA_Overlap(Mol_Set,Mol,Mom,xyzta,xyztb,SOVLPM)
!    
        IF(SOVLPM>Sobmax)THEN
        Sobmax=SOVLPM; natmax=iata; nbtmax=ibtb
          END IF
!
        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
!        
        WRITE(777,*)'---------------------------------------'
        WRITE(777,*)'-----------*OPTIMAL RESULTS*-----------'
        WRITE(777,*)'---------------------------------------'
!
		if(LLORO)then 
        WRITE(777,*)'-------for Molecules: '
        write(777,'(1x,"-->",i4,"-*-",a80)')Mol,Mol_Set(Mol)%Identity
        write(777,'(1x,"-->",i4,"-*-",a80)')Mom,Mol_Set(Mom)%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
        SSmx=Sobmax/Scaling
        WRITE(777,*)' Overlap  scaled             : ',SSmx
        rpm=Sobmax/SQRT(Mol_Set(Mol)%overlap_QSS*Mol_Set(Mom)%overlap_QSS)
        WRITE(777,*)' Overlap Similarity Index A-B: ',rpm
        distpm=SQRT(Mol_Set(Mol)%overlap_QSS+Mol_Set(Mom)%overlap_QSS-2*Sobmax)
        WRITE(777,*)' Overlap Distance Index A-B  : ',distpm
        WRITE(777,*)'---------------------------------------'
!
		Found_Optimal=.True.
!        
        else
        WRITE(777,*)'++++++++++> NO OPTIMAL VALUE FOUND!'
!
        Found_Optimal=.False.
!
            end if
!
!	Form the similarity matrices and preserve optimal TR coordinates 
!						if Zenith .true.
!					for every molecular pair
		if(Zenith)then
        WRITE(313,*)'------------------------------------------------'
        if(Found_Optimal)then
		WRITE(313,*)'-- Optimal Molecular Coordinates Between: '
        else
        WRITE(313,*)'-- Not Moved Molecular Coordinates Between: '
        end if  
        write(313,'(1x,"-->",i4,"-*-",2x,a80)')Mol,Mol_Set(Mol)%Identity
        write(313,'(1x,"-->",i4,"-*-",2x,a80)')Mom,Mol_Set(Mom)%Identity
        WRITE(313,*)'------------------------------------------------'
!
        do ia=1,Mol_Set(Mol)%N_atoms
        write(313,*)Mol_Set(Mol)%Z(ia),(xyzta(ia,i),i=1,3)
        END do
        WRITE(313,*)'------------------------------------------------'
!
        do ib=1,Mol_Set(Mom)%N_atoms
        write(313,*)Mol_Set(Mom)%Z(ib),(xyztb(ib,i),i=1,3) 
        END do
        WRITE(313,*)'------------------------------------------------'
!
        end if !Zenith
!
!    	ZZ(molmom)=Sobmax
!    	CRB_INDX(molmom)=rpm
!
        Molchk=Molchk+1
        if(Molchk/=molmom)then
          Write(*,*)'Internal Inconsistency ! Better stop... enter 1'
          write(*,*)Molchk,molmom
          read(*,*)Icheck
          if(Icheck==1)Stop
          end if
		write(979,*)molmom,mol,mom,Sobmax,rpm
!
        END do  !Mom
        end do  !Mol
!
        close(313)
        close(777)
        close(979)
!
!		create a file with the file name where information is gathered for further file building
!
	  open(777,File='Name_Bound.txt')
      write(777,'(a70)')AUXILY
      close(777)        
!
		allocate(CRB(MM),ZZ(MM))
!
      open(979,File=AUXILY)
!
!		retrieve both QSM and Carbó Index matrices
!
      Do iii=1,MM
        read(979,*,END=321)molmom,mol,mom,Sobmax,rpm
        If(iii/=molmom)then
        write(*,*)' Internal Inconsistency! Better to stop and look!'
        write(*,*)' To stop enter 1'
        write(*,*)iii,molmom
        read(*,*)Icheck
        if(Icheck==1)then; Write(*,*)' Sorry!'; stop; end if
        end if
        ZZ(iii)=Sobmax
        CRB(iii)=rpm
321     end do
        close(979) 
!      
	  Open(111,File=QSM_OUT)
      WRITE(111,'(a80)')General_title
      WRITE(111,*)MoloM,MM
!
 	  write(111,*)(I,ZZ(I),I=1,MM)
!      
      Triangle_Sum=zero
      do I=1,MoloM-1; do J=I+1,MoloM
      MolMom=I+(J*J-J)/2; Triangle_Sum=Triangle_Sum+ZZ(MolMom)
      end do; end do
      Triangle_Sum=Triangle_Sum*2
      write(111,*)' ********>'
      write(111,*)' Z Matrix Triangle_Sum : ',Triangle_Sum
      CLOSE(111)
      deallocate(ZZ)         
!	  
	  OPEN(222,FILE=CRBINDX_OUT)
      WRITE(222,'(a80)')General_title
      WRITE(222,*)MoloM,MM
!
 	  WRITE(222,*)(I,CRB(I),I=1,MM)
!
      Triangle_Sum=zero
      do I=1,MoloM-1; do J=I+1,MoloM
      MolMom=I+(J*J-J)/2; Triangle_Sum=Triangle_Sum+CRB(MolMom)
      end do; end do
      Triangle_Sum=Triangle_Sum*2
      write(222,*)' ********>'
      write(222,*)' Carbó Index Matrix Triangle_Sum : ',Triangle_Sum
   	  CLOSE(222)
!
		deallocate(CRB)
!
!--------------------------------------------------------------------------------------------------------------------
   WRITE(*,*)' Molecular Quantum Similarity Suite'
   write(*,*)' Paired Molecular Quantum Similarity Success!'
   write(*,*)'**********************************************************************'
   write(*,*)'**Do not forget:'
   write(*,*)'  Name_Bound.txt contains the given name to the MQS and CBI Matrices**'
   write(*,*)'**********************************************************************'
   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=70)  	!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
!-----------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------
!__________________________________________________________________________