PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-05 05:33:50 +01:00
Code Releases Activity

commit b4 removing LZ routines

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-16 22:08:04 +01:00
parent 9bf424feab
commit f1b681b90a
16 changed files with 425 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/dft
View File

@ -1,7 +1,7 @@
# Restricted or unrestricted KS calculation # Restricted or unrestricted KS calculation
GOK-RKS GOK-RKS
# exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666 # exchange rung: Hartree = 0, LDA = 1 (RS51,S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666
666 HF 1 RS51
# correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666 # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666
0 HF 0 HF
# quadrature grid SG-n # quadrature grid SG-n
@ -11,4 +11,4 @@
# Ensemble weights: wEns(1),...,wEns(nEns-1) # Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00000 0.00000 0.00000 0.00000
# eKS: maxSCF thresh DIIS n_diis guess_type ortho_type # eKS: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 15 1 1 64 0.0000001 T 5 1 1

3
src/eDFT/GOK_RKS.f90
View File

@ -239,7 +239,7 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
! Build Fock operator ! Build Fock operator
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*Fx(:,:) + Fc(:,:) F(:,:) = Hc(:,:) + J(:,:) + Fx(:,:) + Fc(:,:)
! Check convergence ! Check convergence
@ -276,7 +276,6 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, &
Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex) Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex)
Ex = 0.5d0*Ex
! Correlation energy ! Correlation energy

70
src/eDFT/MFL20_lda_exchange_Levy_Zahariev_shift.f90 Normal file
View File

@ -0,0 +1,70 @@
subroutine MFL20_lda_exchange_Levy_Zahariev_shift(nEns,wEns,nGrid,weight,rho,ExLZ)
! Compute the Marut-Fromager-Loos LDA exchange contribution to Levy-Zahariev shift
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
logical :: LDA_centered = .true.
integer :: iEns
double precision :: ExLZLDA
double precision,allocatable :: aMFL(:,:)
double precision,allocatable :: ExLZeLDA(:)
! Output variables
double precision,intent(out) :: ExLZ
! Allocation
allocate(aMFL(3,nEns),ExLZeLDA(nEns))
! Parameters for weight-dependent LDA correlation functional
! Compute correlation energy for ground, singly-excited and doubly-excited states
do iEns=1,nEns
! call elda_exchange_Levy_Zahariev_shift(nEns,aMFL(:,iEns),nGrid,weight(:),rho(:),ExLZeLDA(iEns))
end do
! LDA-centered functional
ExLZLDA = 0d0
if(LDA_centered) then
! call S51_lda_exchange_Levy_Zahariev_shift(nGrid,weight(:),rho(:),ExLZLDA)
do iEns=1,nEns
ExLZeLDA(iEns) = ExLZeLDA(iEns) + ExLZLDA - ExLZeLDA(1)
end do
end if
! Weight-denpendent functional for ensembles
ExLZ = 0d0
do iEns=1,nEns
ExLZ = ExLZ + wEns(iEns)*ExLZeLDA(iEns)
enddo
end subroutine MFL20_lda_exchange_Levy_Zahariev_shift

7
src/eDFT/RS51_lda_exchange_energy.f90
View File

@ -14,7 +14,8 @@ subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Local variables ! Local variables
integer :: iG integer :: iG
double precision :: alpha,r double precision :: r
double precision :: Cx
! Output variables ! Output variables
@ -22,7 +23,7 @@ subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Cx coefficient for Slater LDA exchange ! Cx coefficient for Slater LDA exchange
alpha = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) Cx = -(3d0/4d0)*(3d0/pi)**(1d0/3d0)
! Compute LDA exchange energy ! Compute LDA exchange energy
@ -33,7 +34,7 @@ subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
if(r > threshold) then if(r > threshold) then
Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) Ex = Ex + weight(iG)*Cx*r**(4d0/3d0)
endif endif

2
src/eDFT/RVWN5_lda_correlation_individual_energy.f90
View File

@ -66,7 +66,7 @@ subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
decdra_p = drsdra*dxdrs*decdx_p decdra_p = drsdra*dxdrs*decdx_p
Ec = Ec + weight(iG)*(ec_p + decdra_p*r)*rI Ec = Ec + weight(iG)*(ec_p*rI + decdra_p*r*rI - decdra_p*r*r)
end if end if

63
src/eDFT/exchange_Levy_Zahariev_shift.f90 Normal file
View File

@ -0,0 +1,63 @@
subroutine exchange_Levy_Zahariev_shift(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,ExLZ)
! Compute the exchange part of the Levy-Zahariev shift
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
! Output variables
double precision,intent(out) :: ExLZ
select case (rung)
! Hartree calculation
case(0)
ExLZ = 0d0
! LDA functionals
case(1)
call lda_exchange_Levy_Zahariev_shift(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),ExLZ)
! GGA functionals
case(2)
call print_warning('!!! Exchange LZ shift NYI for GGAs !!!')
stop
! Hybrid functionals
case(4)
call print_warning('!!! Exchange LZ shift NYI for hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
ExLZ = 0d0
end select
end subroutine exchange_Levy_Zahariev_shift

62
src/eDFT/exchange_derivative_discontinuity.f90 Normal file
View File

@ -0,0 +1,62 @@
subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD)
! Compute the exchange part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
select case (rung)
! Hartree calculation
case(0)
ExDD(:) = 0d0
! LDA functionals
case(1)
call lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),ExDD(:))
! GGA functionals
case(2)
call print_warning('!!! exchange part of the derivative discontinuity NYI for GGAs !!!')
stop
! Hybrid functionals
case(4)
call print_warning('!!! exchange part of derivative discontinuity NYI for hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
ExDD(:) = 0d0
end select
end subroutine exchange_derivative_discontinuity

68
src/eDFT/exchange_individual_energy.f90 Normal file
View File

@ -0,0 +1,68 @@
subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex)
! Compute the exchange energy of individual states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
double precision :: ExLDA
double precision :: ExGGA
double precision :: ExHF
! Output variables
double precision,intent(out) :: Ex
select case (rung)
! Hartree calculation
case(0)
Ex = 0d0
! LDA functionals
case(1)
call lda_exchange_individual_energy(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),rho(:),Ex)
! GGA functionals
case(2)
call print_warning('!!! Individual energies NYI for GGAs !!!')
stop
! Hybrid functionals
case(4)
call print_warning('!!! Individual energies NYI for hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
! call fock_exchange_individual_energy(nEns,wEns(:),nBas,P,FxHF,ExHF)
Ex = ExHF
end select
end subroutine exchange_individual_energy

1
src/eDFT/exchange_potential.f90
View File

@ -35,6 +35,7 @@ subroutine exchange_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,ERI,AO,dAO,
select case (rung) select case (rung)
! Hartree calculation ! Hartree calculation
case(0) case(0)
Fx(:,:) = 0d0 Fx(:,:) = 0d0

2
src/eDFT/fock_exchange_potential.f90
View File

@ -31,4 +31,6 @@ subroutine fock_exchange_potential(nBas,P,ERI,Fx)
enddo enddo
enddo enddo
Fx(:,:) = 0.5d0*Fx(:,:)
end subroutine fock_exchange_potential end subroutine fock_exchange_potential

42
src/eDFT/lda_exchange_Levy_Zahariev_shift.f90 Normal file
View File

@ -0,0 +1,42 @@
subroutine lda_exchange_Levy_Zahariev_shift(DFA,nEns,wEns,nGrid,weight,rho,ExLZ)
! Compute the exchange LDA part of the Levy-Zahariev shift
implicit none
include 'parameters.h'
! Input variables
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Output variables
double precision,intent(out) :: ExLZ
! Select correlation functional
select case (DFA)
case ('S51')
call print_warning('!!! Exchange Levzy-Zahariev shift NYI for S51 !!!')
stop
case ('MFL20')
call MFL20_lda_exchange_Levy_Zahariev_shift(nEns,wEns,nGrid,weight(:),rho(:),ExLZ)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_exchange_Levy_Zahariev_shift

44
src/eDFT/lda_exchange_derivative_discontinuity.f90 Normal file
View File

@ -0,0 +1,44 @@
subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,ExDD)
! Compute the exchange LDA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Select correlation functional
select case (DFA)
case ('S51')
ExDD(:) = 0d0
case ('RMFL20')
! call MFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
ExDD(:) = 0d0
case default
call print_warning('!!! LDA exchange functional not available !!!')
stop
end select
end subroutine lda_exchange_derivative_discontinuity

12
src/eDFT/lda_exchange_energy.f90
View File

@ -24,16 +24,16 @@ subroutine lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,Ex)
! Slater's LDA correlation functional ! Slater's LDA correlation functional
case ('RS51')
call RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Restricted version of Slater's LDA correlation functional
case ('S51') case ('S51')
call S51_lda_exchange_energy(nGrid,weight,rho,Ex) call S51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Restricted version of Slater's LDA correlation functional
case ('RS51')
call RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 exchange functional ! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 exchange functional
case ('RMFL20') case ('RMFL20')

41
src/eDFT/lda_exchange_individual_energy.f90 Normal file
View File

@ -0,0 +1,41 @@
subroutine lda_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ex)
! Compute LDA exchange energy for individual states
implicit none
include 'parameters.h'
! Input variables
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
! Output variables
double precision :: Ex
! Select correlation functional
select case (DFA)
case ('S51')
! call S51_lda_exchange_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ex)
case ('MFL20')
! call MFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_exchange_individual_energy

28
src/eDFT/restricted_individual_energy.f90
View File

@ -81,22 +81,18 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Exchange energy ! Individual exchange energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do iEns=1,nEns do iEns=1,nEns
call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,iEns),ERI(:,:,:,:), & call exchange_energy_individual_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,iEns),FxHF(:,:), &
AO(:,:),dAO(:,:,:),rho(:,iEns),drho(:,:,iEns),Fx(:,:),FxHF(:,:))
call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,iEns),FxHF(:,:), &
rho(:,iEns),drho(:,:,iEns),Ex(iEns)) rho(:,iEns),drho(:,:,iEns),Ex(iEns))
end do end do
Ex(:) = 0.5d0*Ex(:)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Correlation energy ! Indivudual correlation energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do iEns=1,nEns do iEns=1,nEns
@ -110,15 +106,21 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
! Compute Levy-Zahariev shift ! Compute Levy-Zahariev shift
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call restricted_correlation_Levy_Zahariev_shift(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:), & ! call exchange_Levy_Zahariev_shift(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),ExLZ)
ExLZ,EcLZ,ExcLZ)
! call restricted_correlation_Levy_Zahariev_shift(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),EcLZ)
! ExcLZ = ExLZ + ExLZ
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Compute derivative discontinuities ! Compute derivative discontinuities
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), & call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:))
ExDD(:),EcDD(:),ExcDD(:))
call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),EcDD(:))
ExcDD(:) = ExDD(:) + EcDD(:)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Total energy ! Total energy
@ -127,7 +129,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
do iEns=1,nEns do iEns=1,nEns
Exc(iEns) = Ex(iEns) + Ec(iEns) Exc(iEns) = Ex(iEns) + Ec(iEns)
E(iEns) = ENuc + ET(iEns) + EV(iEns) + EJ(iEns) & E(iEns) = ENuc + ET(iEns) + EV(iEns) + EJ(iEns) &
+ Ex(iEns) + Ec(iEns) + ExcLZ + ExcDD(iEns) + Ex(iEns) + Ec(iEns) + ExcDD(iEns)
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -142,7 +144,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
! Dump results ! Dump results
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call print_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:),ExLZ,EcLZ,ExcLZ, & call print_restricted_individual_energy(nEns,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:),ExLZ,EcLZ,ExcLZ, &
ExDD(:),EcDD(:),ExcDD(:),E(:),Om(:)) ExDD(:),EcDD(:),ExcDD(:),E(:),Om(:))
end subroutine restricted_individual_energy end subroutine restricted_individual_energy

8
src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90
View File

@ -26,21 +26,15 @@ subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGr
select case (DFA) select case (DFA)
! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
case ('W38') case ('W38')
Ec(:) = 0d0 Ec(:) = 0d0
! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
case ('VWN5') case ('VWN5')
Ec(:) = 0d0 Ec(:) = 0d0
! Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation) case ('RMFL20')
case ('LF19')
call RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),Ec(:)) call RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),Ec(:))