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

257 lines
9.2 KiB
Fortran
Raw Normal View History

2019-03-13 11:07:31 +01:00
program eDFT
! exchange-correlation density-functional theory calculations
2019-07-09 16:17:10 +02:00
implicit none
2019-03-13 11:07:31 +01:00
include 'parameters.h'
2019-07-09 16:17:10 +02:00
integer :: nNuc,nBas
integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
double precision :: ENuc,Ew,ncent
2019-03-13 11:07:31 +01:00
2019-03-13 11:34:51 +01:00
double precision,allocatable :: ZNuc(:),rNuc(:,:)
2019-03-13 11:07:31 +01:00
integer :: nShell
integer,allocatable :: TotAngMomShell(:)
integer,allocatable :: KShell(:)
double precision,allocatable :: CenterShell(:,:)
double precision,allocatable :: DShell(:,:)
double precision,allocatable :: ExpShell(:,:)
2020-03-25 09:48:58 +01:00
integer,allocatable :: max_ang_mom(:)
double precision,allocatable :: min_exponent(:,:)
double precision,allocatable :: max_exponent(:)
2019-03-13 11:07:31 +01:00
2020-03-18 16:06:29 +01:00
double precision,allocatable :: S(:,:)
double precision,allocatable :: T(:,:)
double precision,allocatable :: V(:,:)
double precision,allocatable :: Hc(:,:)
double precision,allocatable :: X(:,:)
2019-03-13 11:07:31 +01:00
double precision,allocatable :: ERI(:,:,:,:)
2020-03-31 23:33:48 +02:00
double precision,allocatable :: c(:,:)
2019-03-13 11:07:31 +01:00
character(len=8) :: method
2019-03-13 11:07:31 +01:00
integer :: x_rung,c_rung
character(len=12) :: x_DFA ,c_DFA
logical :: LDA_centered = .true.
2020-03-25 11:25:48 +01:00
2019-03-13 11:07:31 +01:00
integer :: SGn
2020-03-25 11:25:48 +01:00
double precision :: radial_precision
integer :: nRad
integer :: nAng
integer :: nGrid
2019-03-13 11:07:31 +01:00
double precision,allocatable :: root(:,:)
double precision,allocatable :: weight(:)
double precision :: aCC_w1(3)
double precision :: aCC_w2(3)
2020-03-25 11:25:48 +01:00
2019-03-13 11:07:31 +01:00
double precision,allocatable :: AO(:,:)
double precision,allocatable :: dAO(:,:,:)
double precision :: start_KS,end_KS,t_KS
2020-03-14 23:00:44 +01:00
double precision :: start_int,end_int,t_int
2019-03-13 11:07:31 +01:00
2020-09-29 11:47:18 +02:00
integer :: nEns
logical :: doNcentered
2019-03-13 11:07:31 +01:00
double precision,allocatable :: wEns(:)
integer :: maxSCF,max_diis
double precision :: thresh
2019-07-09 16:17:10 +02:00
logical :: DIIS
integer :: guess_type
integer :: ortho_type
2019-03-13 11:07:31 +01:00
2020-08-04 12:30:52 +02:00
double precision,allocatable :: occnum(:,:,:)
integer :: Cx_choice
2019-03-13 11:07:31 +01:00
! Hello World
write(*,*)
write(*,*) '******************************************'
write(*,*) '* eDFT: density-functional for ensembles *'
write(*,*) '******************************************'
write(*,*)
!------------------------------------------------------------------------
! Read input information
!------------------------------------------------------------------------
! Read number of atoms, number of electrons of the system
2019-07-09 16:17:10 +02:00
! nC = number of core orbitals
2019-03-13 11:07:31 +01:00
! nO = number of occupied orbitals
! nV = number of virtual orbitals (see below)
2019-07-09 16:17:10 +02:00
! nR = number of Rydberg orbitals
2019-03-13 11:07:31 +01:00
! nBas = number of basis functions (see below)
! = nO + nV
2019-07-09 16:17:10 +02:00
call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
2019-03-13 11:34:51 +01:00
allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
2019-03-13 11:07:31 +01:00
! Read geometry
2019-03-13 11:34:51 +01:00
call read_geometry(nNuc,ZNuc,rNuc,ENuc)
2019-03-13 11:07:31 +01:00
2020-03-25 09:48:58 +01:00
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))
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
! Read basis set information
!------------------------------------------------------------------------
2020-03-25 09:48:58 +01:00
call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent)
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
2019-07-09 16:17:10 +02:00
! DFT options
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
2019-07-09 16:17:10 +02:00
! Allocate ensemble weights
2020-09-11 11:55:04 +02:00
allocate(wEns(maxEns),occnum(nBas,nspin,maxEns))
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)
2019-03-13 11:07:31 +01:00
2019-07-09 16:17:10 +02:00
!------------------------------------------------------------------------
! Read one- and two-electron integrals
!------------------------------------------------------------------------
2019-03-13 11:07:31 +01:00
2019-07-09 16:17:10 +02:00
! Memory allocation for one- and two-electron integrals
2019-03-13 11:07:31 +01:00
2020-03-18 16:06:29 +01:00
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas), &
2020-03-31 23:33:48 +02:00
X(nBas,nBas),ERI(nBas,nBas,nBas,nBas),c(nBas,nBas))
2019-03-13 11:07:31 +01:00
2020-03-23 17:26:12 +01:00
! Read integrals
2019-03-13 11:07:31 +01:00
2020-03-14 23:00:44 +01:00
call cpu_time(start_int)
call system('./GoQCaml')
2020-09-15 10:40:51 +02:00
call read_integrals(nBas,S,T,V,Hc,ERI)
2020-03-14 23:00:44 +01:00
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(*,*)
2019-03-13 11:07:31 +01:00
2019-07-09 16:17:10 +02:00
! Orthogonalization X = S^(-1/2)
2019-03-13 11:07:31 +01:00
2019-07-09 16:17:10 +02:00
call orthogonalization_matrix(ortho_type,nBas,S,X)
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
! Construct quadrature grid
!------------------------------------------------------------------------
2020-03-25 12:56:28 +01:00
call read_grid(SGn,radial_precision,nRad,nAng,nGrid)
2020-03-25 13:03:31 +01:00
! nGrid = nRad*nAng
2020-03-25 10:39:45 +01:00
2020-03-31 12:39:26 +02:00
call allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid)
2019-03-13 11:07:31 +01:00
allocate(root(ncart,nGrid),weight(nGrid))
2020-03-25 13:03:31 +01:00
! call quadrature_grid(nRad,nAng,nGrid,root,weight)
2020-03-23 17:26:12 +01:00
2020-03-31 12:39:26 +02:00
call build_grid(nNuc,ZNuc,rNuc,max_ang_mom,min_exponent,max_exponent, &
2020-03-25 13:03:31 +01:00
radial_precision,nRad,nAng,nGrid,weight,root)
2020-03-23 17:26:12 +01:00
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
! Calculate AO values at grid points
!------------------------------------------------------------------------
allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid))
2019-07-09 16:17:10 +02:00
call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO)
2019-03-13 11:07:31 +01:00
LDA_centered = .true.
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
2020-03-18 11:10:21 +01:00
! Compute GOK-RKS energy
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
2020-03-15 08:23:01 +01:00
if(method == 'GOK-RKS') then
2019-03-13 11:07:31 +01:00
call cpu_time(start_KS)
2020-07-04 14:32:06 +02:00
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), &
S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice)
2020-03-18 11:10:21 +01:00
call cpu_time(end_KS)
t_KS = end_KS - start_KS
2020-03-30 17:50:07 +02:00
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds'
2020-03-18 11:10:21 +01:00
write(*,*)
end if
!------------------------------------------------------------------------
2020-04-26 23:17:27 +02:00
! Compute LIM excitation energies
2020-03-18 11:10:21 +01:00
!------------------------------------------------------------------------
if(method == 'LIM-RKS') then
call cpu_time(start_KS)
2020-04-26 23:17:27 +02:00
call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
2020-09-29 09:43:01 +02:00
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
2020-03-15 08:23:01 +01:00
call cpu_time(end_KS)
t_KS = end_KS - start_KS
2020-03-18 11:10:21 +01:00
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds'
2020-03-15 08:23:01 +01:00
write(*,*)
end if
2020-04-26 23:17:27 +02:00
!------------------------------------------------------------------------
! Compute MOM excitation energies
!------------------------------------------------------------------------
if(method == 'MOM-RKS') then
call cpu_time(start_KS)
call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), &
aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
2020-09-29 09:43:01 +02:00
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered)
2020-04-26 23:17:27 +02:00
call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds'
write(*,*)
end if
2020-03-15 08:23:01 +01:00
!------------------------------------------------------------------------
2020-03-18 11:10:21 +01:00
! Compute GOK-UKS energy (BROKEN)
2020-03-15 08:23:01 +01:00
!------------------------------------------------------------------------
if(method == 'GOK-UKS') then
call cpu_time(start_KS)
call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, &
2020-09-29 09:43:01 +02:00
nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, &
Cx_choice,doNcentered)
2019-03-13 11:07:31 +01:00
call cpu_time(end_KS)
t_KS = end_KS - start_KS
2020-03-15 08:23:01 +01:00
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds'
2019-03-13 11:07:31 +01:00
write(*,*)
2020-03-15 08:23:01 +01:00
end if
!------------------------------------------------------------------------
! Compute N-centered UKS energy (UNBROKEN)
!------------------------------------------------------------------------
if(method == 'eDFT-UKS') then
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, &
2020-09-29 09:43:01 +02:00
nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered)
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
2019-03-13 11:07:31 +01:00
!------------------------------------------------------------------------
! End of eDFT
!------------------------------------------------------------------------
end program eDFT