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

rm LZ

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-16 22:08:33 +01:00
parent f1b681b90a
commit 0a2211fa4c
14 changed files with 0 additions and 1008 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
src/eDFT/MFL20_lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,82 +0,0 @@
subroutine MFL20_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight,rho,EcLZ)
! Compute Marut-Fromager-Loos's LDA 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,nspin)
! Local variables
logical :: LDA_centered = .false.
integer :: iEns
double precision :: EcLZLDA(nsp)
double precision,allocatable :: aMFL(:,:)
double precision,allocatable :: EcLZeLDA(:,:)
! Output variables
double precision,intent(out) :: EcLZ(nsp)
! Allocation
allocate(aMFL(3,nEns),EcLZeLDA(nsp,nEns))
! Parameters for weight-dependent LDA correlation functional
aMFL(1,1) = -0.0238184d0
aMFL(2,1) = +0.00575719d0
aMFL(3,1) = +0.0830576d0
aMFL(1,2) = -0.0282814d0
aMFL(2,2) = +0.00340758d0
aMFL(3,2) = +0.0663967d0
aMFL(1,3) = -0.0144633d0
aMFL(2,3) = -0.0504501d0
aMFL(3,3) = +0.0331287d0
! Compute correlation energy for ground, singly-excited and doubly-excited states
do iEns=1,nEns
call elda_correlation_Levy_Zahariev_shift(nEns,aMFL(:,iEns),nGrid,weight(:),rho(:,:),EcLZeLDA(:,iEns))
end do
! LDA-centered functional
EcLZLDA(:) = 0d0
if(LDA_centered) then
call VWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:,:),EcLZLDA(:))
do iEns=1,nEns
EcLZeLDA(:,iEns) = EcLZeLDA(:,iEns) + EcLZLDA(:)- EcLZeLDA(:,1)
end do
end if
! Weight-denpendent functional for ensembles
EcLZ(:) = 0d0
do iEns=1,nEns
EcLZ(:) = EcLZ(:) + wEns(iEns)*EcLZeLDA(:,iEns)
enddo
end subroutine MFL20_lda_correlation_Levy_Zahariev_shift

70
src/eDFT/MFL20_lda_exchange_Levy_Zahariev_shift.f90
View File

@ -1,70 +0,0 @@
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

78
src/eDFT/RMFL20_lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,78 +0,0 @@
subroutine RMFL20_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight,rho,EcLZ)
! Compute the restricted Marut-Fromager-Loos LDA correlation 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 :: EcLZLDA
double precision,allocatable :: aMFL(:,:)
double precision,allocatable :: EcLZeLDA(:)
! Output variables
double precision,intent(out) :: EcLZ
! Allocation
allocate(aMFL(3,nEns),EcLZeLDA(nEns))
! Parameters for weight-dependent LDA correlation functional
aMFL(1,1) = -0.0238184d0
aMFL(2,1) = +0.00540994d0
aMFL(3,1) = +0.0830766d0
aMFL(1,2) = -0.0144633d0
aMFL(2,2) = -0.0506019d0
aMFL(3,2) = +0.0331417d0
! Compute correlation energy for ground, singly-excited and doubly-excited states
do iEns=1,nEns
call restricted_elda_correlation_Levy_Zahariev_shift(nEns,aMFL(:,iEns),nGrid,weight(:),rho(:),EcLZeLDA(iEns))
end do
! LDA-centered functional
EcLZLDA = 0d0
if(LDA_centered) then
call RVWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:),EcLZLDA)
do iEns=1,nEns
EcLZeLDA(iEns) = EcLZeLDA(iEns) + EcLZLDA - EcLZeLDA(1)
end do
end if
! Weight-denpendent functional for ensembles
EcLZ = 0d0
do iEns=1,nEns
EcLZ = EcLZ + wEns(iEns)*EcLZeLDA(iEns)
enddo
end subroutine RMFL20_lda_correlation_Levy_Zahariev_shift

74
src/eDFT/RVWN5_lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,74 +0,0 @@
subroutine RVWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight,rho,EcLZ)
! Compute the restricted VNW5 LDA correlation contribution to Levy-Zahariev shift
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: r,rs,x
double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
double precision :: dxdrs,dxdx_p,decdx_p
double precision :: drsdra,decdra_p
double precision :: ec_p
! Output variables
double precision,intent(out) :: EcLZ
! Parameters of the functional
a_p = +0.0621814D0/2D0
x0_p = -0.10498d0
b_p = +3.72744d0
c_p = +12.9352d0
! Initialization
EcLZ = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs)
x_p = x*x + b_p*x + c_p
xx0_p = x0_p*x0_p + b_p*x0_p + c_p
q_p = sqrt(4d0*c_p - b_p*b_p)
ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) &
- b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_p = 2d0*x + b_p
decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
- b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
decdra_p = drsdra*dxdrs*decdx_p
EcLZ = EcLZ + weight(iG)*decdra_p*r*r
end if
end do
end subroutine RVWN5_lda_correlation_Levy_Zahariev_shift

199
src/eDFT/VWN5_lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,199 +0,0 @@
subroutine VWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight,rho,EcLZ)
! Compute Wigner's LDA contribution to Levy-Zahariev shift
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: iG
double precision :: ra,rb,r,rs,x,z
double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
double precision :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
double precision :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
double precision :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
double precision :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
double precision :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
double precision :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
! Output variables
double precision,intent(out) :: EcLZ(nsp)
! Parameters of the functional
a_p = +0.0621814D0/2D0
x0_p = -0.10498d0
b_p = +3.72744d0
c_p = +12.9352d0
a_f = +0.0621814D0/4D0
x0_f = -0.325d0
b_f = +7.06042d0
c_f = +18.0578d0
a_a = -1d0/(6d0*pi**2)
x0_a = -0.0047584D0
b_a = +1.13107d0
c_a = +13.0045d0
! Initialization
EcLZ(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
! spin-up contribution
if(ra > threshold) then
r = ra + rb
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r
x = sqrt(rs)
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
x_p = x*x + b_p*x + c_p
x_f = x*x + b_f*x + c_f
x_a = x*x + b_a*x + c_a
xx0_p = x0_p*x0_p + b_p*x0_p + c_p
xx0_f = x0_f*x0_f + b_f*x0_f + c_f
xx0_a = x0_a*x0_a + b_a*x0_a + c_a
q_p = sqrt(4d0*c_p - b_p*b_p)
q_f = sqrt(4d0*c_f - b_f*b_f)
q_a = sqrt(4d0*c_a - b_a*b_a)
ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) &
- b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
- b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) &
- b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
dzdra = (1d0 - z)/r
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
dfzdra = dzdra*dfzdz
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_p = 2d0*x + b_p
dxdx_f = 2d0*x + b_f
dxdx_a = 2d0*x + b_a
decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
- b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
- b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
- b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
decdra_p = drsdra*dxdrs*decdx_p
decdra_f = drsdra*dxdrs*decdx_f
decdra_a = drsdra*dxdrs*decdx_a
decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+ (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
EcLZ(2) = EcLZ(2) + weight(iG)*decdra*r*r
end if
! spin-down contribution
if(rb > threshold) then
r = ra + rb
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r
x = sqrt(rs)
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
x_p = x*x + b_p*x + c_p
x_f = x*x + b_f*x + c_f
x_a = x*x + b_a*x + c_a
xx0_p = x0_p*x0_p + b_p*x0_p + c_p
xx0_f = x0_f*x0_f + b_f*x0_f + c_f
xx0_a = x0_a*x0_a + b_a*x0_a + c_a
q_p = sqrt(4d0*c_p - b_p*b_p)
q_f = sqrt(4d0*c_f - b_f*b_f)
q_a = sqrt(4d0*c_a - b_a*b_a)
ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) &
- b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
- b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) &
- b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
dzdrb = -(1d0 + z)/r
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
dfzdrb = dzdrb*dfzdz
drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_p = 2d0*x + b_p
dxdx_f = 2d0*x + b_f
dxdx_a = 2d0*x + b_a
decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
- b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
- b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
- b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
decdrb_p = drsdrb*dxdrs*decdx_p
decdrb_f = drsdrb*dxdrs*decdx_f
decdrb_a = drsdrb*dxdrs*decdx_a
decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+ (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
EcLZ(2) = EcLZ(2) + weight(iG)*decdrb*r*r
end if
end do
end subroutine VWN5_lda_correlation_Levy_Zahariev_shift

56
src/eDFT/W38_lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,56 +0,0 @@
subroutine W38_lda_correlation_Levy_Zahariev_shift(nGrid,weight,rho,EcLZ)
! Compute Wigner's LDA contribution to Levy-Zahariev shift
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: iG
double precision :: ra,rb,r
double precision :: a,d,epsc
double precision :: dFcdra,dFcdrb
! Output variables
double precision,intent(out) :: EcLZ(nsp)
! Coefficients for Wigner's LDA correlation
a = 0.04918d0
d = 0.349d0
! Compute LDA correlation matrix in the AO basis
EcLZ(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
if(r > threshold) then
epsc = ra*rb/(r + d*r**(2d0/3d0))
dFcdra = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 2d0/r + 1d0/ra)
dFcdrb = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 2d0/r + 1d0/rb)
EcLZ(2) = EcLZ(2) - weight(iG)*r*r*(dFcdra + dFcdrb)
endif
enddo
EcLZ(:) = -4d0*a*EcLZ(:)
end subroutine W38_lda_correlation_Levy_Zahariev_shift

69
src/eDFT/correlation_Levy_Zahariev_shift.f90
View File

@ -1,69 +0,0 @@
subroutine correlation_Levy_Zahariev_shift(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,EcLZ)
! Compute the correlation 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,nspin)
double precision,intent(in) :: drho(ncart,nGrid,nspin)
! Local variables
double precision :: EcLZLDA(nsp),EcLZGGA(nsp)
double precision :: aC
! Output variables
double precision,intent(out) :: EcLZ(nsp)
select case (rung)
! Hartree calculation
case(0)
EcLZ(:) = 0d0
! LDA functionals
case(1)
call lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),EcLZ(:))
! 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
aC = 0.81d0
EcLZ(:) = EcLZLDA(:) + aC*(EcLZGGA(:) - EcLZLDA(:))
! Hartree-Fock calculation
case(666)
EcLZ(:) = 0d0
end select
end subroutine correlation_Levy_Zahariev_shift

62
src/eDFT/elda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,62 +0,0 @@
subroutine elda_correlation_Levy_Zahariev_shift(nEns,aLF,nGrid,weight,rho,EcLZ)
! Compute Levy-Zahariev LDA correlation shift of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: aLF(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: iG
double precision :: ra,rb,r,ec_p
double precision :: dFcdra,dFcdrb
! Output variables
double precision,intent(out) :: EcLZ
! Compute Levy-Zahariev eLDA correlation shift
EcLZ = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
if(ra > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdra = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
dFcdra = dFcdra/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdra = ec_p*dFcdra/6d0
dFcdra = ec_p + dFcdra
EcLZ = EcLZ - weight(iG)*r*r*dFcdra
end if
if(rb > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdrb = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
dFcdrb = dFcdrb/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdrb = ec_p*dFcdrb/6d0
dFcdrb = ec_p + dFcdrb
EcLZ = EcLZ - weight(iG)*r*r*dFcdrb
end if
end do
end subroutine elda_correlation_Levy_Zahariev_shift

63
src/eDFT/exchange_Levy_Zahariev_shift.f90
View File

@ -1,63 +0,0 @@
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

51
src/eDFT/lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,51 +0,0 @@
subroutine lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns,nGrid,weight,rho,EcLZ)
! Compute the lda correlation 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,nspin)
! Output variables
double precision,intent(out) :: EcLZ(nsp)
! Select correlation functional
select case (DFA)
! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
case ('W38')
call W38_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:,:),EcLZ(:))
! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
case ('VWN5')
call VWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:,:),EcLZ(:))
! Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation)
case ('MFL20')
call MFL20_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight(:),rho(:,:),EcLZ(:))
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_correlation_Levy_Zahariev_shift

42
src/eDFT/lda_exchange_Levy_Zahariev_shift.f90
View File

@ -1,42 +0,0 @@
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

69
src/eDFT/restricted_correlation_Levy_Zahariev_shift.f90
View File

@ -1,69 +0,0 @@
subroutine restricted_correlation_Levy_Zahariev_shift(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,EcLZ)
! Compute the correlation 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
double precision :: EcLZLDA,EcLZGGA
double precision :: aC
! Output variables
double precision,intent(out) :: EcLZ
select case (rung)
! Hartree calculation
case(0)
EcLZ = 0d0
! LDA functionals
case(1)
call restricted_lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns(:),nGrid,weight(:),rho(:),EcLZ)
! 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
aC = 0.81d0
EcLZ = EcLZLDA + aC*(EcLZGGA - EcLZLDA)
! Hartree-Fock calculation
case(666)
EcLZ = 0d0
end select
end subroutine restricted_correlation_Levy_Zahariev_shift

48
src/eDFT/restricted_elda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,48 +0,0 @@
subroutine restricted_elda_correlation_Levy_Zahariev_shift(nEns,aLF,nGrid,weight,rho,EcLZ)
! Compute the restricted Levy-Zahariev LDA correlation shift of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: aLF(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: r,ec_p
double precision :: dFcdr
! Output variables
double precision,intent(out) :: EcLZ
! Compute Levy-Zahariev eLDA correlation shift
EcLZ = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdr = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
dFcdr = dFcdr/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdr = ec_p*dFcdr/6d0
dFcdr = ec_p + dFcdr
EcLZ = EcLZ - weight(iG)*r*r*dFcdr
end if
end do
end subroutine restricted_elda_correlation_Levy_Zahariev_shift

45
src/eDFT/restricted_lda_correlation_Levy_Zahariev_shift.f90
View File

@ -1,45 +0,0 @@
subroutine restricted_lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns,nGrid,weight,rho,EcLZ)
! Compute the lda correlation 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) :: EcLZ
! Select correlation functional
select case (DFA)
! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
case ('RVWN5')
call RVWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:),EcLZ)
! Marut-Fromager-Loos weight-dependent correlation functional
case ('RMFL20')
call RMFL20_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight(:),rho(:),EcLZ)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine restricted_lda_correlation_Levy_Zahariev_shift