subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,level_shift,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 logical,intent(in) :: level_shift 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 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,level_shift,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 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,level_shift,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