!---------------------------------------------------------------------------------------------------
! Stereographic projections of density functions computed from ASA approximation
! Data taken from MQSS Pack of programs output
!
! The program also generates a 2D Density & EMP grids in an specified Z plane:
!     -sp_rho_2D.dat           : Density 2D Map
!     -sp_potential_2D.dat     : EMP 2D Map
!     -sp_potential_2D_neg.dat : EMP Negative values (others are set to zero)
!---------------------------------------------------------------------------------------------------
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(:)
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
double precision MaxEP,MinEP
double precision EPotential ! Electrostatic Potential

! Molecular file directory
character (len=99) InputDir

! 2D GRAPH REPRESENTATIONS
double precision deltaX,deltaY
double precision XmaxD,XminD,YmaxD,YminD
double precision XmaxP,XminP,YmaxP,YminP
double precision ZPlane ! Z Plane level
double precision Margin ! Percentage of margin factor beyond minimal & maximal x-y cartesian coordinates
integer np2D ! Number of graphical points per dimension

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

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

open(unit=1,file="sp-2.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="sp-2.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,*)
   read(2,*)
   read(2,*) ZPlane ! Z Plane level
   read(2,*) Margin ! Percentage of margin factor beyond minimal & maximal x-y cartesian coordinates
   read(2,*) np2D   ! Number of graphical points per dimension

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!"
   stop
   close(1)
end if
open(unit=2,file=MQSS_file,status="old")
   ! Header
   do i=1,11
      read(2,*)
   end do

   ! Molecule name. IT ONLY READS THE FIRST ONE!  <=========
   read(2,*) molecule_name
   write(*,'(" Molecule name: ''",a,"''")') TRIM(molecule_name)
   write(1,'(" Molecule name: ''",a,"''")') TRIM(molecule_name)
   read(2,*)
   
   read(2,'(a999)') textline
   write(*,*)
   write(1,*)
   write(*,'(" CHECK THIS!! : ",a)') TRIM(textline) 
   write(1,'(" CHECK THIS!! : ",a)') TRIM(textline) 
   write(1,*)
   write(*,*)

   read(2,*)
   read(2,*)
   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))
   allocate(nASA(n)) ! Number of ASA functions per atom
   allocate(ASACoeff(MaxASA,n),ASAAlpha(MaxASA,n)) ! ASA parameters
   allocate(Charge(n)) ! Mulliken Charge

   ! Iterates 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
      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
   AtSymb(i)=Symbol(AtNum(i))
end do

! Geometrical Center (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(*,*)

! Set sphere center at the 'Center' of the molecule
if (centerCom) then
   write(1,*) "The sphere center is placed at the Center of Mass."
   write(*,*) "The sphere center is placed at the Center of Mass."
   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 the new origin."
write(*,*) "This center defines the new 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,'("                       Y ->",2g14.6)') Yi,Yf
write(*,'("                       Y ->",2g14.6)') Yi,Yf
write(1,*)
write(*,*)

write(1,'(" Z-Plane in 2D representations:",g14.6)') ZPlane ! Z Plane level
write(*,'(" Z-Plane in 2D representations:",g14.6)') ZPlane ! Z Plane level
write(1,'(" Margin factor",f7.2,"%")') Margin ! Percentage of margin factor beyond minimal & maximal x-y cartesian coordinates
write(*,'(" Margin factor",f7.2,"%")') Margin ! Percentage of margin factor beyond minimal & maximal x-y cartesian coordinates
write(1,'(" Number of graphical points per dimension:",i5)') np2D ! Number of graphical points per dimension
write(*,'(" Number of graphical points per dimension:",i5)') np2D ! Number of graphical points per dimension
write(1,*)
write(*,*)

! 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:",)') 
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.5," vs.",i3,": Total net charge:",f9.5)') suma,SUM(AtNum(:)),1.0d0*SUM(AtNum(:))-suma
write(1,'(" ",9x,f9.5," vs.",i3,": Total net charge:",f9.5)') 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! Found inconsistence in ASA normalization!"
      write(*,*) "WARNING! Found inconsistence in ASA normalization!"
   end if
end do
write(1,*)

! Sets new origin: the center of the sphere. SCenter(:) will become (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 (Center sphere = Origin) & total number of gaussians:"
write(*,*) "New coordinates (Center sphere = 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,*)

! ============================ Projection starts here ======================================
open(unit=8,file="sp_potential.dat",status="unknown")
open(unit=7,file="sp2_rho.dat",status="unknown")

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

MinEP=+9.99d99
MaxEP=-9.99d99

MinDen=+9.99d99
MaxDen=-9.99d99

do i=1,360 ! Loop over points in the sphere (Degrees)
   theta=i*deg_to_rad
   st=sin(theta)
   ct=cos(theta)
   do j=0,180
      phi=j*deg_to_rad
      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)

         ! Put 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)

                     ! Electrostatic potential function value
                     EPotential=EPotential_from_ASA(n,A,AtNum,ASAAlpha,ASACoeff,nASA,MaxASA,x,y,z)

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

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

                     MinEP=Min(MinEP,EPotential)
                     MaxEP=Max(MaxEP,EPotential)

                  end if
               end if
            end if
        end if

      end if

   end do
end do

close(8)
close(7)
! ==========================================================================================

write(1,*) "STEREOGRAPHIC PROJECTION:"
write(*,*) "STEREOGRAPHIC PROJECTION:"
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:",2g14.6)') MinDen,MaxDen
write(*,'(" Minimal & Maximal Density values:",2g14.6)') MinDen,MaxDen
write(1,'("      Minimal & Maximal EP values:",2g14.6)') MinEp,MaxEp
write(*,'("      Minimal & Maximal EP values:",2g14.6)') MinEp,MaxEp
write(1,*)
write(*,*)

! ============================ 2D Projection starts here ======================================

write(*,*) "PERFORMING 2D CALCULATIONS..."
write(1,*) "PERFORMING 2D CALCULATIONS..."

MinXp=minval(A(1,:))
MaxXp=maxval(A(1,:))
deltaX=MaxXp-MinXp
if (deltaX<0.10d0) then
   deltaX=1.0d0
   MinXp=MinXp-deltaX
   MaxXp=MaxXp+deltaX
else
   deltaX=deltaX*Margin/100
   MinXp=MinXp-deltaX
   MaxXp=MaxXp+deltaX
end if
deltaX=(MaxXp-MinXp)/np2D

MinYp=minval(A(2,:))
MaxYp=maxval(A(2,:))
deltaY=MaxYp-MinYp
if (deltaY<0.10d0) then
   deltaY=1.0d0
   MinYp=MinYp-deltaY
   MaxYp=MaxYp+deltaY
else
   deltaY=deltaY*Margin/100
   MinYp=MinYp-deltaY
   MaxYp=MaxYp+deltaY
end if
deltaY=(MaxYp-MinYp)/np2D

write(1,'(" X range: [",g14.6,",",g14.6," ]")') MinXp,MaxXp
write(*,'(" X range: [",g14.6,",",g14.6," ]")') MinXp,MaxXp
write(1,'(" Y range: [",g14.6,",",g14.6," ]")') MinYp,MaxYp
write(*,'(" Y range: [",g14.6,",",g14.6," ]")') MinYp,MaxYp

MinEP=+9.99d99
MaxEP=-9.99d99

MinDen=+9.99d99
MaxDen=-9.99d99

open(unit=7,file="sp_rho_2D.dat",status="unknown")
open(unit=8,file="sp_potential_2D.dat",status="unknown")
open(unit=9,file="sp_potential_2D_neg.dat",status="unknown")

do i=0,np2D
   x=MinXp+i*deltaX

   do j=0,np2D
      y=MinYp+j*deltaY     

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

      ! Electrostatic potential function value
      EPotential=EPotential_from_ASA(n,A,AtNum,ASAAlpha,ASACoeff,nASA,MaxASA,x,y,ZPlane)

      write(7,'(" ",2g14.6,g20.10)') x,y,rho
      write(8,'(" ",3g14.6,g20.10)') x,y,EPotential
      if (EPotential<0.0d0) then;
         write(9,'(" ",3g14.6,g20.10)') x,y,EPotential
      else
         write(9,'(" ",3g14.6,g20.10)') x,y,0.0d0
      end if

      if (MinDen>rho) then
         MinDen=rho
         XminD=x
         YminD=y
      end if
      if (MaxDen<rho) then
         MaxDen=rho
         XmaxD=x
         YmaxD=y
      end if

      if (MinEP>EPotential) then
         MinEP=EPotential
         XminP=x
         YminP=y
      end if
      if (MaxEP<EPotential) then
         MaxEP=EPotential
         XmaxP=x
         YmaxP=y
      end if

   end do
end do

close(9)
close(8)
close(7)
! ==========================================================================================

write(1,*) "2D PROJECTION:"
write(*,*) "2D PROJECTION:"
write(1,'(" Minimal & Maximal X values:",2g14.6)') MinXp,MaxXp
write(*,'(" Minimal & Maximal 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 Density value:",g14.6," X-Y Position:",2g14.6)') MinDen,XminD,YminD
write(1,'(" Maximal Density value:",g14.6," X-Y Position:",2g14.6)') MaxDen,XmaxD,YmaxD
write(*,'(" Minimal Density value:",g14.6," X-Y Position:",2g14.6)') MinDen,XminD,YminD
write(*,'(" Maximal Density value:",g14.6," X-Y Position:",2g14.6)') MaxDen,XmaxD,YmaxD
write(1,*)
write(*,*)
write(1,'(" Minimal EP value:",g14.6," X-Y Position:",2g14.6)') MinEP,XminP,YminP
write(1,'(" Maximal EP value:",g14.6," X-Y Position:",2g14.6)') MaxEP,XmaxP,YmaxP
write(*,'(" Minimal EP value:",g14.6," X-Y Position:",2g14.6)') MinEP,XminP,YminP
write(*,'(" Maximal EP value:",g14.6," X-Y Position:",2g14.6)') MaxEP,XmaxP,YmaxP
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
!---------------------------------------------------------------------------------------------------
!-!
!---------------------------------------------------------------------------------------------------
! It's assumed that the ASA integration for an atom returns its electronic charge.
!---------------------------------------------------------------------------------------------------
double precision function EPotential_from_ASA(n,A,AtNum,ASAAlpha,ASACoeff,nASA,MaxASA,x,y,z)
implicit none

double precision Pi
parameter (Pi=3.1415926535897932385d0)

integer n ! Number of Atoms
integer nASA(n) ! Number of Gaussians in ASA expansion per atom
integer MaxASA  ! Maximal number of Gaussians in ASA expansion per atom
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
integer AtNum(n) ! Atomic numbers

! Auxiliar
integer i,j
double precision r2,NP,EP
double precision F(0:36)

NP=0.0d0 ! Nuclear part
EP=0.0d0 ! Electronic part

do i=1,n
   r2=(x-A(1,i))**2+(y-A(2,i))**2+(z-A(3,i))**2
   if (r2<=1.0d-8) then
      write(*,*) "A point collapses with nucleus in nuclear electrostatic potential calculation!"
      write(1,*) "A point collapses with nucleus in nuclear electrostatic potential calculation!"
      close(1)
      close(8)
      close(7)
      stop
   else
      NP=NP+1.0d0*AtNum(i)/dsqrt(r2)
      do j=1,nASA(i)
         call Fm_integrals(0,ASAAlpha(j,i)*r2,F,1.0d0)
         EP=EP+ASACoeff(j,i)*F(0)/ASAAlpha(j,i)
      end do
   end if
end do
EP=EP*2*Pi

EPotential_from_ASA=NP-EP

END
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
include "fm.sub"
!---------------------------------------------------------------------------------------------------
