!---------------------------------------------------------------------------------------------------
! Stereographic projections of density functions computed from ASA approximation.
! Data taken from MQSS Pack of programs output.
! Reads the first molecule found in MQSS output file.
! Output files: stereop.out & sp.dat. Last one is the set of projected (X,Y,Z) points.
!---------------------------------------------------------------------------------------------------
program STEREOGRAPHIC_PROJECTION_MQSS
implicit none

integer MaxASA ! Maximal number of ASA functions per atom
double precision pi,deg_to_rad
double precision AU_to_A

parameter (MaxASA=9)
parameter (Pi=3.1415926535897932385d0)
parameter (deg_to_rad=Pi/180)
parameter (AU_to_A=0.5291772108)

character (len=2) :: Symbol(103) ! Atomic symbols
data Symbol /"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",&
             "Cs","Ba","La", &
                  "Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu", &
                            "Hf","Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn", &
             "Fr","Ra","Ac", &
                            "Th","Pa","U ","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr"/

! Molecule
integer n ! Number of atoms
double precision, allocatable :: A(:,:)   ! Centers
integer, allocatable :: AtNum(:) ! Atomic Number
character (len=2) , allocatable :: AtSymb(:) ! Atomic Symbol
integer, allocatable :: nASA(:) ! Number of ASA functions per atom
double precision, allocatable :: ASACoeff(:,:),ASAAlpha(:,:) ! ASA parameters
double precision, allocatable :: Charge(:) ! (Mulliken) Charge
double precision com(3) ! Center of mass

character (len=99) MQSS_file ! Input data file from MQSS Programs Pack
character (len=99) molecule_name

! Sphere definition
double precision SCenter(3)
double precision Radius
logical centerCom

! Angles, coordinates & functions
double precision theta,phi,st,ct,sp,cp
double precision x,y,z
double precision Xp,Yp,Xi,Xf,Yi,Yf
double precision MinXp,MaxXp,MinYp,MaxYp
double precision rho
double precision MaxDen,MinDen
integer ScaleFactor ! Integer [1,2,5,10,...] Scale factor to increase resolution (may generate huge files!)

! Molecular file directory
character (len=99) InputDir

! Auxiliar variables
integer i,j,k
double precision suma,tpc
double precision Rho_from_ASA,EPotential_from_ASA
logical found
character (len=999) textline

!=========================================================

open(unit=1,file="stereop.out",status="unknown")
write(1,*) "****************************************************"
write(*,*) "****************************************************"
write(1,*) "  STEREOGRAPHIC PROJECTION with ASA from MQSS Pack"
write(*,*) "  STEREOGRAPHIC PROJECTION with ASA from MQSS Pack"
write(1,*) "****************************************************"
write(*,*) "****************************************************"
write(1,*)
write(*,*)

! Input file
open(unit=2,file="stereop.in",status="old")
   read(2,*) InputDir
   read(2,*) MQSS_file
   MQSS_file=TRIM(InputDir)//"\"//TRIM(MQSS_file)

   read(2,*)
   read(2,*) ! Sphere definition
   read(2,*) SCenter(:)
   read(2,*) Radius
   if (Radius<0.0d0) then
      centerCom=.true.
      Radius=-Radius
   else
      centerCom=.false.
   end if

   read(2,*)
   read(2,*) ! Plane limits in file (in sphere radius units)
   read(2,*) Xi,Xf ! for X
   read(2,*) Yi,Yf ! for Y
   read(2,*) ScaleFactor

close(2)

write(1,'(" MQSS Input file: ''",a,"''")') TRIM(MQSS_file)
write(*,'(" MQSS Input file: ''",a,"''")') TRIM(MQSS_file)
write(1,*)
write(*,*)

! Reads data from MQSS file
inquire(file=MQSS_file,exist=found)
if (.not. found) then
   write(*,*) "MQSS Input file not found!"
   write(1,*) "MQSS Input file not found!"
   close(1)
   stop
end if
open(unit=2,file=MQSS_file,status="old")

   ! Molecule name. IT ONLY READS THE FIRST ONE!
   read(2,'(a99)') molecule_name
   write(*,'(" Molecule name: ''",a,"''")') TRIM(molecule_name)
   write(1,'(" Molecule name: ''",a,"''")') TRIM(molecule_name)
   
   read(2,'(a999)') textline ! ASA basis set
   write(*,*) TRIM(textline)
   write(1,*) TRIM(textline)
   write(*,*)
   write(1,*)

   read(2,*)
   read(2,*)

   read(2,'(a999)') textline ! Name
   write(*,*) TRIM(textline)
   write(1,*) TRIM(textline)
   write(*,*)
   write(1,*)

   read(2,*) n ! Number of atoms
   write(1,'(" Number of atoms:",i3)') n
   write(*,'(" Number of atoms:",i3)') n
   write(1,*)
   write(*,*)

   read(2,*) ! Number of electrons

   ! Allocates data
   allocate(A(3,n))    ! Centers
   allocate(AtNum(n))  ! Atomic Number
   allocate(AtSymb(n)) ! Atomic Symbol
   allocate(nASA(n))   ! Number of ASA functions per atom
   allocate(ASACoeff(MaxASA,n),ASAAlpha(MaxASA,n)) ! ASA parameters
   allocate(Charge(n)) ! (Mulliken) Charge

   ! Iterates in file in order to read atomic properties & data
   do i=1,n
      read(2,*) j        ! Dummy: atom number
      read(2,*) AtNum(i) ! Atomic number
      read(2,*) Charge(i)! Atomic charge
      read(2,*) A(:,i)   ! Atomic coordinates (ASSUMED TO BE A.U.)
      read(2,*) nASA(i)  ! Number of ASA functions
      if (nASA(i)>MaxASA) then
         write(1,*) "Maximal number of ASA functions per atoms has been exceeded!"
         write(*,*) "Maximal number of ASA functions per atoms has been exceeded!"
         close(2)
         close(1)
         stop
      end if
      read(2,*) ASAAlpha(1:nASA(i),i)
      read(2,*) ASACoeff(1:nASA(i),i)
   end do

close(2)

! Symbols
do i=1,n
   if (AtNum(i)<1 .or. AtNum(i)>103) then
      write(1,*) "Invalid atomic number!"
      write(*,*) "Invalid atomic number!"
      close(1)
      stop
   end if
   AtSymb(i)=Symbol(AtNum(i))
end do

! Geometrical Center (IT IS NOT ESTRICTLY a Center of mass)
do i=1,3
   com(i)=SUM(A(i,:))/n
end do
write(1,'(" ''Geometrical Center'': (",2(f7.3,","),f7.3,")")') com(:)
write(*,'(" ''Geometrical Center'': (",2(f7.3,","),f7.3,")")') com(:)
write(1,*)
write(*,*)

! Sets sphere center at the Geometrical Center of the molecule
if (centerCom) then
   write(1,*) "The sphere center is placed at the molecular geometrical center."
   write(*,*) "The sphere center is placed at the molecular geometrical center."
   SCenter(:)=com(:)
end if

write(1,'(" Sphere Radius =",f8.4,"   Original center = (",2(f8.4," ,"),f8.4,")")') Radius,SCenter(:)
write(*,'(" Sphere Radius =",f8.4,"   Original center = (",2(f8.4," ,"),f8.4,")")') Radius,SCenter(:)
write(1,*) "This center defines now the new coordinates origin."
write(*,*) "This center defines now the new coordinates origin."
write(1,*)
write(*,*)

write(1,'(" Plane limits in file: X ->",2g14.6)') Xi,Xf
write(*,'(" Plane limits in file: X ->",2g14.6)') Xi,Xf
write(1,'("    (output data)      Y ->",2g14.6)') Yi,Yf
write(*,'("    (output data)      Y ->",2g14.6)') Yi,Yf
write(1,'(" Scale factor:",i3)') ScaleFactor
write(*,'(" Scale factor:",i3)') ScaleFactor
write(1,*)
write(*,*)

! Creates simple XYZ files
open(unit=9,file="molecule.xyz",status="unknown")
   write(9,*) n
   write(9,*) "XYZ file"
   do i=1,n
      write(9,'(" ",a2,3f7.3)') AtSymb(i),A(:,i)*AU_to_A
   end do
close(9)
open(unit=9,file="molecule_sphere.xyz",status="unknown")
   write(9,*) n+1
   write(9,*) "XYZ file"
   do i=1,n
      write(9,'(" ",a2,3f7.3)') AtSymb(i),A(:,i)*AU_to_A
   end do
   write(9,'(" ",a2,3f7.3)') "X ",SCenter(:)*AU_to_A
close(9)

write(1,'(" Mulliken charges & Original coordinates (AU):",)') 
suma=0.0d0
do i=1,n
   write(1,'(" ",i6,i3,f9.5,1x,a2," (",2(f7.3,","),f7.3,")")') i,AtNum(i),Charge(i),AtSymb(i),A(:,i)
   suma=suma+Charge(i)
end do
write(*,'("  Mulliken Charges:",f9.2," vs.",i3,": Total net charge:",f9.2)') suma,SUM(AtNum(:)),1.0d0*SUM(AtNum(:))-suma
write(1,'(" ",9x,f9.2," vs.",i5,": Total net charge:",f9.2)') suma,SUM(AtNum(:)),1.0d0*SUM(AtNum(:))-suma
write(1,*)

write(1,'(" Sum of original ASA coefficients:",)') 
suma=0.0d0
do i=1,n
   write(1,'(" ",i6,i3,1x,a2,f9.5)') i,AtNum(i),AtSymb(i),SUM(ASACoeff(1:nASA(i),i))
   suma=suma+SUM(ASACoeff(1:nASA(i),i))
end do
write(*,'("           Charges:",f9.5," vs.",i3,": Total net charge:",f9.5)') suma,SUM(AtNum(:)),1.0d0*SUM(AtNum(:))-suma
write(1,'(" ",12x,f9.5," vs.",i3,": Total net charge:",f9.5)') suma,SUM(AtNum(:)),1.0d0*SUM(AtNum(:))-suma
write(1,*)
write(*,*)
if (abs(1.0d0*SUM(AtNum(:))-suma)>0.0010d0) then
   write(1,*) "Warning! Charged molecule or invalid ASA coefficients!"
   write(*,*) "Warning! Charged molecule or invalid ASA coefficients!"
end if

! ASA coefficients transformation
write(1,*) "ASA coefficients transformation:"
do i=1,n
   write(1,'(" Atom #",i3,"(",a2,"):")') i,AtSymb(i)
   do j=1,nASA(i)
      suma=ASACoeff(j,i)
      ASACoeff(j,i)=ASACoeff(j,i)*(ASAAlpha(j,i)/Pi)**1.50d0
      write(1,'(" ",i8,f9.5," -->",f9.5)') j,suma,ASACoeff(j,i)
   end do
   write(1,*)
end do

! ASA Integration Test
write(1,'(" ASA integration test:",)')
do i=1,n
   suma=0.0d0
   do j=1,nASA(i)
      suma=suma+ASACoeff(j,i)*(Pi/ASAAlpha(j,i))**1.50d0
   end do
   write(1,'("   #",i3," Mulliken Charge=",f9.5,":",f9.5,1x,a2)') i,Charge(i),suma,AtSymb(i)
   if (abs(suma-Charge(i))>0.0010d0) then
      write(1,*) "WARNING! An inconsistence has been found in ASA normalization!"
      write(*,*) "WARNING! An inconsistence has been found in ASA normalization!"
   end if
end do
write(1,*)

! Sets new origin: the center of the sphere. SCenter(:) will become the (0,0,0) point
do i=1,n
   A(:,i)=A(:,i)-SCenter(:)
end do

! New coordinates & total number of Gaussians
write(1,*) "New coordinates (Sphere center at Origin) & total number of Gaussians:"
write(*,*) "New coordinates (Sphere center at Origin) & total number of Gaussians:"
j=0
do i=1,n
   j=j+nASA(i)
end do
write(1,'(" Total number of Gaussians:",i3)') j
write(*,'(" Total number of Gaussians:",i3)') j
do i=1,n
   write(1,'(" ",i6," :  A = (",2(f8.4," ,"),f8.4,")  Z=",i2," (",a2,")  Q=",f6.2,"  Net Q=",f6.2)') i,A(:,i),AtNum(i),AtSymb(i),Charge(i),1.0d0*AtNum(i)-Charge(i)
   write(*,'(" ",i6," :  A = (",2(f8.4," ,"),f8.4,")  Z=",i2," (",a2,")  Q=",f6.2,"  Net Q=",f6.2)') i,A(:,i),AtNum(i),AtSymb(i),Charge(i),1.0d0*AtNum(i)-Charge(i)
   do j=1,nASA(i)
      write(1,'(" ",i6," Coeff =",f12.4,"  Alpha =",f12.4)') j,ASACoeff(j,i),ASAAlpha(j,i)
   end do
   write(1,*)
end do
write(1,*)
write(*,*)

! Distances matrix
write(1,*) "Distances matrix (AU):"
write(1,'(" ",8x,999(6x,a2,6x))') AtSymb(:)
do i=1,n
   write(1,'(" ",i5,1x,a2,999g14.6)') i,AtSymb(i),(dsqrt(SUM((A(:,i)-A(:,j))**2)),j=1,n)
end do
write(1,*)
write(1,*) "Distances matrix (A):"
write(1,'(" ",8x,999(6x,a2,6x))') AtSymb(:)
do i=1,n
   write(1,'(" ",i5,1x,a2,999g14.6)') i,AtSymb(i),(dsqrt(SUM((A(:,i)-A(:,j))**2))*AU_to_A,j=1,n)
end do
write(1,*)

write(1,*) "Distances from sphere center (AU):"
do i=1,n
   write(1,'(" ",i5,1x,a2,g14.6)') i,AtSymb(i),dsqrt(SUM((A(:,i))**2))
end do
write(1,*)
write(1,*) "Distances from sphere center (A):"
do i=1,n
   write(1,'(" ",i5,1x,a2,g14.6)') i,AtSymb(i),dsqrt(SUM((A(:,i))**2))*AU_to_A
end do
write(1,*)

! ============================ Stereographic Projection starts here ======================================

open(unit=7,file="sp.dat",status="unknown")

MinXp=+9.99d99
MaxXp=-9.99d99
MinYp=+9.99d99
MaxYp=-9.99d99

MinDen=+9.99d99
MaxDen=-9.99d99

write(*,*)
write(*,'("+",f7.1,"%")') 0.0
tpc=100.0/360/ScaleFactor

do i=0,360*ScaleFactor ! Loop over points in the sphere (a fraction of Degree)

   if (mod(i,10)==0) write(*,'("+",f7.1)') tpc*i

   theta=i*deg_to_rad/ScaleFactor
   st=sin(theta)
   ct=cos(theta)

   do j=0,180*ScaleFactor
      phi=j*deg_to_rad/ScaleFactor
      sp=sin(phi)
      cp=cos(phi)

      ! Coordinates on the sphere surface (origin at its center)
      x=Radius*st*cp
      y=Radius*st*sp
      z=Radius*ct

      ! Projection in plane
      if (abs(Radius-z)>1.0d-8) then
         Xp=x/(Radius-z)
         Yp=y/(Radius-z)

         MinXp=min(MinXp,Xp)
         MaxXp=max(MaxXp,Xp)
         MinYp=min(MinYp,Yp)
         MaxYp=max(MaxYp,Yp)

         ! Puts point in file if 2D cartessian limits are fullfilled
         if (Xp>=Xi) then
            if (Xp<=Xf) then
               if (Yp>=Yi) then
                  if (Yp<=Yf) then

                     ! Rho function value
                     rho=Rho_from_ASA(n,A,ASAAlpha,ASACoeff,nASA,MaxASA,x,y,z)

                     write(7,'(" ",2g14.6,g20.10)') Xp,Yp,rho

                     MinDen=Min(MinDen,rho)
                     MaxDen=Max(MaxDen,rho)

                  end if
               end if
            end if
        end if

      end if

   end do
end do

close(7)

write(1,'(" Minimal & Maximal projected X values:",2g14.6)') MinXp,MaxXp
write(*,'(" Minimal & Maximal projected X values:",2g14.6)') MinXp,MaxXp
write(1,'("                             Y values:",2g14.6)') MinYp,MaxYp
write(*,'("                             Y values:",2g14.6)') MinYp,MaxYp
write(1,*)
write(*,*)
write(1,'(" Minimal & Maximal Density values found:",2g14.6)') MinDen,MaxDen
write(*,'(" Minimal & Maximal Density values found:",2g14.6)') MinDen,MaxDen
write(1,*)
write(*,*)

write(1,*)
write(*,*)
write(1,*) "------------------------------------ Normal Termination ------------------------------------"
write(*,*) "------------------------------------ Normal Termination ------------------------------------"
close(1)
END
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
double precision function Rho_from_ASA(n,A,ASAAlpha,ASACoeff,nASA,MaxASA,x,y,z)
implicit none

integer n ! Number of atoms
integer nASA(n) ! Number of Gaussians in ASA expansion
integer MaxASA ! Maximal number of Gaussians in ASA expansion
double precision A(3,n)   ! Centers
double precision ASAAlpha(MaxASA,n) ! ASA Exponents
double precision ASACoeff(MaxASA,n) ! ASA Coefficients
double precision x,y,z    ! Point Coordinates

! Auxiliar
integer i,j
double precision r2

Rho_from_ASA=0.0d0
do i=1,n
   r2=(x-A(1,i))**2+(y-A(2,i))**2+(z-A(3,i))**2
   do j=1,nASA(i)
      Rho_from_ASA=Rho_from_ASA+ASACoeff(j,i)*dexp(-ASAAlpha(j,i)*r2)
   end do
end do

END
!---------------------------------------------------------------------------------------------------
