remove old dft grid

src/QuAcK/QuAcK.f90

@@ -204,7 +204,6 @@ program QuAcK
   allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell),DShell(maxShell,maxK), & 
            ExpShell(maxShell,maxK),max_ang_mom(nNuc),min_exponent(nNuc,maxL+1),max_exponent(nNuc))
 
-
 !------------------------------------------------------------------------
 ! Read basis set information
 !------------------------------------------------------------------------
@@ -315,8 +314,9 @@ program QuAcK
   if(doKS) then
 
     call cpu_time(start_KS)
-!   call KS(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & 
+    call eDFT(nNuc,ZNuc,rNuc,nBas,nEl,nC,nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
-!            nBas,nO,S,T,V,Hc,ERI_AO,dipole_int,X,EUHF,eHF,cHF,PHF)
+              max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int)
+             
     call cpu_time(end_KS)
 
     t_KS = end_KS - start_KS

src/eDFT/eDFT.f90

@@ -1,4 +1,5 @@
-program eDFT
+subroutine eDFT(nNuc,ZNuc,rNuc,nBas,nEl,nC,nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
+              max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int)
 
 ! exchange-correlation density-functional theory calculations
 
@@ -7,28 +8,43 @@ program eDFT
   implicit none
   include 'parameters.h'
 
-  integer                       :: nNuc,nBas
+! Input variables
-  integer                       :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
-  double precision              :: ENuc,Ew,ncent
 
-  double precision,allocatable  :: ZNuc(:),rNuc(:,:)
+  integer,intent(in)            :: nNuc
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nEl(nspin)
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  double precision,intent(in)   :: ENuc,Ew
+
+  double precision,intent(in)   :: ZNuc(nNuc)
+  double precision,intent(in)   :: rNuc(nNuc,ncart)
+
+  integer,intent(in)            :: nBas,nShell
+  double precision,intent(in)   :: CenterShell(maxShell,ncart)
+  integer,intent(in)            :: TotAngMomShell(maxShell)
+  integer,intent(in)            :: KShell(maxShell)
+  double precision,intent(in)   :: DShell(maxShell,maxK)
+  double precision,intent(in)   :: ExpShell(maxShell,maxK)
+  integer,intent(in)            :: max_ang_mom(nNuc)
+  double precision,intent(in)   :: min_exponent(nNuc,maxL+1)
+  double precision,intent(in)   :: max_exponent(nNuc)
+
+
+  double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: Hc(nBas,nBas)
+  double precision,intent(in)   :: X(nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+
+! Local variables
 
-  integer                       :: nShell
-  integer,allocatable           :: TotAngMomShell(:)
-  integer,allocatable           :: KShell(:)
-  double precision,allocatable  :: CenterShell(:,:)
-  double precision,allocatable  :: DShell(:,:)
-  double precision,allocatable  :: ExpShell(:,:)
-  integer,allocatable           :: max_ang_mom(:)
   double precision,allocatable  :: min_exponent(:,:)
   double precision,allocatable  :: max_exponent(:)
-
-  double precision,allocatable  :: S(:,:)
-  double precision,allocatable  :: T(:,:)
-  double precision,allocatable  :: V(:,:)
-  double precision,allocatable  :: Hc(:,:)
-  double precision,allocatable  :: X(:,:)
-  double precision,allocatable  :: ERI(:,:,:,:)
   double precision,allocatable  :: c(:,:)
 
   character(len=8)              :: method
@@ -77,40 +93,11 @@ program eDFT
 
 ! Libxc version
 
-! call xc_f90_version(vmajor, vminor, vmicro)
+  call xc_f90_version(vmajor, vminor, vmicro)
-! write(*,'("Libxc version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
+  write(*,'("Libxc version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
 
 ! call xcinfo()
 
-!------------------------------------------------------------------------
-! Read input information
-!------------------------------------------------------------------------
-
-! Read number of atoms, number of electrons of the system
-! nC   = number of core orbitals
-! nO   = number of occupied orbitals
-! nV   = number of virtual orbitals (see below)
-! nR   = number of Rydberg orbitals 
-! nBas = number of basis functions (see below)
-!      = nO + nV
-
-  call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
-  allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
-
-! Read geometry
-
-  call read_geometry(nNuc,ZNuc,rNuc,ENuc)
-
-  allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell),DShell(maxShell,maxK), & 
-           ExpShell(maxShell,maxK),max_ang_mom(nNuc),min_exponent(nNuc,maxL+1),max_exponent(nNuc))
-
-!------------------------------------------------------------------------
-! Read basis set information
-!------------------------------------------------------------------------
-
-  call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & 
-                  max_ang_mom,min_exponent,max_exponent)
-
 !------------------------------------------------------------------------
 ! DFT options
 !------------------------------------------------------------------------
@@ -121,45 +108,15 @@ program eDFT
   call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & 
                         maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice)
 
-!------------------------------------------------------------------------
-! Read one- and two-electron integrals
-!------------------------------------------------------------------------
-
-! Memory allocation for one- and two-electron integrals
-
-  allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas), &
-           X(nBas,nBas),ERI(nBas,nBas,nBas,nBas),c(nBas,nBas))
-
- ! Read integrals
-
-  call cpu_time(start_int)
-
-  call system('./GoQCaml')
-  call read_integrals(nBas,S,T,V,Hc,ERI)
-
-  call cpu_time(end_int)
-
-  t_int = end_int - start_int
-  write(*,*)
-  write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds'
-  write(*,*)
-
-! Orthogonalization X = S^(-1/2)
-
-  call orthogonalization_matrix(ortho_type,nBas,S,X)
-
 !------------------------------------------------------------------------
 ! Construct quadrature grid
 !------------------------------------------------------------------------
   call read_grid(SGn,radial_precision,nRad,nAng,nGrid)
-! nGrid = nRad*nAng
 
   call allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid)
 
   allocate(root(ncart,nGrid),weight(nGrid))
 
-! call quadrature_grid(nRad,nAng,nGrid,root,weight)
-
   call build_grid(nNuc,ZNuc,rNuc,max_ang_mom,min_exponent,max_exponent, & 
                   radial_precision,nRad,nAng,nGrid,weight,root)
 

src/eDFT/quadrature_grid.f90

@@ -1,79 +0,0 @@
-subroutine quadrature_grid(nRad,nAng,nGrid,root,weight)
-
-! Build roots and weights of quadrature grid
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)                 :: nRad
-  integer,intent(in)                 :: nAng
-  integer,intent(in)                 :: nGrid
-
-! Local variables
-
-  integer                            :: i,j,k
-  double precision                   :: scale
-  double precision,allocatable       :: Radius(:)
-  double precision,allocatable       :: RadWeight(:)
-  double precision,allocatable       :: XYZ(:,:)
-  double precision,allocatable       :: XYZWeight(:)
-
-! Output variables
-
-  double precision,intent(out)       :: root(ncart,nGrid)
-  double precision,intent(out)       :: weight(nGrid)
-
-! Memory allocation
-
-  allocate(Radius(nRad),RadWeight(nRad),XYZ(ncart,nAng),XYZWeight(nAng))
-
-! Findthe radial grid
-
-  scale = 1d0
-  call EulMac(Radius,RadWeight,nRad,scale)
-
-  write(*,20)
-  write(*,30)
-  write(*,20)
-  do i = 1,nRad
-     write(*,40) i,Radius(i),RadWeight(i)
-  end do
-  write(*,20)
-  write(*,*)
-
-! Find the angular grid
-
-  call Lebdev(XYZ,XYZWeight,nAng)
-
-  write(*,20)
-  write(*,50)
-  write(*,20)
-  do j = 1,nAng
-     write(*,60) j,(XYZ(k,j),k=1,ncart), XYZWeight(j)
-  end do
-  write(*,20)
-
-! Form the roots and weights
-
-  k = 0
-  do i=1,nRad
-    do j=1,nAng
-       k = k + 1
-       root(:,k) = Radius(i)*XYZ(:,j)
-       weight(k) = RadWeight(i)*XYZWeight(j)
-    enddo
-  enddo
- 
-! Print grids
-
-20 format(T2,58('-'))
-30 format(T20,'Radial Quadrature',/, &
-          T6,'I',T26,'Radius',T50,'Weight')
-40 format(T3,I4,T18,F17.10,T35,F25.10)
-50 format(T20,'Angular Quadrature',/, &
-          T6,'I',T19,'X',T29,'Y',T39,'Z',T54,'Weight')
-60 format(T3,I4,T13,3F10.5,T50,F10.5)
-
-end subroutine quadrature_grid