10
1
mirror of https://github.com/pfloos/quack synced 2024-12-24 13:23:38 +01:00
QuAcK/src/eDFT/eDFT.f90

199 lines
6.8 KiB
Fortran

subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,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,dipole_int,Ew,eKS,cKS,PKS,Vxc)
! exchange-correlation density-functional theory calculations
! use xc_f90_lib_m
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
double precision,intent(in) :: thresh
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
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
integer,intent(in) :: 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
character(len=8) :: method
integer :: x_rung,c_rung
integer :: x_DFA,c_DFA
logical :: LDA_centered = .true.
integer :: SGn
double precision :: radial_precision
integer :: nRad
integer :: nAng
integer :: nGrid
double precision,allocatable :: root(:,:)
double precision,allocatable :: weight(:)
integer :: nCC
double precision,allocatable :: aCC(:,:)
double precision,allocatable :: AO(:,:)
double precision,allocatable :: dAO(:,:,:)
double precision :: start_KS,end_KS,t_KS
double precision :: start_int,end_int,t_int
integer :: nEns
logical :: doNcentered
double precision,allocatable :: wEns(:)
double precision,allocatable :: occnum(:,:,:)
integer :: Cx_choice
integer :: i,vmajor,vminor,vmicro
integer :: iBas,iEns,ispin
integer :: icart,iGrid
! Output variables
double precision,intent(out) :: Ew
double precision,intent(out) :: eKS(nBas,nspin)
double precision,intent(out) :: cKS(nBas,nBas,nspin)
double precision,intent(out) :: PKS(nBas,nBas,nspin)
double precision,intent(out) :: Vxc(nBas,nspin)
! Hello World
write(*,*)
write(*,*) '******************************************'
write(*,*) '* eDFT: density-functional for ensembles *'
write(*,*) '******************************************'
write(*,*)
!------------------------------------------------------------------------
! DFT options
!------------------------------------------------------------------------
! Allocate ensemble weights and MO coefficients
allocate(wEns(maxEns),aCC(maxCC,maxEns-1),occnum(nBas,nspin,maxEns))
call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,nCC,aCC, &
doNcentered,occnum,Cx_choice)
!------------------------------------------------------------------------
! Construct quadrature grid
!------------------------------------------------------------------------
if(SGn == -1) then
write(*,*) '*** Quadrature grid on atomic sites ! ***'
write(*,*)
nGrid = nNuc
allocate(root(ncart,nGrid),weight(nGrid))
do icart=1,ncart
do iGrid=1,nGrid
root(icart,iGrid) = rNuc(iGrid,icart)
end do
end do
weight(:) = 1d0
else
call read_grid(SGn,radial_precision,nRad,nAng,nGrid)
call allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid)
allocate(root(ncart,nGrid),weight(nGrid))
call build_grid(nNuc,ZNuc,rNuc,max_ang_mom,min_exponent,max_exponent, &
radial_precision,nRad,nAng,nGrid,weight,root)
end if
!------------------------------------------------------------------------
! Calculate AO values at grid points
!------------------------------------------------------------------------
allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid))
call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO)
!------------------------------------------------------------------------
! Compute UKS energy
!------------------------------------------------------------------------
if(method == 'UKS') then
! Reset occupation numbers for conventional UKS calculation
occnum(:,:,:) = 0d0
do ispin=1,nspin
do iBas=1,nO(ispin)
do iEns=1,nEns
occnum(iBas,ispin,iEns) = 1d0
end do
end do
end do
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute UKS energy for ensembles
!------------------------------------------------------------------------
if(method == 'eDFT-UKS') then
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! End of eDFT
!------------------------------------------------------------------------
end subroutine eDFT