LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

remove eDFT

This commit is contained in:
Pierre-Francois Loos 2023-07-03 14:41:46 +02:00
parent 989e0b9bd8
commit 5aff287039
102 changed files with 2 additions and 8318 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/GF/dSigmaC_GF2.f90
View File

@ -19,7 +19,9 @@ double precision function dSigmaC_GF2(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,ERI)
integer :: i,j,a,b
double precision :: eps
! Initialize
dSigmaC_GF2 = 0d0
! Occupied part of the correlation self-energy
do i=nC+1,nO

BIN
src/eDFT.tgz Normal file
View File

Binary file not shown.

101
src/eDFT/AO_values_grid.f90
View File

@ -1,101 +0,0 @@
subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
nGrid,root,AO,dAO)
! Compute values of the AOs and their derivatives with respect to the cartesian coordinates
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,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)
double precision,intent(in) :: root(ncart,nGrid)
integer,intent(in) :: nGrid
! Local variables
integer :: atot,nShellFunction,a(ncart)
integer,allocatable :: ShellFunction(:,:)
double precision :: rASq,xA,yA,zA,norm_coeff,prim
integer :: iSh,iShF,iK,iG,iBas
! Output variables
double precision,intent(out) :: AO(nBas,nGrid)
double precision,intent(out) :: dAO(ncart,nBas,nGrid)
! Initialization
iBas = 0
AO(:,:) = 0d0
dAO(:,:,:) = 0d0
!------------------------------------------------------------------------
! Loops over shells
!------------------------------------------------------------------------
do iSh=1,nShell
atot = TotAngMomShell(iSh)
nShellFunction = (atot*atot + 3*atot + 2)/2
allocate(ShellFunction(1:nShellFunction,1:3))
call generate_shell(atot,nShellFunction,ShellFunction)
do iShF=1,nShellFunction
iBas = iBas + 1
a(:) = ShellFunction(iShF,:)
do iG=1,nGrid
xA = root(1,iG) - CenterShell(iSh,1)
yA = root(2,iG) - CenterShell(iSh,2)
zA = root(3,iG) - CenterShell(iSh,3)
! Calculate distance for exponential
rASq = xA**2 + yA**2 + zA**2
!------------------------------------------------------------------------
! Loops over contraction degrees
!-------------------------------------------------------------------------
do iK=1,KShell(iSh)
! Calculate the exponential part
prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)
AO(iBas,iG) = AO(iBas,iG) + prim
prim = -2d0*ExpShell(iSh,iK)*prim
dAO(:,iBas,iG) = dAO(:,iBas,iG) + prim
enddo
dAO(1,iBas,iG) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(1,iBas,iG)
if(a(1) > 0) dAO(1,iBas,iG) = dAO(1,iBas,iG) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
dAO(2,iBas,iG) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(2,iBas,iG)
if(a(2) > 0) dAO(2,iBas,iG) = dAO(2,iBas,iG) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iBas,iG)
dAO(3,iBas,iG) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(3,iBas,iG)
if(a(3) > 0) dAO(3,iBas,iG) = dAO(3,iBas,iG) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iBas,iG)
! Calculate polynmial part
AO(iBas,iG) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
enddo
enddo
deallocate(ShellFunction)
enddo
!------------------------------------------------------------------------
! End loops over shells
!------------------------------------------------------------------------
end subroutine AO_values_grid

48
src/eDFT/B88_gga_exchange_energy.f90
View File

@ -1,48 +0,0 @@
subroutine B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute Becke's 88 GGA exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
integer :: iG
double precision :: b
double precision :: r,g,x
! Output variables
double precision :: Ex
! Coefficients for B88 GGA exchange functional
b = 0.0042d0
! Compute GGA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
x = sqrt(g)/r**(4d0/3d0)
Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - b*x**2/(1d0 + 6d0*b*x*asinh(x)))
end if
end do
end subroutine B88_gga_exchange_energy

73
src/eDFT/B88_gga_exchange_potential.f90
View File

@ -1,73 +0,0 @@
subroutine B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute Becke's GGA exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: b
double precision :: vAO,gAO
double precision :: r,g,x,dxdr,dxdg,f
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Coefficients for B88 GGA exchange functional
b = 0.0042d0
! Compute GGA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
x = sqrt(g)/r**(4d0/3d0)
dxdr = - 4d0*sqrt(g)/(3d0*r**(7d0/3d0))/x
dxdg = + 1d0/(2d0*sqrt(g)*r**(4d0/3d0))/x
f = b*x**2/(1d0 + 6d0*b*x*asinh(x))
Fx(mu,nu) = Fx(mu,nu) + vAO*( &
4d0/3d0*r**(1d0/3d0)*(CxLSDA - f) &
- 2d0*r**(4d0/3d0)*dxdr*f &
+ r**(4d0/3d0)*dxdr*(6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) )
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*r**(4d0/3d0)*dxdg*( &
- 2d0*f + (6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) )
end if
end do
end do
end do
end subroutine B88_gga_exchange_potential

93
src/eDFT/C16_lda_correlation_energy.f90
View File

@ -1,93 +0,0 @@
subroutine C16_lda_correlation_energy(nGrid,weight,rho,Ec)
! Compute unrestricted Chachiyo's LDA correlation energy
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
double precision :: a_p,b_p,ec_p
double precision :: a_f,b_f,ec_f
double precision :: z,fz,ec_z
! Output variables
double precision :: Ec(nsp)
! Coefficients for Chachiyo's LDA correlation
a_p = (log(2d0) - 1d0)/(2d0*pi**2)
b_p = 20.4562557d0
a_f = (log(2d0) - 1d0)/(4d0*pi**2)
b_f = 27.4203609d0
! Compute LDA correlation energy
Ec(:) = 0d0
do iG=1,nGrid
! Spin-up and spin-down densities
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
! Total density
r = ra + rb
! Spin-up part contribution
if(ra > threshold) then
rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
Ec(1) = Ec(1) + weight(iG)*ec_f*ra
endif
! Opposite-spin contribution
if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
ec_p = a_p*log(1d0 + b_p/rs + b_p/rs**2)
ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
z = (ra-rb)/r
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
ec_z = ec_p + (ec_f - ec_p)*fz
Ec(2) = Ec(2) + weight(iG)*ec_z*r
endif
! Spin-down contribution
if(rb > threshold) then
rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
Ec(3) = Ec(3) + weight(iG)*ec_f*rb
endif
enddo
Ec(2) = Ec(2) - Ec(1) - Ec(3)
end subroutine C16_lda_correlation_energy

131
src/eDFT/C16_lda_correlation_potential.f90
View File

@ -1,131 +0,0 @@
subroutine C16_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
! Compute unrestricted LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: mu,nu,iG
double precision :: ra,rb,r,rs
double precision :: a_p,b_p,ec_p,decdrs_p,decdra_p,decdrb_p
double precision :: a_f,b_f,ec_f,decdrs_f,decdra_f,decdrb_f
double precision :: ec_z,decdra_z,decdrb_z
double precision :: z,dzdra,dzdrb,fz,dfzdz,dfzdra,dfzdrb
double precision :: drsdra,drsdrb,dFcdra,dFcdrb
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Coefficients for Chachiyo's LDA correlation
a_p = (log(2d0) - 1d0)/(2d0*pi**2)
b_p = 20.4562557d0
a_f = (log(2d0) - 1d0)/(4d0*pi**2)
b_f = 27.4203609d0
! Compute LDA correlation matrix in the AO basis
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
! Spin-up and spin-down densities
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
! Total density
r = ra + rb
! Spin-up part contribution
if(ra > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
ec_p = a_p*log(1d0 + b_p/rs + b_p/rs**2)
ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
z = (ra-rb)/r
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
ec_z = ec_p + (ec_f - ec_p)*fz
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)
decdrs_p = - a_p/rs**2*(b_p + 2d0*b_p/rs)/(1d0 + b_p/rs + b_p/rs**2)
decdrs_f = - a_f/rs**2*(b_f + 2d0*b_f/rs)/(1d0 + b_f/rs + b_f/rs**2)
decdra_p = drsdra*decdrs_p
decdra_f = drsdra*decdrs_f
decdra_z = decdra_p + (decdra_f - decdra_p)*fz + (ec_f - ec_p)*dfzdra
dFcdra = decdra_z*r + ec_z
Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdra
endif
! Spin-down part contribution
if(rb > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
ec_p = a_p*log(1d0 + b_p/rs + b_p/rs**2)
ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
z = (ra-rb)/r
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
ec_z = ec_p + (ec_f - ec_p)*fz
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)
decdrs_p = - a_p/rs**2*(b_p + 2d0*b_p/rs)/(1d0 + b_p/rs + b_p/rs**2)
decdrs_f = - a_f/rs**2*(b_f + 2d0*b_f/rs)/(1d0 + b_f/rs + b_f/rs**2)
decdrb_p = drsdrb*decdrs_p
decdrb_f = drsdrb*decdrs_f
decdrb_z = decdrb_p + (decdrb_f - decdrb_p)*fz + (ec_f - ec_p)*dfzdrb
dFcdrb = decdrb_z*r + ec_z
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdrb
endif
enddo
enddo
enddo
end subroutine C16_lda_correlation_potential

100
src/eDFT/CC_B88_gga_exchange_energy.f90
View File

@ -1,100 +0,0 @@
subroutine CC_B88_gga_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,&
rho,drho,Cx_choice,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
integer :: iG
double precision :: b
double precision :: r,g,x
double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2
double precision :: Fx1,Fx2,Cx
! Output variables
double precision :: Ex
! Coefficients for B88 GGA exchange functional
b = 0.0042d0
! Defining enhancements factor for weight-dependent functionals
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
d1 = aCC(4,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
d2 = aCC(4,2)
w1 = wEns(2)
Fx1 = 1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4
w2 = wEns(3)
Fx2 = 1d0 + a2*w2 + b2*w2**2 + c2*w2**3 + d2*w2**4
select case (Cx_choice)
case(1)
Cx = Fx1
case(2)
Cx = Fx2
case(3)
Cx = Fx2*Fx1
case default
Cx = 1.d0
end select
! Compute GIC-GGA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
x = sqrt(g)/r**(4d0/3d0)
Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - b*x**2/(1d0 + 6d0*b*x*asinh(x)))
end if
end do
Ex = Cx*Ex
end subroutine CC_B88_gga_exchange_energy

125
src/eDFT/CC_B88_gga_exchange_potential.f90
View File

@ -1,125 +0,0 @@
subroutine CC_B88_gga_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,&
AO,dAO,rho,drho,Cx_choice,doNcentered,Fx)
! Compute the unrestricted version of the curvature-corrected exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
integer :: mu,nu,iG
double precision :: b
double precision :: vAO,gAO
double precision :: r,g,x,dxdr,dxdg,f
double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2
double precision :: Fx1,Fx2,Cx
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Coefficients for B88 GGA exchange functional
b = 0.0042d0
! Defining enhancements factor for weight-dependent functionals
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
d1 = aCC(4,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
d2 = aCC(4,2)
w1 = wEns(2)
Fx1 = 1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4
w2 = wEns(3)
Fx2 = 1d0 + a2*w2 + b2*w2**2 + c2*w2**3 + d2*w2**4
select case (Cx_choice)
case(1)
Cx = Fx1
case(2)
Cx = Fx2
case(3)
Cx = Fx2*Fx1
case default
Cx = 1.d0
end select
! Compute GGA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
x = sqrt(g)/r**(4d0/3d0)
dxdr = - 4d0*sqrt(g)/(3d0*r**(7d0/3d0))/x
dxdg = + 1d0/(2d0*sqrt(g)*r**(4d0/3d0))/x
f = b*x**2/(1d0 + 6d0*b*x*asinh(x))
Fx(mu,nu) = Fx(mu,nu) + vAO*( &
4d0/3d0*r**(1d0/3d0)*(CxLSDA - f) &
- 2d0*r**(4d0/3d0)*dxdr*f &
+ r**(4d0/3d0)*dxdr*(6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) )
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*r**(4d0/3d0)*dxdg*( &
- 2d0*f + (6d0*b*x*asinh(x) + 6d0*b*x**2/sqrt(1d0+x**2))*f/(1d0 + 6d0*b*x*asinh(x)) )
end if
end do
end do
end do
Fx(:,:) = Cx*Fx(:,:)
end subroutine CC_B88_gga_exchange_potential

170
src/eDFT/CC_lda_exchange_derivative_discontinuity.f90
View File

@ -1,170 +0,0 @@
subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,&
doNcentered,kappa,ExDD)
! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
! Local variables
integer :: iEns,jEns
integer :: iG
double precision :: r
double precision,allocatable :: dExdw(:)
double precision,external :: Kronecker_delta
double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2
double precision :: dCxdw1,dCxdw2
! Output variables
double precision,intent(out) :: ExDD(nEns)
! External variable
double precision,external :: electron_number
! Memory allocation
allocate(dExdw(nEns))
! Defining enhancements factor for weight-dependent functionals
if (doNcentered) then
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
d1 = aCC(4,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
d2 = aCC(4,2)
w1 = wEns(2)
w2 = wEns(3)
select case (Cx_choice)
case(1)
dCxdw1 = a1 + 2.d0*b1*w1 + 3.d0*c1*w1**2 + 4.d0*d1*w1**3
dCxdw2 = 0.d0
case(2)
dCxdw1 = 0.d0
dCxdw2 = a2 + 2.d0*b2*w2 + 3.d0*c2*w2**2 + 4.d0*d2*w2**3
case(3)
dCxdw1 = (a1 + 2.d0*b1*w1 + 3.d0*c1*w1**2 + 4.d0*d1*w1**3) &
* (1d0 + a2*w2 + b2*w2**2 + c2*w2**3 + d2*w2**4)
dCxdw2 = (1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4) &
* (a2 + 2.d0*b2*w2 + 3.d0*c2*w2**2 + 4.d0*d2*w2**3)
case default
dCxdw1 = 0d0
dCxdw2 = 0d0
end select
else
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
w1 = wEns(2)
w2 = wEns(3)
select case (Cx_choice)
case(1)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0)))
dCxdw2 = 0.d0
case(2)
dCxdw1 = 0.d0
dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
case(3)
dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) &
* (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2))
dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) &
* (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0)))
case default
dCxdw1 = 0d0
dCxdw2 = 0d0
end select
end if
dCxdw1 = CxLSDA*dCxdw1
dCxdw2 = CxLSDA*dCxdw2
dExdw(:) = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
if(r > threshold) then
dExdw(1) = 0d0
dExdw(2) = dExdw(2) + weight(iG)*dCxdw1*r**(4d0/3d0)
dExdw(3) = dExdw(3) + weight(iG)*dCxdw2*r**(4d0/3d0)
end if
end do
ExDD(:) = 0d0
do iEns=1,nEns
do jEns=2,nEns
if(doNcentered) then
ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - kappa(iEns)*wEns(jEns))*dExdw(jEns)
else
ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns)
end if
end do
end do
end subroutine CC_lda_exchange_derivative_discontinuity

110
src/eDFT/CC_lda_exchange_energy.f90
View File

@ -1,110 +0,0 @@
subroutine CC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
integer :: iG
double precision :: r
double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2
double precision :: Fx1,Fx2,Cx
! Output variables
double precision :: Ex
! Defining enhancements factor for weight-dependent functionals
if(doNcentered) then
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
d1 = aCC(4,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
d2 = aCC(4,2)
w1 = wEns(2)
Fx1 = 1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4
w2 = wEns(3)
Fx2 = 1d0 + a2*w2 + b2*w2**2 + c2*w2**3 + d2*w2**4
else
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
endif
select case (Cx_choice)
case(1)
Cx = CxLSDA*Fx1
case(2)
Cx = CxLSDA*Fx2
case(3)
Cx = CxLSDA*Fx2*Fx1
case default
Cx = CxLSDA
end select
! Compute GIC-LDA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) Ex = Ex + weight(iG)*Cx*r**(1d0/3d0)*r
enddo
end subroutine CC_lda_exchange_energy

131
src/eDFT/CC_lda_exchange_individual_energy.f90
View File

@ -1,131 +0,0 @@
subroutine CC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,LZx,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
integer :: iG,iEns,ispin
double precision :: r,rI
double precision :: e,dedr
double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2
double precision :: Fx1,Fx2,Cx
! Output variables
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
! Defining enhancements factor for weight-dependent functionals
if(doNcentered) then
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
d1 = aCC(4,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
d2 = aCC(4,2)
w1 = wEns(2)
Fx1 = 1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4
w2 = wEns(3)
Fx2 = 1d0 + a2*w2 + b2*w2**2 + c2*w2**3 + d2*w2**4
else
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
endif
select case (Cx_choice)
case(1)
Cx = CxLSDA*Fx1
case(2)
Cx = CxLSDA*Fx2
case(3)
Cx = CxLSDA*Fx2*Fx1
case default
Cx = CxLSDA
end select
! Compute LDA exchange matrix in the AO basis
Ex(:,:) = 0d0
LZx(:) = 0d0
do ispin=1,nspin
do iG=1,nGrid
r = max(0d0,rhow(iG,ispin))
if(r > threshold) then
e = Cx*r**(+1d0/3d0)
dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
LZx(ispin) = LZx(ispin) - weight(iG)*dedr*r*r
do iEns=1,nEns
rI = max(0d0,rho(iG,ispin,iEns))
if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(e+dedr*r)*rI
end do
endif
enddo
enddo
end subroutine CC_lda_exchange_individual_energy

119
src/eDFT/CC_lda_exchange_potential.f90
View File

@ -1,119 +0,0 @@
subroutine CC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx)
! Compute the unrestricted version of the curvature-corrected exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
integer :: mu,nu,iG
double precision :: r,vAO
double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2
double precision :: Fx1,Fx2,Cx
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Defining enhancements factor for weight-dependent functionals
if(doNcentered) then
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
d1 = aCC(4,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
d2 = aCC(4,2)
w1 = wEns(2)
Fx1 = 1d0 + a1*w1 + b1*w1**2 + c1*w1**3 + d1*w1**4
w2 = wEns(3)
Fx2 = 1d0 + a2*w2 + b2*w2**2 + c2*w2**3 + d2*w2**4
else
! Parameters for first state
a1 = aCC(1,1)
b1 = aCC(2,1)
c1 = aCC(3,1)
! Parameters for second state
a2 = aCC(1,2)
b2 = aCC(2,2)
c2 = aCC(3,2)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)
w2 = wEns(3)
Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)
endif
select case (Cx_choice)
case(1)
Cx = CxLSDA*Fx1
case(2)
Cx = CxLSDA*Fx2
case(3)
Cx = CxLSDA*Fx2*Fx1
case default
Cx = CxLSDA
end select
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*Cx*r**(1d0/3d0)
endif
enddo
enddo
enddo
end subroutine CC_lda_exchange_potential

48
src/eDFT/G96_gga_exchange_energy.f90
View File

@ -1,48 +0,0 @@
subroutine G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute Gill's 96 GGA exchange energy
implicit none
include 'parameters.h'
! Input variables
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
integer :: iG
double precision :: beta
double precision :: r,g
! Output variables
double precision :: Ex
! Coefficients for G96 GGA exchange functional
beta = 1d0/137d0
! Compute GGA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
Ex = Ex + weight(iG)*r**(4d0/3d0)*(CxLSDA - beta*g**(3d0/4d0)/r**2)
end if
end do
end subroutine G96_gga_exchange_energy

64
src/eDFT/G96_gga_exchange_potential.f90
View File

@ -1,64 +0,0 @@
subroutine G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute Gill's GGA exchange poential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: beta
double precision :: r,g,vAO,gAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Coefficients for G96 GGA exchange functional
beta = +1d0/137d0
! Compute GGA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(CxLSDA - beta*g**(3d0/4d0)/r**2) &
+ 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0))
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) - 2d0*gAO*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0)
endif
enddo
enddo
enddo
end subroutine G96_gga_exchange_potential

73
src/eDFT/LYP_gga_correlation_energy.f90
View File

@ -1,73 +0,0 @@
subroutine LYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
! Compute unrestricted LYP GGA correlation energy
implicit none
include 'parameters.h'
! Input variables
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
integer :: iG
double precision :: ra,rb,r
double precision :: ga,gab,gb,g
double precision :: a,b,c,d
double precision :: Cf,omega,delta
! Output variables
double precision :: Ec(nsp)
! Parameters of the functional
a = 0.04918d0
b = 0.132d0
c = 0.2533d0
d = 0.349d0
Cf = 3d0/10d0*(3d0*pi**2)**(2d0/3d0)
! Initialization
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
if(r > threshold) then
ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1)
gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2)
gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2)
g = ga + 2d0*gab + gb
omega = exp(-c*r**(-1d0/3d0))/(1d0 + d*r**(-1d0/3d0))*r**(-11d0/3d0)
delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0))
Ec(2) = Ec(2) - weight(iG)*4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r &
- weight(iG)*a*b*omega*ra*rb*( &
2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g &
- (5d0/2d0 - delta/18d0)*(ga + gb) &
- (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) ) &
- weight(iG)*a*b*omega*( &
- 2d0*r**2/3d0*g &
+ (2d0*r**2/3d0 - ra**2)*gb &
+ (2d0*r**2/3d0 - rb**2)*ga )
end if
end do
end subroutine LYP_gga_correlation_energy

156
src/eDFT/LYP_gga_correlation_potential.f90
View File

@ -1,156 +0,0 @@
subroutine LYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute LYP correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(ncart,nGrid,nspin)
! Local variables
integer :: mu,nu,iG
double precision :: vAO,gaAO,gbAO
double precision :: ra,rb,r
double precision :: ga,gab,gb,g
double precision :: dfdra,dfdrb
double precision :: dfdga,dfdgab,dfdgb
double precision :: dodra,dodrb,dddra,dddrb
double precision :: a,b,c,d
double precision :: Cf,omega,delta
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Prameter of the functional
a = 0.04918d0
b = 0.132d0
c = 0.2533d0
d = 0.349d0
Cf = 3d0/10d0*(3d0*pi**2)**(2d0/3d0)
! Compute matrix elements in the AO basis
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
if(r > threshold) then
ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1)
gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2)
gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2)
g = ga + 2d0*gab + gb
omega = exp(-c*r**(-1d0/3d0))/(1d0 + d*r**(-1d0/3d0))*r**(-11d0/3d0)
delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0))
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
dodra = (d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) + c/(3d0*r**(4d0/3d0)) - 11d0/(3d0*r))*omega
dodrb = dodra
dddra = - c/3d0*r**(-4d0/3d0) &
+ d**2/(3d0*(1d0 + d*r**(-1d0/3d0))**2)*r**(-5d0/3d0) &
- d/(3d0*(1d0 + d*r**(-1d0/3d0)))*r**(-4d0/3d0)
dddrb = dddra
dfdra = - 4d0*a/(1d0 + d*r**(-1d0/3d0))*rb/r &
- 4d0/3d0*a*d/(1d0 + d*r**(-1d0/3d0))**2*ra*rb/r**(7d0/3d0) &
+ 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r**2 &
- a*b*omega*rb*( &
+ 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g &
- (5d0/2d0 - delta/18d0)*(ga + gb) &
- (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) &
- 4d0/3d0*r/rb*g &
+ (4d0/3d0*r/rb - 2d0*ra/rb)*gb &
+ 4d0/3d0*r/rb*ga ) &
- a*b*omega*ra*rb*( &
+ 8d0/3d0*2d0**(11d0/3d0)*Cf*ra**(5d0/3d0) &
- 7d0*dddra/18d0*g &
+ dddra/18d0*(ga + gb) &
- dddra/9d0*(ra/r*ga + rb/r*gb) &
- (delta - 11d0)/(9d0*r)*(-ra/r*ga - rb/r*gb + ga) ) &
- a*b*dodra*ra*rb*( &
+ 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g &
- (5d0/2d0 - delta/18d0)*(ga + gb) &
- (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) ) &
- a*b*dodra*( &
- 2d0*r**2/3d0*g &
+ (2d0*r**2/3d0 - ra**2)*gb &
+ (2d0*r**2/3d0 - rb**2)*ga )
dfdrb = - 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra/r &
- 4d0/3d0*a*d/(1d0 + d*r**(-1d0/3d0))**2*ra*rb/r**(7d0/3d0) &
+ 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r**2 &
- a*b*omega*ra*( &
+ 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g &
- (5d0/2d0 - delta/18d0)*(ga + gb) &
- (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) &
- 4d0/3d0*r/ra*g &
+ (4d0/3d0*r/ra - 2d0*rb/ra)*ga &
+ 4d0/3d0*r/ra*gb ) &
- a*b*omega*ra*rb*( &
+ 8d0/3d0*2d0**(11d0/3d0)*Cf*rb**(5d0/3d0) &
- 7d0*dddrb/18d0*g &
+ dddrb/18d0*(ga + gb) &
- dddrb/9d0*(ra/r*ga + rb/r*gb) &
- (delta - 11d0)/(9d0*r)*(-ra/r*ga - rb/r*gb + gb) ) &
- a*b*dodrb*ra*rb*( &
+ 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g &
- (5d0/2d0 - delta/18d0)*(ga + gb) &
- (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) ) &
- a*b*dodrb*( &
- 2d0*r**2/3d0*g &
+ (2d0*r**2/3d0 - ra**2)*gb &
+ (2d0*r**2/3d0 - rb**2)*ga )
Fc(mu,nu,1) = Fc(mu,nu,1) + vAO*dfdra
Fc(mu,nu,2) = Fc(mu,nu,2) + vAO*dfdrb
gaAO = drho(1,iG,1)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG,1)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG,1)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gaAO = weight(iG)*gaAO
gbAO = drho(1,iG,2)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG,2)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG,2)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gbAO = weight(iG)*gbAO
dfdga = -a*b*omega*(-rb**2 + ra*rb*(1d0/9d0 - (delta-11d0)/9d0*ra/r - delta/3d0))
dfdgab = -a*b*omega*(-4d0/3d0*r**2 + 2d0*ra*rb*(47d0/18d0 - 7d0*delta/18d0))
dfdgb = -a*b*omega*(-ra**2 + ra*rb*(1d0/9d0 - (delta-11d0)/9d0*rb/r - delta/3d0))
Fc(mu,nu,1) = Fc(mu,nu,1) + 2d0*gaAO*dfdga + gbAO*dfdgab
Fc(mu,nu,2) = Fc(mu,nu,2) + 2d0*gbAO*dfdgb + gaAO*dfdgab
end if
end do
end do
end do
end subroutine LYP_gga_correlation_potential

172
src/eDFT/PBE_gga_correlation_energy.f90
View File

@ -1,172 +0,0 @@
subroutine PBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
! Compute unrestricted PBE GGA correlation energy
implicit none
include 'parameters.h'
! Input variables
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
integer :: iG
double precision :: ra,rb,r,rs,z
double precision :: ga,gab,gb,g
double precision :: a,b,c,d
double precision :: gam,beta
double precision :: A_p,a1_p,b1_p,b2_p,b3_p,b4_p
double precision :: A_f,a1_f,b1_f,b2_f,b3_f,b4_f
double precision :: A_a,a1_a,b1_a,b2_a,b3_a,b4_a
double precision :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
double precision :: H,kf,ks,t,phi
! Output variables
double precision :: Ec(nsp)
! Parameters for PW92
A_p = 0.031091d0
a1_p = 0.21370d0
b1_p = 7.5957d0
b2_p = 3.5876d0
b3_p = 1.6382d0
b4_p = 0.49294d0
A_f = 0.015545d0
a1_f = 0.20548d0
b1_f = 14.1189d0
b2_f = 6.1977d0
b3_f = 3.3662d0
b4_f = 0.62517d0
A_a = 0.016887d0
a1_a = 0.11125d0
b1_a = 10.357d0
b2_a = 3.6231d0
b3_a = 0.88026d0
b4_a = 0.49671d0
! Parameters PBE
gam = (1d0 - log(2d0))/pi**2
beta = 0.066725d0
! Initialization
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
z = (ra - rb)/r
! alpha-alpha contribution
if(ra > threshold) then
rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1)
kf = (3d0*pi**2*ra)**(1d0/3d0)
ks = sqrt(4d0*kf/pi)
phi = 1d0
t = sqrt(ga)/(2d0*phi*ks*ra)
A = beta/gam/(exp(-ec_f/(gam*phi**3)) - 1d0)
H = gam*phi**3*log(1d0 + beta/gam*t**2*((1d0 + A*t**2)/(1d0 + A*t**2 + A**2*t**4)))
Ec(1) = Ec(1) + weight(iG)*(ec_f + H)*ra
end if
r = ra + rb
! alpha-beta contribution
if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
ec_p = b1_p*sqrt(rs) + b2_p*rs + b3_p*rs**(3d0/2d0) + b4_p*rs**2
ec_p = -2d0*A_p*(1d0 + a1_p*rs)*log(1d0 + 1d0/(2d0*A_p*ec_p))
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
ec_a = b1_a*sqrt(rs) + b2_a*rs + b3_a*rs**(3d0/2d0) + b4_a*rs**2
ec_a = +2d0*A_a*(1d0 + a1_a*rs)*log(1d0 + 1d0/(2d0*A_a*ec_a))
ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1)
gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2)
gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2)
g = ga + 2d0*gab + gb
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
kf = (3d0*pi**2*r)**(1d0/3d0)
ks = sqrt(4d0*kf/pi)
phi = ((1d0 + z)**(2d0/3d0) + (1d0 - z)**(2d0/3d0))/2d0
t = sqrt(g)/(2d0*phi*ks*r)
A = beta/gam/(exp(-ec_p/(gam*phi**3)) - 1d0)
H = gam*phi**3*log(1d0 + beta/gam*t**2*((1d0 + A*t**2)/(1d0 + A*t**2 + A**2*t**4)))
Ec(2) = Ec(2) - weight(iG)*(ec_p + H)*r
end if
! beta-beta contribution
if(rb > threshold) then
rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2)
kf = (3d0*pi**2*rb)**(1d0/3d0)
ks = sqrt(4d0*kf/pi)
phi = 1d0
t = sqrt(gb)/(2d0*phi*ks*rb)
A = beta/gam/(exp(-ec_f/(gam*phi**3)) - 1d0)
H = gam*phi**3*log(1d0 + beta/gam*t**2*((1d0 + A*t**2)/(1d0 + A*t**2 + A**2*t**4)))
Ec(3) = Ec(3) + weight(iG)*(ec_f + H)*rb
end if
end do
Ec(2) = Ec(2) - Ec(1) - Ec(3)
end subroutine PBE_gga_correlation_energy

88
src/eDFT/PBE_gga_correlation_potential.f90
View File

@ -1,88 +0,0 @@
subroutine PBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute LYP correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(ncart,nGrid,nspin)
! Local variables
integer :: mu,nu,iG
double precision :: vAO,gaAO,gbAO
double precision :: ra,rb,r
double precision :: ga,gab,gb,g
double precision :: dfdra,dfdrb
double precision :: dfdga,dfdgab,dfdgb
double precision :: dodra,dodrb,dddra,dddrb
double precision :: a,b,c,d
double precision :: Cf,omega,delta
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Prameter of the functional
! Compute matrix elements in the AO basis
call PW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
if(r > threshold) then
ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1)
gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2)
gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2)
g = ga + 2d0*gab + gb
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
dfdra = 0d0
dfdrb = 0d0
Fc(mu,nu,1) = Fc(mu,nu,1) + vAO*dfdra
Fc(mu,nu,2) = Fc(mu,nu,2) + vAO*dfdrb
gaAO = drho(1,iG,1)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG,1)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG,1)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gaAO = weight(iG)*gaAO
gbAO = drho(1,iG,2)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG,2)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG,2)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gbAO = weight(iG)*gbAO
dfdga = 0d0
dfdgab = 0d0
dfdgb = 0d0
Fc(mu,nu,1) = Fc(mu,nu,1) + 2d0*gaAO*dfdga + gbAO*dfdgab
Fc(mu,nu,2) = Fc(mu,nu,2) + 2d0*gbAO*dfdgb + gaAO*dfdgab
end if
end do
end do
end do
end subroutine PBE_gga_correlation_potential

49
src/eDFT/PBE_gga_exchange_energy.f90
View File

@ -1,49 +0,0 @@
subroutine PBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute PBE GGA exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
integer :: iG
double precision :: mupbe,kappa
double precision :: r,g,s2
! Output variables
double precision :: Ex
! Coefficients for PBE exchange functional
mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0
kappa = 0.804d0
! Compute GGA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
s2 = g/r**(8d0/3d0)
Ex = Ex + weight(iG)*CxLSDA*r**(4d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa))
end if
end do
end subroutine PBE_gga_exchange_energy

67
src/eDFT/PBE_gga_exchange_potential.f90
View File

@ -1,67 +0,0 @@
subroutine PBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute PBE GGA exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: mupbe,kappa
double precision :: r,g,s2,vAO,gAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Coefficients for PBE exchange functional
mupbe = ((1d0/2d0**(1d0/3d0))/(2d0*(3d0*pi**2)**(1d0/3d0)))**2*0.21951d0
kappa = 0.804d0
! Compute GGA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
s2 = g/r**(8d0/3d0)
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*4d0/3d0*CxLSDA*r**(1d0/3d0)*(1d0 + kappa - kappa/(1d0 + mupbe*s2/kappa)) &
- vAO*8d0/3d0*CxLSDA*r**(1d0/3d0)*mupbe*s2/(1d0 + mupbe*s2/kappa)**2
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) + 2d0*gAO*CxLSDA*r**(-4d0/3d0)*mupbe/(1d0 + mupbe*s2/kappa)**2
end if
end do
end do
end do
end subroutine PBE_gga_exchange_potential

120
src/eDFT/PW92_lda_correlation_energy.f90
View File

@ -1,120 +0,0 @@
subroutine PW92_lda_correlation_energy(nGrid,weight,rho,Ec)
! Compute unrestricted PW92 LDA correlation energy
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,z
double precision :: A_p,a1_p,b1_p,b2_p,b3_p,b4_p
double precision :: A_f,a1_f,b1_f,b2_f,b3_f,b4_f
double precision :: A_a,a1_a,b1_a,b2_a,b3_a,b4_a
double precision :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
! Output variables
double precision :: Ec(nsp)
! Parameters of the functional
A_p = 0.031091d0
a1_p = 0.21370d0
b1_p = 7.5957d0
b2_p = 3.5876d0
b3_p = 1.6382d0
b4_p = 0.49294d0
A_f = 0.015545d0
a1_f = 0.20548d0
b1_f = 14.1189d0
b2_f = 6.1977d0
b3_f = 3.3662d0
b4_f = 0.62517d0
A_a = 0.016887d0
a1_a = 0.11125d0
b1_a = 10.357d0
b2_a = 3.6231d0
b3_a = 0.88026d0
b4_a = 0.49671d0
! Initialization
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
z = (ra - rb)/r
! alpha-alpha contribution
if(ra > threshold) then
rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
Ec(1) = Ec(1) + weight(iG)*ec_f*ra
end if
! alpha-beta contribution
if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
ec_p = b1_p*sqrt(rs) + b2_p*rs + b3_p*rs**(3d0/2d0) + b4_p*rs**2
ec_p = -2d0*A_p*(1d0 + a1_p*rs)*log(1d0 + 1d0/(2d0*A_p*ec_p))
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
ec_a = b1_a*sqrt(rs) + b2_a*rs + b3_a*rs**(3d0/2d0) + b4_a*rs**2
ec_a = +2d0*A_a*(1d0 + a1_a*rs)*log(1d0 + 1d0/(2d0*A_a*ec_a))
ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
Ec(2) = Ec(2) + weight(iG)*ec_z*r
end if
! beta-beta contribution
if(rb > threshold) then
rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
Ec(3) = Ec(3) + weight(iG)*ec_f*rb
end if
end do
Ec(2) = Ec(2) - Ec(1) - Ec(3)
end subroutine PW92_lda_correlation_energy

185
src/eDFT/PW92_lda_correlation_potential.f90
View File

@ -1,185 +0,0 @@
subroutine PW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
! Compute unrestricted PW92 LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: mu,nu,iG
double precision :: ra,rb,r,rs,z,t,dt
double precision :: A_p,a1_p,b1_p,b2_p,b3_p,b4_p
double precision :: A_f,a1_f,b1_f,b2_f,b3_f,b4_f
double precision :: A_a,a1_a,b1_a,b2_a,b3_a,b4_a
double precision :: dfzdz,decdrs_p,decdrs_f,decdrs_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 :: Fc(nBas,nBas,nspin)
! Parameters of the functional
A_p = 0.031091d0
a1_p = 0.21370d0
b1_p = 7.5957d0
b2_p = 3.5876d0
b3_p = 1.6382d0
b4_p = 0.49294d0
A_f = 0.015545d0
a1_f = 0.20548d0
b1_f = 14.1189d0
b2_f = 6.1977d0
b3_f = 3.3662d0
b4_f = 0.62517d0
A_a = 0.016887d0
a1_a = 0.11125d0
b1_a = 10.357d0
b2_a = 3.6231d0
b3_a = 0.88026d0
b4_a = 0.49671d0
! Initialization
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
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
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
ec_p = b1_p*sqrt(rs) + b2_p*rs + b3_p*rs**(3d0/2d0) + b4_p*rs**2
ec_p = -2d0*A_p*(1d0 + a1_p*rs)*log(1d0 + 1d0/(2d0*A_p*ec_p))
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
ec_a = b1_a*sqrt(rs) + b2_a*rs + b3_a*rs**(3d0/2d0) + b4_a*rs**2
ec_a = +2d0*A_a*(1d0 + a1_a*rs)*log(1d0 + 1d0/(2d0*A_a*ec_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)
t = b1_p*sqrt(rs) + b2_p*rs + b3_p*rs**(3d0/2d0) + b4_p*rs**2
dt = 0.5d0*b1_p*sqrt(rs) + b2_p*rs + 1.5d0*b3_p*rs**(3d0/2d0) + 2d0*b4_p*rs**2
decdrs_p = (1d0 + a1_p*rs)/rs*dt/(t**2*(1d0 + 1d0/(2d0*A_p*t))) &
- 2d0*A_p*a1_p*log(1d0 + 1d0/(2d0*A_p*t))
t = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
dt = 0.5d0*b1_f*sqrt(rs) + b2_f*rs + 1.5d0*b3_f*rs**(3d0/2d0) + 2d0*b4_f*rs**2
decdrs_f = (1d0 + a1_f*rs)/rs*dt/(t**2*(1d0 + 1d0/(2d0*A_f*t))) &
- 2d0*A_f*a1_f*log(1d0 + 1d0/(2d0*A_f*t))
t = b1_a*sqrt(rs) + b2_a*rs + b3_a*rs**(3d0/2d0) + b4_a*rs**2
dt = 0.5d0*b1_a*sqrt(rs) + b2_a*rs + 1.5d0*b3_a*rs**(3d0/2d0) + 2d0*b4_a*rs**2
decdrs_a = (1d0 + a1_a*rs)/rs*dt/(t**2*(1d0 + 1d0/(2d0*A_a*t))) &
- 2d0*A_a*a1_a*log(1d0 + 1d0/(2d0*A_a*t))
decdra_p = drsdra*decdrs_p
decdra_f = drsdra*decdrs_f
decdra_a = drsdra*decdrs_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
Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdra*r)
end if
! spin-down contribution
if(rb > threshold) then
r = ra + rb
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
ec_p = b1_p*sqrt(rs) + b2_p*rs + b3_p*rs**(3d0/2d0) + b4_p*rs**2
ec_p = -2d0*A_p*(1d0 + a1_p*rs)*log(1d0 + 1d0/(2d0*A_p*ec_p))
ec_f = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
ec_f = -2d0*A_f*(1d0 + a1_f*rs)*log(1d0 + 1d0/(2d0*A_f*ec_f))
ec_a = b1_a*sqrt(rs) + b2_a*rs + b3_a*rs**(3d0/2d0) + b4_a*rs**2
ec_a = +2d0*A_a*(1d0 + a1_a*rs)*log(1d0 + 1d0/(2d0*A_a*ec_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)
t = b1_p*sqrt(rs) + b2_p*rs + b3_p*rs**(3d0/2d0) + b4_p*rs**2
dt = 0.5d0*b1_p*sqrt(rs) + b2_p*rs + 1.5d0*b3_p*rs**(3d0/2d0) + 2d0*b4_p*rs**2
decdrs_p = (1d0 + a1_p*rs)/rs*dt/(t**2*(1d0 + 1d0/(2d0*A_p*t))) &
- 2d0*A_p*a1_p*log(1d0 + 1d0/(2d0*A_p*t))
t = b1_f*sqrt(rs) + b2_f*rs + b3_f*rs**(3d0/2d0) + b4_f*rs**2
dt = 0.5d0*b1_f*sqrt(rs) + b2_f*rs + 1.5d0*b3_f*rs**(3d0/2d0) + 2d0*b4_f*rs**2
decdrs_f = (1d0 + a1_f*rs)/rs*dt/(t**2*(1d0 + 1d0/(2d0*A_f*t))) &
- 2d0*A_f*a1_f*log(1d0 + 1d0/(2d0*A_f*t))
t = b1_a*sqrt(rs) + b2_a*rs + b3_a*rs**(3d0/2d0) + b4_a*rs**2
dt = 0.5d0*b1_a*sqrt(rs) + b2_a*rs + 1.5d0*b3_a*rs**(3d0/2d0) + 2d0*b4_a*rs**2
decdrs_a = (1d0 + a1_a*rs)/rs*dt/(t**2*(1d0 + 1d0/(2d0*A_a*t))) &
- 2d0*A_a*a1_a*log(1d0 + 1d0/(2d0*A_a*t))
decdrb_p = drsdrb*decdrs_p
decdrb_f = drsdrb*decdrs_f
decdrb_a = drsdrb*decdrs_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
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
end if
end do
end do
end do
end subroutine PW92_lda_correlation_potential

34
src/eDFT/S51_lda_exchange_energy.f90
View File

@ -1,34 +0,0 @@
subroutine S51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Compute Slater's LDA exchange energy
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
! Output variables
double precision :: Ex
! Compute LDA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) Ex = Ex + weight(iG)*CxLSDA*r**(1d0/3d0)*r
enddo
end subroutine S51_lda_exchange_energy

61
src/eDFT/S51_lda_exchange_individual_energy.f90
View File

@ -1,61 +0,0 @@
subroutine S51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex)
! Compute the restricted version of Slater's LDA exchange individual energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
! Local variables
integer :: iG
integer :: iEns
integer :: ispin
double precision :: r
double precision :: rI
double precision :: e
double precision :: dedr
! Output variables
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
LZx(:) = 0d0
Ex(:,:) = 0d0
do ispin=1,nspin
do iG=1,nGrid
r = max(0d0,rhow(iG,ispin))
if(r > threshold) then
e = CxLSDA*r**(+1d0/3d0)
dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0)
LZx(ispin) = LZx(ispin) - weight(iG)*dedr*r*r
do iEns=1,nEns
rI = max(0d0,rho(iG,ispin,iEns))
if(rI > threshold) Ex(ispin,iEns) = Ex(ispin,iEns) + weight(iG)*(e + dedr*r)*rI
end do
endif
enddo
enddo
end subroutine S51_lda_exchange_individual_energy

45
src/eDFT/S51_lda_exchange_potential.f90
View File

@ -1,45 +0,0 @@
subroutine S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
! Compute Slater's LDA exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: r,vAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLSDA*r**(1d0/3d0)
endif
enddo
enddo
enddo
end subroutine S51_lda_exchange_potential

392
src/eDFT/UKS.f90
View File

@ -1,392 +0,0 @@
subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,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,c,Pw,Vxc)
! Perform unrestricted Kohn-Sham calculation for ensembles
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: x_rung,c_rung
integer,intent(in) :: x_DFA,c_DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
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)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
integer :: xc_rung
logical :: LDA_centered = .false.
integer :: nSCF
integer :: nBasSq
integer :: n_diis
integer :: nO(nspin,nEns)
double precision :: Conv
double precision :: rcond(nspin)
double precision :: ET(nspin)
double precision :: EV(nspin)
double precision :: EH(nsp)
double precision :: Ex(nspin)
double precision :: Ec(nsp)
double precision :: dipole(ncart)
double precision,allocatable :: cp(:,:,:)
double precision,allocatable :: J(:,:,:)
double precision,allocatable :: F(:,:,:)
double precision,allocatable :: Fp(:,:,:)
double precision,allocatable :: Fx(:,:,:)
double precision,allocatable :: FxHF(:,:,:)
double precision,allocatable :: Fc(:,:,:)
double precision,allocatable :: err(:,:,:)
double precision,allocatable :: err_diis(:,:,:)
double precision,allocatable :: F_diis(:,:,:)
double precision,external :: trace_matrix
double precision,external :: electron_number
double precision,allocatable :: rhow(:,:)
double precision,allocatable :: drhow(:,:,:)
double precision :: nEl(nspin)
double precision,allocatable :: P(:,:,:,:)
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: drho(:,:,:,:)
integer :: ispin,iEns,iBas
! Output variables
double precision,intent(out) :: Ew
double precision,intent(out) :: eKS(nBas,nspin)
double precision,intent(out) :: Pw(nBas,nBas,nspin)
double precision,intent(out) :: c(nBas,nBas,nspin)
double precision,intent(out) :: Vxc(nBas,nspin)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'* Unrestricted Kohn-Sham calculation *'
write(*,*)'* *** for ensembles *** *'
write(*,*)'************************************************'
write(*,*)
! Useful stuff
nBasSq = nBas*nBas
!------------------------------------------------------------------------
! Rung of Jacob's ladder
!------------------------------------------------------------------------
xc_rung = max(x_rung,c_rung)
! Memory allocation
allocate(cp(nBas,nBas,nspin),J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), &
Fx(nBas,nBas,nspin),FxHF(nBas,nBas,nspin),Fc(nBas,nBas,nspin),err(nBas,nBas,nspin), &
rhow(nGrid,nspin),drhow(ncart,nGrid,nspin), &
err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin), &
P(nBas,nBas,nspin,nEns),rho(nGrid,nspin,nEns),drho(ncart,nGrid,nspin,nEns))
! Guess coefficients and eigenvalues
nO(:,:) = 0
do iEns=1,nEns
do ispin=1,nspin
nO(ispin,iEns) = int(sum(occnum(:,ispin,iEns)))
end do
end do
do ispin=1,nspin
call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
end do
! Mix guess for UHF solution in singlet states
if(mix) then
write(*,*) '!! guess mixing disabled in UKS !!'
write(*,*)
end if
! Initialization
nSCF = 0
conv = 1d0
nEl(:) = 0d0
Ex(:) = 0d0
Ec(:) = 0d0
Fx(:,:,:) = 0d0
FxHF(:,:,:) = 0d0
Fc(:,:,:) = 0d0
n_diis = 0
F_diis(:,:,:) = 0d0
err_diis(:,:,:) = 0d0
rcond(:) = 1d0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'------------------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|'
write(*,*)'------------------------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
!------------------------------------------------------------------------
! Compute density matrix
!------------------------------------------------------------------------
call density_matrix(nBas,nEns,c(:,:,:),P(:,:,:,:),occnum(:,:,:))
! Weight-dependent density matrix
Pw(:,:,:) = 0d0
do iEns=1,nEns
Pw(:,:,:) = Pw(:,:,:) + wEns(iEns)*P(:,:,:,iEns)
end do
!------------------------------------------------------------------------
! Compute one-electron density and its gradient if necessary
!------------------------------------------------------------------------
do ispin=1,nspin
do iEns=1,nEns
call density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),rho(:,ispin,iEns))
end do
end do
! Weight-dependent one-electron density
rhow(:,:) = 0d0
do iEns=1,nEns
rhow(:,:) = rhow(:,:) + wEns(iEns)*rho(:,:,iEns)
end do
if(xc_rung > 1) then
! Ground state density
do ispin=1,nspin
do iEns=1,nEns
call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns))
end do
end do
! Weight-dependent one-electron density
drhow(:,:,:) = 0d0
do iEns=1,nEns
drhow(:,:,:) = drhow(:,:,:) + wEns(iEns)*drho(:,:,:,iEns)
end do
end if
!------------------------------------------------------------------------
! Compute Hxc potential and Fock operator
!------------------------------------------------------------------------
! Compute Hartree potential
do ispin=1,nspin
call hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin))
end do
! Compute exchange potential
do ispin=1,nspin
call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
Pw(:,:,ispin),ERI,AO,dAO,rhow(:,ispin),drhow(:,:,ispin), &
Cx_choice,doNcentered,Fx(:,:,ispin),FxHF(:,:,ispin))
end do
! Compute correlation potential
call correlation_potential(c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rhow,drhow,Fc)
! Build Fock operator
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + Fx(:,:,ispin) + Fc(:,:,ispin)
end do
! Check convergence
do ispin=1,nspin
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin))
end do
if(nSCF > 1) Conv = maxval(abs(err(:,:,:)))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
do ispin=1,nspin
if(rcond(ispin) > 1d-15) then
call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, &
err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
else
n_diis = 0
end if
end do
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(level_shift,nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin))
end do
end if
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
!------------------------------------------------------------------------
! Compute KS energy
!------------------------------------------------------------------------
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),T(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),V(:,:)))
end do
! Hartree energy
call hartree_energy(nBas,Pw,J,EH)
! Exchange energy
do ispin=1,nspin
call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin), &
Cx_choice,doNcentered,Ex(ispin))
end do
! Correlation energy
call correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
! Total energy
Ew = sum(ET(:)) + sum(EV(:)) + sum(EH(:)) + sum(Ex(:)) + sum(Ec(:))
! Check the grid accuracy by computing the number of electrons
do ispin=1,nspin
nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin))
end do
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',Conv,'|',sum(nEl(:)),'|'
end do
write(*,*)'------------------------------------------------------------------------------------------'
! print*,'Ensemble energy:',Ew + ENuc,'au'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
end if
! Compute final KS energy
call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_UKS(nBas,nEns,occnum,S,wEns,eKS,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
! Compute Vxc for post-HF calculations
call xc_potential(nBas,c,Fx,Fc,Vxc)
!------------------------------------------------------------------------
! Compute individual energies from ensemble energy
!------------------------------------------------------------------------
call individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew)
end subroutine UKS

137
src/eDFT/VWN3_lda_correlation_energy.f90
View File

@ -1,137 +0,0 @@
subroutine VWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
! Compute unrestricted VWN3 LDA correlation energy
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 :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
! Output variables
double precision :: Ec(nsp)
! Parameters of the functional
a_p = +0.0621814d0/2d0
x0_p = -0.409286d0
b_p = +13.0720d0
c_p = +42.7198d0
a_f = +0.0621814d0/4d0
x0_f = -0.743294d0
b_f = +20.1231d0
c_f = +101.578d0
a_a = -1d0/(6d0*pi**2)
x0_a = -0.0047584D0
b_a = 1.13107d0
c_a = 13.0045d0
! Initialization
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
z = (ra - rb)/r
! alpha-alpha contribution
if(ra > threshold) then
rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
x = sqrt(rs)
x_f = x*x + b_f*x + c_f
xx0_f = x0_f*x0_f + b_f*x0_f + c_f
q_f = sqrt(4d0*c_f - b_f*b_f)
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(1) = Ec(1) + weight(iG)*ec_f*ra
end if
! alpha-beta contribution
if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
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
Ec(2) = Ec(2) + weight(iG)*ec_z*r
end if
! beta-beta contribution
if(rb > threshold) then
rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
x = sqrt(rs)
x_f = x*x + b_f*x + c_f
xx0_f = x0_f*x0_f + b_f*x0_f + c_f
q_f = sqrt(4d0*c_f - b_f*b_f)
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(3) = Ec(3) + weight(iG)*ec_f*rb
end if
end do
Ec(2) = Ec(2) - Ec(1) - Ec(3)
end subroutine VWN3_lda_correlation_energy

181
src/eDFT/VWN3_lda_correlation_individual_energy.f90
View File

@ -1,181 +0,0 @@
subroutine VWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
! Compute VWN3 LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
logical,intent(in) :: doNcentered
! Local variables
integer :: iG
integer :: iEns
double precision :: ra,rb,r,raI,rbI,rI,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,dzdrb,dfzdra,dfzdrb,drsdr,decdr_p,decdr_f,decdr_a,decdra,decdrb
double precision :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
! Output variables
double precision,intent(out) :: LZc(nsp)
double precision,intent(out) :: Ec(nsp,nEns)
! Parameters of the functional
a_p = +0.0621814d0/2d0
x0_p = -0.409286d0
b_p = +13.0720d0
c_p = +42.7198d0
a_f = +0.0621814d0/4d0
x0_f = -0.743294d0
b_f = +20.1231d0
c_f = +101.578d0
a_a = -1d0/(6d0*pi**2)
x0_a = -0.0047584D0
b_a = +1.13107d0
c_a = +13.0045d0
! Initialization
LZc(:) = 0d0
Ec(:,:) = 0d0
do iG=1,nGrid
ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2))
r = ra + rb
! up-down contribution
if(r > threshold) then
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
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
drsdr = - (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*( 2d0/(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*( 2d0/(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*( 2d0/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
decdr_p = drsdr*dxdrs*decdx_p
decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a
dzdra = + (1d0 - z)/r
dfzdra = dzdra*dfzdz
decdra = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+ (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
dzdrb = - (1d0 + z)/r
dfzdrb = dzdrb*dfzdz
decdrb = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+ (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
! spin-up contribution
if(ra > threshold) then
LZc(1) = LZc(1) - weight(iG)*decdra*ra*ra
do iEns=1,nEns
raI = max(0d0,rho(iG,1,iEns))
if(raI > threshold) then
Ec(1,iEns) = Ec(1,iEns) + weight(iG)*(ec_z + decdra*ra)*raI
Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdra*rb*raI
end if
end do
end if
! spin-down contribution
if(rb > threshold) then
LZc(3) = LZc(3) - weight(iG)*decdrb*rb*rb
do iEns=1,nEns
rbI = max(0d0,rho(iG,2,iEns))
if(rbI > threshold) then
Ec(3,iEns) = Ec(3,iEns) + weight(iG)*(ec_z + decdrb*rb)*rbI
Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdrb*ra*rbI
end if
end do
end if
end if
end do
end subroutine VWN3_lda_correlation_individual_energy

196
src/eDFT/VWN3_lda_correlation_potential.f90
View File

@ -1,196 +0,0 @@
subroutine VWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
! Compute unrestricted VWN3 LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: mu,nu,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 :: Fc(nBas,nBas,nspin)
! Parameters of the functional
a_p = +0.0621814d0/2d0
x0_p = -0.409286d0
b_p = +13.0720d0
c_p = +42.7198d0
a_f = +0.0621814d0/4d0
x0_f = -0.743294d0
b_f = +20.1231d0
c_f = +101.578d0
a_a = -1d0/(6d0*pi**2)
x0_a = -0.0047584D0
b_a = +1.13107d0
c_a = +13.0045d0
! Initialization
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
z = (ra - rb)/r
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
! spin-up contribution
if(ra > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs)
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
Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdra*r)
end if
! spin-down contribution
if(rb > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs)
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
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
end if
end do
end do
end do
end subroutine VWN3_lda_correlation_potential

137
src/eDFT/VWN5_lda_correlation_energy.f90
View File

@ -1,137 +0,0 @@
subroutine VWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
! Compute unrestricted VWN5 LDA correlation energy
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 :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
! Output variables
double precision :: Ec(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
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
z = (ra - rb)/r
! alpha-alpha contribution
if(ra > threshold) then
rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
x = sqrt(rs)
x_f = x*x + b_f*x + c_f
xx0_f = x0_f*x0_f + b_f*x0_f + c_f
q_f = sqrt(4d0*c_f - b_f*b_f)
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(1) = Ec(1) + weight(iG)*ec_f*ra
end if
! alpha-beta contribution
if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
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
Ec(2) = Ec(2) + weight(iG)*ec_z*r
end if
! beta-beta contribution
if(rb > threshold) then
rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
x = sqrt(rs)
x_f = x*x + b_f*x + c_f
xx0_f = x0_f*x0_f + b_f*x0_f + c_f
q_f = sqrt(4d0*c_f - b_f*b_f)
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(3) = Ec(3) + weight(iG)*ec_f*rb
end if
end do
Ec(2) = Ec(2) - Ec(1) - Ec(3)
end subroutine VWN5_lda_correlation_energy

184
src/eDFT/VWN5_lda_correlation_individual_energy.f90
View File

@ -1,184 +0,0 @@
subroutine VWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec)
! Compute VWN5 LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
! Local variables
integer :: iG
integer :: iEns
double precision :: ra,rb,r,raI,rbI,rI,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,dzdrb,dfzdra,dfzdrb,drsdr,decdr_p,decdr_f,decdr_a,decdra,decdrb
double precision :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz
! Output variables
double precision,intent(out) :: LZc(nsp)
double precision,intent(out) :: Ec(nsp,nEns)
! 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
LZc(:) = 0d0
Ec(:,:) = 0d0
do iG=1,nGrid
ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2))
r = ra + rb
! up-down contribution
if(r > threshold) then
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
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
drsdr = - (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*( 2d0/(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*( 2d0/(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*( 2d0/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
decdr_p = drsdr*dxdrs*decdx_p
decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a
dzdra = + (1d0 - z)/r
dfzdra = dzdra*dfzdz
decdra = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+ (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
dzdrb = - (1d0 + z)/r
dfzdrb = dzdrb*dfzdz
decdrb = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+ (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
! spin-up contribution
if(ra > threshold) then
LZc(1) = LZc(1) - weight(iG)*decdra*ra*ra
do iEns=1,nEns
raI = max(0d0,rho(iG,1,iEns))
if(raI > threshold) then
Ec(1,iEns) = Ec(1,iEns) + weight(iG)*(ec_z + decdra*ra)*raI
if(rb > threshold) Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdra*rb*raI
end if
end do
end if
! up-down contribution
if(ra > threshold .and. rb > threshold) LZc(2) = LZc(2) -weight(iG)*(decdra + decdrb)*ra*rb
! spin-down contribution
if(rb > threshold) then
LZc(3) = LZc(3) - weight(iG)*decdrb*rb*rb
do iEns=1,nEns
rbI = max(0d0,rho(iG,2,iEns))
if(rbI > threshold) then
Ec(3,iEns) = Ec(3,iEns) + weight(iG)*(ec_z + decdrb*rb)*rbI
if(ra > threshold) Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdrb*ra*rbI
end if
end do
end if
end if
end do
end subroutine VWN5_lda_correlation_individual_energy

193
src/eDFT/VWN5_lda_correlation_potential.f90
View File

@ -1,193 +0,0 @@
subroutine VWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
! Compute unrestricted VWN5 LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: mu,nu,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 :: Fc(nBas,nBas,nspin)
! 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
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
z = (ra - rb)/r
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
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))
! spin-up contribution
if(ra > threshold) then
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*( 2d0/(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*( 2d0/(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*( 2d0/(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
Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdra*r)
end if
! spin-down contribution
if(rb > threshold) then
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*( 2d0/(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*( 2d0/(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*( 2d0/(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
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
end if
end do
end do
end do
end subroutine VWN5_lda_correlation_potential

52
src/eDFT/W38_lda_correlation_energy.f90
View File

@ -1,52 +0,0 @@
subroutine W38_lda_correlation_energy(nGrid,weight,rho,Ec)
! Compute the unrestricted version of the Wigner's LDA correlation energy
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
! Output variables
double precision :: Ec(nsp)
! Coefficients for Wigner's LDA correlation
a = 0.04918d0
d = 0.349d0
! Compute LDA correlation energy
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
if(ra > threshold .or. rb > threshold) then
r = ra + rb
epsc = ra*rb/(r + d*r**(2d0/3d0))
Ec(2) = Ec(2) + weight(iG)*epsc
endif
enddo
Ec(2) = -4d0*a*Ec(2)
end subroutine W38_lda_correlation_energy

62
src/eDFT/W38_lda_correlation_individual_energy.f90
View File

@ -1,62 +0,0 @@
subroutine W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
! Compute the unrestricted version of the Wigner's LDA individual energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: iG
double precision :: ra,rb,r
double precision :: raI,rbI,rI
double precision :: a,d,epsc
double precision :: dFcdra,dFcdrb
! Output variables
double precision,intent(out) :: Ec(nsp)
! Coefficients for Wigner's LDA correlation
a = 0.04918d0
d = 0.349d0
! Compute LDA correlation individual energy
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2))
raI = max(0d0,rho(iG,1))
rbI = max(0d0,rho(iG,2))
r = ra + rb
rI = raI + rbI
if(r > threshold .or. rI > threshold) then
epsc = ra*rb/(r + d*r**(2d0/3d0))
dFcdra = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/ra)
dFcdrb = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/rb)
Ec(2) = Ec(2) + weight(iG)*rI*0.5d0*(dFcdra + dFcdrb)
endif
enddo
Ec(2) = -4d0*a*Ec(2)
end subroutine W38_lda_correlation_individual_energy

76
src/eDFT/W38_lda_correlation_potential.f90
View File

@ -1,76 +0,0 @@
subroutine W38_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
! Compute the unrestricted version of the Wigner's LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: mu,nu,iG
double precision :: ra,rb,r
double precision :: a,d,ec
double precision :: dFcdr
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Coefficients for Wigner's LDA correlation
a = 0.04918d0
d = 0.349d0
! Compute LDA correlation matrix in the AO basis
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
! Spin-up part contribution
if(ra > threshold) then
r = ra + rb
ec = ra*rb/(r + d*r**(2d0/3d0))
dFcdr = ec*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/ra)
Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdr
endif
! Spin-down part contribution
if(rb > threshold) then
r = ra + rb
ec = ra*rb/(r + d*r**(2d0/3d0))
dFcdr = ec*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/rb)
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdr
endif
enddo
enddo
enddo
Fc(:,:,:) = -4d0*a*Fc(:,:,:)
end subroutine W38_lda_correlation_potential

57
src/eDFT/allocate_grid.f90
View File

@ -1,57 +0,0 @@
subroutine allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid)
! Allocate quadrature grid with numgrid (Radovan Bast)
use numgrid
use, intrinsic :: iso_c_binding, only: c_ptr
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
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 :: radial_precision
integer,intent(in) :: nAng
! Local variables
integer :: iNuc
integer :: min_num_angular_points
integer :: max_num_angular_points
type(c_ptr) :: context
! Output variables
integer,intent(out) :: nGrid
! Set useful variables
min_num_angular_points = nAng
max_num_angular_points = nAng
! Get total number of grid points
nGrid = 0
do iNuc=1,nNuc
context = numgrid_new_atom_grid(radial_precision,min_num_angular_points,max_num_angular_points, &
int(ZNuc(iNuc)),max_exponent(iNuc),max_ang_mom(iNuc), &
min_exponent(iNuc,1:max_ang_mom(iNuc)+1))
nGrid = nGrid + numgrid_get_num_grid_points(context)
call numgrid_free_atom_grid(context)
end do
end subroutine allocate_grid

55
src/eDFT/auxiliary_energy.f90
View File

@ -1,55 +0,0 @@
subroutine auxiliary_energy(nBas,nEns,eps,occnum,Eaux)
! Compute the auxiliary KS energies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables
integer :: iEns,iBas
integer :: ispin
integer :: p
double precision,allocatable :: nEl(:)
! Output variables
double precision,intent(out) :: Eaux(nspin,nEns)
! Memory allocation
allocate(nEl(nEns))
! Compute the number of electrons
nEl(:) = 0d0
do iEns=1,nEns
do iBas=1,nBas
nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
end do
end do
! Compute auxiliary energies for each state of the ensemble based on occupation numbers
Eaux(:,:) = 0d0
do iEns=1,nEns
do ispin=1,nspin
do p=1,nBas
Eaux(ispin,iEns) = Eaux(ispin,iEns) + occnum(p,ispin,iEns)*eps(p,ispin)
end do
end do
end do
end subroutine auxiliary_energy

107
src/eDFT/build_grid.f90
View File

@ -1,107 +0,0 @@
subroutine build_grid(nNuc,ZNuc,rNuc,max_ang_mom,min_exponent,max_exponent, &
radial_precision,nRad,nAng,nGrid,weight,root)
! Compute quadrature grid with numgrid (Radovan Bast)
use numgrid
use, intrinsic :: iso_c_binding, only: c_ptr
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
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) :: radial_precision
integer,intent(in) :: nAng
integer,intent(in) :: nGrid
! Local variables
logical :: dump_grid = .false.
integer :: iNuc
integer :: iG
integer :: min_num_angular_points
integer :: max_num_angular_points
integer :: num_points
integer :: num_radial_points
integer :: center_index
type(c_ptr) :: context
! Output variables
integer,intent(out) :: nRad
double precision,intent(out) :: root(ncart,nGrid)
double precision,intent(out) :: weight(nGrid)
! Set useful variables
min_num_angular_points = nAng
max_num_angular_points = nAng
!------------------------------------------------------------------------
! Main loop over atoms
!------------------------------------------------------------------------
iG = 0
nRad = 0
do iNuc=1,nNuc
context = numgrid_new_atom_grid(radial_precision,min_num_angular_points,max_num_angular_points, &
int(ZNuc(iNuc)),max_exponent(iNuc),max_ang_mom(iNuc), &
min_exponent(iNuc,1:max_ang_mom(iNuc)+1))
center_index = iNuc - 1
num_points = numgrid_get_num_grid_points(context)
num_radial_points = numgrid_get_num_radial_grid_points(context)
call numgrid_get_grid(context,nNuc,center_index,rNuc(:,1),rNuc(:,2),rNuc(:,3),int(ZNuc(:)), &
root(1,iG+1:iG+num_points),root(2,iG+1:iG+num_points),root(3,iG+1:iG+num_points), &
weight(iG+1:iG+num_points))
iG = iG + num_points
nRad = nRad + num_radial_points
call numgrid_free_atom_grid(context)
end do
!------------------------------------------------------------------------
! End main loop over atoms
!------------------------------------------------------------------------
! Print grid
write(*,*)
write(*,'(A30,E10.1)') 'Radial precision = ',radial_precision
write(*,'(A30,I10)') 'Number of radial points = ',nRad
write(*,'(A30,I10)') 'Number of angular points = ',nAng
write(*,'(A30,I10)') 'Total number of points = ',nGrid
write(*,*)
if(dump_grid) then
write(*,*) ' ***********************'
write(*,*) ' *** QUADRATURE GRID ***'
write(*,*) ' ***********************'
write(*,'(A10,3X,3A15)') 'Grid point','X','Y','Z'
do iG=1,nGrid
write(*,'(I10,3X,4F15.10)') iG,weight(iG),root(:,iG)
end do
write(*,*)
end if
end subroutine build_grid

59
src/eDFT/correlation_derivative_discontinuity.f90
View File

@ -1,59 +0,0 @@
subroutine correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
! Compute the correlation part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,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,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
! Local variables
! Output variables
double precision,intent(out) :: Ec(nsp,nEns)
select case (rung)
! Hartree calculation
case(0)
Ec(:,:) = 0d0
! LDA functionals
case(1)
call lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! GGA functionals
case(2)
call gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! MGGA functionals
case(3)
call mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! Hybrid functionals
case(4)
call hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
end select
end subroutine correlation_derivative_discontinuity

59
src/eDFT/correlation_energy.f90
View File

@ -1,59 +0,0 @@
subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
! Compute the unrestricted version of the correlation energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,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
! Output variables
double precision,intent(out) :: Ec(nsp)
select case (rung)
! Hartree calculation
case(0)
Ec(:) = 0d0
! LDA functionals
case(1)
call lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec)
! GGA functionals
case(2)
call gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
! MGGA functionals
case(3)
call mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
! Hybrid functionals
case(4)
call hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
end select
end subroutine correlation_energy

62
src/eDFT/correlation_individual_energy.f90
View File

@ -1,62 +0,0 @@
subroutine correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight, &
rhow,drhow,rho,drho,LZc,Ec)
! Compute the correlation energy of individual states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
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,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns)
! Output variables
double precision,intent(out) :: LZc(nsp)
double precision,intent(out) :: Ec(nsp,nEns)
select case (rung)
! Hartree calculation
case(0)
LZc(:) = 0d0
! LDA functionals
case(1)
call lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec)
! GGA functionals
case(2)
call print_warning('!!! Individual energies NYI for GGAs !!!')
! MGGA functionals
case(3)
call print_warning('!!! Individual energies NYI for MGGAs !!!')
! Hybrid functionals
case(4)
call hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZc,Ec)
end select
end subroutine correlation_individual_energy

68
src/eDFT/correlation_potential.f90
View File

@ -1,68 +0,0 @@
subroutine correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute the correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(ncart,nGrid,nspin)
! Local variables
double precision,allocatable :: FcLDA(:,:,:)
double precision,allocatable :: FcGGA(:,:,:)
double precision :: aC
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Memory allocation
select case (rung)
! Hartree calculation
case(0)
Fc(:,:,:) = 0d0
! LDA functionals
case(1)
call lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc)
! GGA functionals
case(2)
call gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! MGGA functionals
case(3)
call mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Hybrid functionals
case(4)
call hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
end select
end subroutine correlation_potential

38
src/eDFT/density.f90
View File

@ -1,38 +0,0 @@
subroutine density(nGrid,nBas,P,AO,rho)
! Calculate one-electron density
implicit none
include 'parameters.h'
! Input variables
double precision,parameter :: thresh = 1d-15
integer,intent(in) :: nGrid
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
! Local variables
integer :: iG,mu,nu
! Output variables
double precision,intent(out) :: rho(nGrid)
rho(:) = 0d0
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
enddo
enddo
enddo
do iG=1,nGrid
rho(iG) = max(0d0,rho(iG))
enddo
end subroutine density

48
src/eDFT/density_matrix.f90
View File

@ -1,48 +0,0 @@
subroutine density_matrix(nBas,nEns,c,P,occnum)
! Calculate density matrices
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables
integer :: ispin
integer :: iEns
integer :: q
integer :: mu,nu
! Output variables
double precision,intent(out) :: P(nBas,nBas,nspin,nEns)
! Compute density matrix for each state of the ensemble based on occupation numbers
P(:,:,:,:) = 0d0
do iEns=1,nEns
do ispin=1,nspin
do mu=1,nBas
do nu=1,nBas
do q=1,nBas
P(mu,nu,ispin,iEns) = P(mu,nu,ispin,iEns) &
+ occnum(q,ispin,iEns)*c(mu,q,ispin)*c(nu,q,ispin)
end do
end do
end do
end do
end do
end subroutine density_matrix

201
src/eDFT/eDFT.f90
View File

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

69
src/eDFT/elda_correlation_energy.f90
View File

@ -1,69 +0,0 @@
subroutine elda_correlation_energy(aLF,nGrid,weight,rho,Ec)
! Compute LDA correlation energy of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: aLF(3)
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
! Output variables
double precision,intent(out) :: Ec(nsp)
! Compute eLDA correlation energy
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
r = ra + rb
! Spin-up contribution
if(ra > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*ra**(-1d0/6d0) + aLF(3)*ra**(-1d0/3d0))
Ec(1) = Ec(1) + weight(iG)*ec_p*ra
end if
! Opposite-spin contribution
if(r > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
Ec(2) = Ec(2) + weight(iG)*ec_p*r
end if
! Spin-down contribution
if(rb > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*rb**(-1d0/6d0) + aLF(3)*rb**(-1d0/3d0))
Ec(3) = Ec(3) + weight(iG)*ec_p*rb
end if
end do
Ec(2) = Ec(2) - Ec(1) - Ec(3)
end subroutine elda_correlation_energy

57
src/eDFT/elda_correlation_individual_energy.f90
View File

@ -1,57 +0,0 @@
subroutine elda_correlation_individual_energy(aLF,nGrid,weight,rhow,rho,Ec)
! Compute LDA correlation individual energy of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: aLF(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: iG
double precision :: ra, rb, r
double precision :: raI,rbI,rI
double precision :: ec_p,dFcdr
! Output variables
double precision,intent(out) :: Ec(nsp)
! Compute eLDA correlation potential
Ec(:) = 0d0
do iG=1,nGrid
ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2))
raI = max(0d0,rho(iG,1))
rbI = max(0d0,rho(iG,2))
r = ra + rb
rI = raI + rbI
if(r > threshold .or. rI > 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*r)
dFcdr = ec_p + dFcdr*r
Ec(2) = Ec(2) + weight(iG)*rI*dFcdr
end if
end do
end subroutine elda_correlation_individual_energy

70
src/eDFT/elda_correlation_potential.f90
View File

@ -1,70 +0,0 @@
subroutine elda_correlation_potential(aLF,nGrid,weight,nBas,AO,rho,Fc)
! Compute LDA correlation energy of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: aLF(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Local variables
integer :: mu,nu,iG
double precision :: ra,rb,r,ec_p
double precision :: dFcdra,dFcdrb
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Compute eLDA correlation potential
Fc(:,:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
ra = max(0d0,rho(iG,1))
rb = max(0d0,rho(iG,2))
if(ra > threshold) then
r = ra + rb
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*r)
dFcdra = ec_p + dFcdra*r
Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdra
endif
if(rb > threshold) then
r = ra + rb
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*r)
dFcdrb = ec_p + dFcdrb*r
Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdrb
endif
end do
end do
end do
end subroutine elda_correlation_potential

20
src/eDFT/electron_number.f90
View File

@ -1,20 +0,0 @@
function electron_number(nGrid,w,rho) result(nEl)
! Compute the number of electrons via integration of the one-electron density
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: w(nGrid)
double precision,intent(in) :: rho(nGrid)
! Output variables
double precision :: nEl
nEl = dot_product(w,rho)
end function electron_number

67
src/eDFT/exchange_derivative_discontinuity.f90
View File

@ -1,67 +0,0 @@
subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,drhow,&
Cx_choice,doNcentered,kappa,ExDD)
! Compute the exchange part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
!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(:),nCC,aCC,nGrid,weight(:),&
rhow(:),Cx_choice,doNcentered,kappa,ExDD(:))
! GGA functionals
case(2)
call gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:))
! MGGA functionals
case(3)
call mgga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:))
! Hybrid functionals
case(4)
call hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),&
rhow(:),Cx_choice,doNcentered,ExDD(:))
end select
end subroutine exchange_derivative_discontinuity

69
src/eDFT/exchange_energy.f90
View File

@ -1,69 +0,0 @@
subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
rho,drho,Cx_choice,doNcentered,Ex)
! Compute the exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
! Output variables
double precision,intent(out) :: Ex
select case (rung)
! Hartree calculation
case(0)
Ex = 0d0
! LDA functionals
case(1)
call lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
! GGA functionals
case(2)
call gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex)
! MGGA functionals
case(3)
call mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
! Hybrid functionals
case(4)
call hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
rho,drho,Cx_choice,doNcentered,Ex)
end select
end subroutine exchange_energy

71
src/eDFT/exchange_individual_energy.f90
View File

@ -1,71 +0,0 @@
subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
ERI,Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex)
! Compute the exchange individual energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Pw(nBas,nBas,nspin)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: P(nBas,nBas,nspin,nEns)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Output variables
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
select case (rung)
! Hartree calculation
case(0)
Ex = 0d0
! LDA functionals
case(1)
call lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,&
rhow,rho,Cx_choice,doNcentered,LZx,Ex)
! GGA functionals
case(2)
call gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! MGGA functionals
case(3)
call mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! Hybrid functionals
case(4)
call hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,P,rho,drho,LZx,Ex)
end select
end subroutine exchange_individual_energy

80
src/eDFT/exchange_potential.f90
View File

@ -1,80 +0,0 @@
subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF)
! Compute the exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
double precision,allocatable :: FxLDA(:,:)
double precision,allocatable :: FxGGA(:,:)
double precision :: cX,aX
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
double precision,intent(out) :: FxHF(nBas,nBas)
! Memory allocation
select case (rung)
! Hartree calculation
case(0)
Fx(:,:) = 0d0
! LDA functionals
case(1)
call lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,&
Cx_choice,doNcentered,Fx)
! GGA functionals
case(2)
call gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,&
Cx_choice,Fx)
! MGGA functionals
case(3)
call mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Hybrid functionals
case(4)
call hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF)
end select
end subroutine exchange_potential

25
src/eDFT/fock_exchange_energy.f90
View File

@ -1,25 +0,0 @@
subroutine fock_exchange_energy(nBas,P,Fx,Ex)
! Compute the (exact) Fock exchange energy
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: Fx(nBas,nBas)
! Local variables
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: Ex
! Compute HF exchange energy
Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx))
end subroutine fock_exchange_energy

46
src/eDFT/fock_exchange_individual_energy.f90
View File

@ -1,46 +0,0 @@
subroutine fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex)
! Compute the HF individual energy in the unrestricted formalism
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
double precision,intent(in) :: Pw(nBas,nBas,nspin)
double precision,intent(in) :: P(nBas,nBas,nspin,nEns)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: Fx(:,:,:)
double precision,external :: trace_matrix
integer :: iEns
integer :: ispin
! Output variables
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
! Compute HF exchange matrix
allocate(Fx(nBas,nBas,nspin))
do ispin=1,nspin
call fock_exchange_potential(nBas,Pw(:,:,ispin),ERI,Fx(:,:,ispin))
LZx(ispin) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,ispin),Fx(:,:,ispin)))
do iEns=1,nEns
Ex(ispin,iEns) = - 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin,iEns),Fx(:,:,ispin)))
end do
end do
end subroutine fock_exchange_individual_energy

34
src/eDFT/fock_exchange_potential.f90
View File

@ -1,34 +0,0 @@
subroutine fock_exchange_potential(nBas,P,ERI,Fx)
! Compute the Fock exchange potential
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Compute HF exchange matrix
Fx(:,:) = 0d0
do si=1,nBas
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
enddo
enddo
enddo
enddo
end subroutine fock_exchange_potential

32
src/eDFT/generate_shell.f90
View File

@ -1,32 +0,0 @@
subroutine generate_shell(atot,nShellFunction,ShellFunction)
! Generate shells for a given total angular momemtum
implicit none
! Input variables
integer,intent(in) :: atot,nShellFunction
! Local variables
integer :: ax,ay,az,ia
! Output variables
integer,intent(out) :: ShellFunction(nShellFunction,3)
ia = 0
do ax=atot,0,-1
do az=0,atot
ay = atot - ax - az
if(ay >= 0) then
ia = ia + 1
ShellFunction(ia,1) = ax
ShellFunction(ia,2) = ay
ShellFunction(ia,3) = az
endif
enddo
enddo
end subroutine generate_shell

44
src/eDFT/gga_correlation_derivative_discontinuity.f90
View File

@ -1,44 +0,0 @@
subroutine gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! Compute the correlation GGA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
! Local variables
double precision :: aC
! Output variables
double precision,intent(out) :: Ec(nsp,nEns)
! Select correlation functional
select case (DFA)
case (1)
Ec(:,:) = 0d0
case (2)
Ec(:,:) = 0d0
case default
call print_warning('!!! GGA correlation functional not available !!!')
stop
end select
end subroutine gga_correlation_derivative_discontinuity

44
src/eDFT/gga_correlation_energy.f90
View File

@ -1,44 +0,0 @@
subroutine gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
! Compute unrestricted GGA correlation energy
implicit none
include 'parameters.h'
! Input variables
integer,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
integer :: iG
double precision :: ra,rb,ga,gb
! Output variables
double precision :: Ec(nsp)
select case (DFA)
case (1)
call LYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
case (2)
call PBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
case default
call print_warning('!!! GGA correlation energy not available !!!')
stop
end select
end subroutine gga_correlation_energy

46
src/eDFT/gga_correlation_potential.f90
View File

@ -1,46 +0,0 @@
subroutine gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute unrestricted GGA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(3,nGrid,nspin)
! Local variables
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Select GGA exchange functional
select case (DFA)
case (1)
call LYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
case (2)
call PBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
case default
call print_warning('!!! GGA correlation potential not available !!!')
stop
end select
end subroutine gga_correlation_potential

48
src/eDFT/gga_exchange_derivative_discontinuity.f90
View File

@ -1,48 +0,0 @@
subroutine gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD)
! Compute the exchange GGA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,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 exchange functional
select case (DFA)
case (1)
ExDD(:) = 0d0
case (2)
ExDD(:) = 0d0
case (3)
ExDD(:) = 0d0
case default
call print_warning('!!! GGA exchange derivative discontinuity not available !!!')
stop
end select
end subroutine gga_exchange_derivative_discontinuity

53
src/eDFT/gga_exchange_energy.f90
View File

@ -1,53 +0,0 @@
subroutine gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex)
! Select GGA exchange functional for energy calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
double precision,intent(in) :: drho(ncart,nGrid)
! Output variables
double precision :: Ex
select case (DFA)
case (1)
call G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
case (2)
call B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
case (3)
call PBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
case (4)
call CC_B88_gga_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,&
Cx_choice,Ex)
case default
call print_warning('!!! GGA exchange energy not available !!!')
stop
end select
end subroutine gga_exchange_energy

36
src/eDFT/gga_exchange_individual_energy.f90
View File

@ -1,36 +0,0 @@
subroutine gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! Compute GGA exchange energy for individual states
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns)
! Output variables
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional
select case (DFA)
case default
call print_warning('!!! GGA exchange individual energy not available !!!')
stop
end select
end subroutine gga_exchange_individual_energy

57
src/eDFT/gga_exchange_potential.f90
View File

@ -1,57 +0,0 @@
subroutine gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,&
rho,drho,Cx_choice,Fx)
! Select GGA exchange functional for potential calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
integer,intent(in) :: Cx_choice
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Select GGA exchange functional
select case (DFA)
case (1)
call G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
case (2)
call B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
case (3)
call PBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
case (4)
call CC_B88_gga_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,&
Cx_choice,Fx)
case default
call print_warning('!!! GGA exchange potential not available !!!')
stop
end select
end subroutine gga_exchange_potential

45
src/eDFT/gradient_density.f90
View File

@ -1,45 +0,0 @@
subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho)
! Calculate gradient of the one-electron density
implicit none
include 'parameters.h'
! Input variables
double precision,parameter :: thresh = 1d-15
integer,intent(in) :: nGrid
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
! Local variables
integer :: ixyz,iG,mu,nu
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: drho(ncart,nGrid)
drho(:,:) = 0d0
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
do ixyz=1,ncart
drho(ixyz,iG) = drho(ixyz,iG) &
+ P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
enddo
enddo
enddo
enddo
! do iG=1,nGrid
! do ixyz=1,ncart
! if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh
! enddo
! enddo
end subroutine gradient_density

29
src/eDFT/hartree_energy.f90
View File

@ -1,29 +0,0 @@
subroutine hartree_energy(nBas,P,J,EH)
! Compute the unrestricted version of the Hartree energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas,nspin)
double precision,intent(in) :: J(nBas,nBas,nspin)
! Local variables
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: EH(nsp)
! Compute the components of the Hartree energy
EH(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EH(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) &
+ 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1)))
EH(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
end subroutine hartree_energy

55
src/eDFT/hartree_individual_energy.f90
View File

@ -1,55 +0,0 @@
subroutine hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH)
! Compute the hartree contribution to the individual energies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
double precision,intent(in) :: Pw(nBas,nBas,nspin)
double precision,intent(in) :: P(nBas,nBas,nspin,nEns)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: J(:,:,:)
double precision,external :: trace_matrix
integer :: iEns
integer :: ispin
! Output variables
double precision,intent(out) :: LZH(nsp)
double precision,intent(out) :: EH(nsp,nEns)
! Compute HF exchange matrix
allocate(J(nBas,nBas,nspin))
LZH(:) = 0.d0
EH(:,:) = 0.d0
do ispin=1,nspin
call hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin))
end do
LZH(1) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
LZH(2) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
LZH(3) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
do iEns=1,nEns
EH(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1)))
EH(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
+ trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1)))
EH(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2)))
end do
end subroutine hartree_individual_energy

33
src/eDFT/hartree_potential.f90
View File

@ -1,33 +0,0 @@
subroutine hartree_potential(nBas,P,ERI,J)
! Compute the unrestricted version of the Hartree potential
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: J(nBas,nBas)
J(:,:) = 0d0
do si=1,nBas
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
enddo
enddo
enddo
enddo
end subroutine hartree_potential

46
src/eDFT/hybrid_correlation_derivative_discontinuity.f90
View File

@ -1,46 +0,0 @@
subroutine hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! Compute the correlation hybrid part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
! Local variables
! Output variables
double precision,intent(out) :: Ec(nsp,nEns)
! Select correlation functional
select case (DFA)
case (1)
Ec(:,:) = 0d0
case (2)
Ec(:,:) = 0d0
case (3)
Ec(:,:) = 0d0
case default
call print_warning('!!! Hybrid correlation functional not available !!!')
stop
end select
end subroutine hybrid_correlation_derivative_discontinuity

58
src/eDFT/hybrid_correlation_energy.f90
View File

@ -1,58 +0,0 @@
subroutine hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
! Compute the unrestricted version of the correlation energy for hybrid functionals
implicit none
include 'parameters.h'
! Input variables
integer,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 :: EcLDA(nsp)
double precision :: EcGGA(nsp)
double precision :: aC
! Output variables
double precision,intent(out) :: Ec(nsp)
select case (DFA)
case(1)
Ec(:) = 0d0
case(2)
aC = 0.81d0
call lda_correlation_energy(3,nEns,wEns,nGrid,weight,rho,EcLDA)
call gga_correlation_energy(1,nEns,wEns,nGrid,weight,rho,drho,EcGGA)
Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:))
case(3)
call gga_correlation_energy(1,nEns,wEns,nGrid,weight,rho,drho,Ec)
case(4)
call gga_correlation_energy(2,nEns,wEns,nGrid,weight,rho,drho,Ec)
case default
call print_warning('!!! Hybrid correlation energy not available !!!')
stop
end select
end subroutine hybrid_correlation_energy

42
src/eDFT/hybrid_correlation_individual_energy.f90
View File

@ -1,42 +0,0 @@
subroutine hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight, &
rhow,drhow,rho,drho,LZc,Ec)
! Compute the hybrid correlation energy for individual states
implicit none
include 'parameters.h'
! Input variables
integer,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,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nEns)
! Output variables
double precision :: LZc(nsp)
double precision :: Ec(nsp,nEns)
! Select correlation functional
select case (DFA)
case (1)
LZc(:) = 0d0
Ec(:,:) = 0d0
case default
call print_warning('!!! Hybrid correlation individual energy not available !!!')
stop
end select
end subroutine hybrid_correlation_individual_energy

69
src/eDFT/hybrid_correlation_potential.f90
View File

@ -1,69 +0,0 @@
subroutine hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute the correlation potential for hybrid functionals
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(ncart,nGrid,nspin)
! Local variables
double precision,allocatable :: FcLDA(:,:,:)
double precision,allocatable :: FcGGA(:,:,:)
double precision :: aC
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Memory allocation
select case (DFA)
case(1)
Fc(:,:,:) = 0d0
case(2)
allocate(FcLDA(nBas,nBas,nspin),FcGGA(nBas,nBas,nspin))
aC = 0.81d0
call lda_correlation_potential(3,nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA)
call gga_correlation_potential(1,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA)
Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:))
case(3)
allocate(FcGGA(nBas,nBas,nspin))
call gga_correlation_potential(1,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
case(4)
allocate(FcGGA(nBas,nBas,nspin))
call gga_correlation_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
case default
call print_warning('!!! Hybrid correlation potential not available !!!')
stop
end select
end subroutine hybrid_correlation_potential

53
src/eDFT/hybrid_exchange_derivative_discontinuity.f90
View File

@ -1,53 +0,0 @@
subroutine hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
Cx_choice,doNcentered,ExDD)
! Compute the exchange part of the derivative discontinuity for hybrid functionals
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Select exchange functional
select case (DFA)
case (1)
ExDD(:) = 0d0
case (2)
ExDD(:) = 0d0
case (3)
ExDD(:) = 0d0
case default
call print_warning('!!! Hybrid exchange derivative discontinuity not available !!!')
stop
end select
end subroutine hybrid_exchange_derivative_discontinuity

77
src/eDFT/hybrid_exchange_energy.f90
View File

@ -1,77 +0,0 @@
subroutine hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
rho,drho,Cx_choice,doNcentered,Ex)
! Compute the exchange energy for hybrid functionals
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
double precision :: ExLDA,ExGGA,ExHF
double precision :: a0,aX
! Output variables
double precision,intent(out) :: Ex
select case (DFA)
case (1)
call fock_exchange_energy(nBas,P,FxHF,Ex)
case (2)
a0 = 0.20d0
aX = 0.72d0
call lda_exchange_energy(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,&
rho,Cx_choice,doNcentered,ExLDA)
call gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = ExLDA &
+ a0*(ExHF - ExLDA) &
+ aX*(ExGGA - ExLDA)
case (3)
call gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = 0.5d0*ExHF + 0.5d0*ExGGA
case (4)
call gga_exchange_energy(3,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = 0.25d0*ExHF + 0.75d0*ExGGA
case default
call print_warning('!!! Hybrid exchange energy not available !!!')
stop
end select
end subroutine hybrid_exchange_energy

46
src/eDFT/hybrid_exchange_individual_energy.f90
View File

@ -1,46 +0,0 @@
subroutine hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow, &
P,rho,drho,LZx,Ex)
! Compute the hybrid exchange energy for individual states
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns)
integer,intent(in) :: nBas
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Pw(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas,nEns)
! Output variables
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional
select case (DFA)
case (1)
call fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex)
case default
call print_warning('!!! Hybrid exchange individual energy not available !!!')
stop
end select
end subroutine hybrid_exchange_individual_energy

91
src/eDFT/hybrid_exchange_potential.f90
View File

@ -1,91 +0,0 @@
subroutine hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF)
! Compute the exchange potential for hybrid functionals
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
double precision,allocatable :: FxLDA(:,:)
double precision,allocatable :: FxGGA(:,:)
double precision :: a0
double precision :: aX
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
double precision,intent(out) :: FxHF(nBas,nBas)
! Memory allocation
select case (DFA)
case(1)
call fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxHF(:,:)
case(2)
allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas))
a0 = 0.20d0
aX = 0.72d0
call lda_exchange_potential(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight, &
nBas,AO,rho,Cx_choice,doNcentered,FxLDA)
call gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxLDA(:,:) &
+ a0*(FxHF(:,:) - FxLDA(:,:)) &
+ aX*(FxGGA(:,:) - FxLDA(:,:))
case(3)
allocate(FxGGA(nBas,nBas))
call gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = 0.5d0*FxHF(:,:) + 0.5d0*FxGGA(:,:)
case(4)
allocate(FxGGA(nBas,nBas))
call gga_exchange_potential(3,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = 0.25d0*FxHF(:,:) + 0.75d0*FxGGA(:,:)
case default
call print_warning('!!! Hybrid exchange potential not available !!!')
stop
end select
end subroutine hybrid_exchange_potential

241
src/eDFT/individual_energy.f90
View File

@ -1,241 +0,0 @@
subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO, &
T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew)
! Compute unrestricted individual energies as well as excitation energies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: x_rung,c_rung
integer,intent(in) :: x_DFA,c_DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: eKS(nBas,nspin)
double precision,intent(in) :: Pw(nBas,nBas,nspin)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: P(nBas,nBas,nspin,nEns)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns)
double precision,intent(inout):: J(nBas,nBas,nspin)
double precision,intent(inout):: Fx(nBas,nBas,nspin)
double precision,intent(inout):: FxHF(nBas,nBas,nspin)
double precision,intent(inout):: Fc(nBas,nBas,nspin)
double precision,intent(in) :: Ew
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
double precision :: ET(nspin,nEns)
double precision :: EV(nspin,nEns)
double precision :: EH(nsp,nEns)
double precision :: Ex(nspin,nEns)
double precision :: Ec(nsp,nEns)
double precision :: LZH(nsp)
double precision :: LZx(nspin)
double precision :: LZc(nsp)
double precision :: Eaux(nspin,nEns)
double precision :: ExDD(nspin,nEns)
double precision :: EcDD(nsp,nEns)
double precision :: OmH(nEns)
double precision :: Omx(nEns)
double precision :: Omc(nEns)
double precision :: Omaux(nEns)
double precision :: OmxDD(nEns)
double precision :: OmcDD(nEns)
double precision,external :: trace_matrix
integer :: ispin,iEns,iBas
double precision,allocatable :: nEl(:)
double precision,allocatable :: kappa(:)
double precision :: E(nEns)
double precision :: Om(nEns)
double precision,external :: electron_number
! Compute scaling factor for N-centered ensembles
allocate(nEl(nEns),kappa(nEns))
nEl(:) = 0d0
do iEns=1,nEns
do iBas=1,nBas
do ispin=1,nspin
nEl(iEns) = nEl(iEns) + occnum(iBas,ispin,iEns)
end do
end do
kappa(iEns) = nEl(iEns)/nEl(1)
end do
!------------------------------------------------------------------------
! Kinetic energy
!------------------------------------------------------------------------
do ispin=1,nspin
do iEns=1,nEns
ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
end do
end do
!------------------------------------------------------------------------
! Potential energy
!------------------------------------------------------------------------
do iEns=1,nEns
do ispin=1,nspin
EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
end do
end do
!------------------------------------------------------------------------
! Individual Hartree energy
!------------------------------------------------------------------------
LZH(:) = 0d0
EH(:,:) = 0d0
call hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH)
!------------------------------------------------------------------------
! Individual exchange energy
!------------------------------------------------------------------------
LZx(:) = 0d0
Ex(:,:) = 0d0
call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, &
Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex)
!------------------------------------------------------------------------
! Individual correlation energy
!------------------------------------------------------------------------
LZc(:) = 0d0
Ec(:,:) = 0d0
call correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, &
rhow,drhow,rho,drho,LZc,Ec)
!------------------------------------------------------------------------
! Compute auxiliary energies
!------------------------------------------------------------------------
call auxiliary_energy(nBas,nEns,eKS,occnum,Eaux)
!------------------------------------------------------------------------
! Compute derivative discontinuities
!------------------------------------------------------------------------
do ispin=1,nspin
call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nCC,aCC,nGrid,weight, &
rhow(:,ispin),drhow(:,:,ispin),Cx_choice,&
doNcentered,kappa,ExDD(ispin,:))
end do
call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)
! Scaling derivative discontinuity for N-centered ensembles
! if(doNcentered) then
! do iEns=1,nEns
! ExDD(:,iEns) = (1d0 - kappa(iEns))*ExDD(:,iEns)
! EcDD(:,iEns) = (1d0 - kappa(iEns))*EcDD(:,iEns)
! end do
! end if
!------------------------------------------------------------------------
! Total energy
!------------------------------------------------------------------------
if(doNcentered) then
do iEns=1,nEns
E(iEns) = sum(Eaux(:,iEns)) &
+ kappa(iEns)*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))) &
+ sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
end do
else
do iEns=1,nEns
E(iEns) = sum(Eaux(:,iEns)) &
+ sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) &
+ sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
end do
end if
print*,'LZshift=',sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))
! do iEns=1,nEns
! E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) &
! + sum(EH(:,iEns)) + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) &
! + sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) &
! + sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
! end do
!------------------------------------------------------------------------
! Excitation energies
!------------------------------------------------------------------------
do iEns=1,nEns
Om(iEns) = E(iEns) - E(1)
OmH(iEns) = sum(EH(:,iEns)) - sum(EH(:,1))
Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1))
Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1))
Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1))
OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1))
OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1))
end do
if(doNcentered) then
do iEns=1,nEns
OmH(iEns) = OmH(iEns) + (kappa(iEns) - kappa(1))*sum(LZH(:))
Omx(iEns) = Omx(iEns) + (kappa(iEns) - kappa(1))*sum(LZx(:))
Omc(iEns) = Omc(iEns) + (kappa(iEns) - kappa(1))*sum(LZc(:))
Omaux(iEns) = Omaux(iEns) &
+ (kappa(iEns) - kappa(1))*(sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)))
end do
end if
!------------------------------------------------------------------------
! Dump results
!------------------------------------------------------------------------
call print_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, &
Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD)
end subroutine individual_energy

52
src/eDFT/lda_correlation_derivative_discontinuity.f90
View File

@ -1,52 +0,0 @@
subroutine lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! Compute the correlation LDA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
! Local variables
double precision :: aC
! Output variables
double precision,intent(out) :: Ec(nsp,nEns)
! Select correlation functional
select case (DFA)
case (1)
Ec(:,:) = 0d0
case (2)
Ec(:,:) = 0d0
case (3)
Ec(:,:) = 0d0
case (4)
Ec(:,:) = 0d0
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_correlation_derivative_discontinuity

52
src/eDFT/lda_correlation_energy.f90
View File

@ -1,52 +0,0 @@
subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec)
! Select the unrestrited version of the LDA correlation functional
implicit none
include 'parameters.h'
! Input variables
integer,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) :: Ec(nsp)
! Select correlation functional
select case (DFA)
case (1)
call W38_lda_correlation_energy(nGrid,weight,rho,Ec)
case (2)
call PW92_lda_correlation_energy(nGrid,weight,rho,Ec)
case (3)
call VWN3_lda_correlation_energy(nGrid,weight,rho,Ec)
case (4)
call VWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
case (5)
call C16_lda_correlation_energy(nGrid,weight,rho,Ec)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_correlation_energy

51
src/eDFT/lda_correlation_individual_energy.f90
View File

@ -1,51 +0,0 @@
subroutine lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec)
! Compute LDA correlation energy for individual states
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
integer,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,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
! Output variables
double precision :: LZc(nsp)
double precision :: Ec(nsp,nEns)
! Select correlation functional
select case (DFA)
case (1)
! call W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec)
case (2)
! call PW92_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec)
case (3)
! call VWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec)
case (4)
call VWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_correlation_individual_energy

56
src/eDFT/lda_correlation_potential.f90
View File

@ -1,56 +0,0 @@
subroutine lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc)
! Select LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Select correlation functional
select case (DFA)
! Hartree-Fock
case (1)
call W38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
case (2)
call PW92_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
case (3)
call VWN3_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
case (4)
call VWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
case (5)
call C16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine lda_correlation_potential

51
src/eDFT/lda_exchange_derivative_discontinuity.f90
View File

@ -1,51 +0,0 @@
subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
Cx_choice,doNcentered,kappa,ExDD)
! Compute the exchange LDA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa(nEns)
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Select exchange functional
select case (DFA)
case (1)
ExDD(:) = 0d0
case (2)
call CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
Cx_choice,doNcentered,kappa,ExDD)
case default
call print_warning('!!! LDA exchange derivative discontinuity not available !!!')
stop
end select
end subroutine lda_exchange_derivative_discontinuity

46
src/eDFT/lda_exchange_energy.f90
View File

@ -1,46 +0,0 @@
subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
! Select LDA exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Output variables
double precision,intent(out) :: Ex
! Select correlation functional
select case (DFA)
case (1)
call S51_lda_exchange_energy(nGrid,weight,rho,Ex)
case (2)
call CC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
case default
call print_warning('!!! LDA exchange energy not available !!!')
stop
end select
end subroutine lda_exchange_energy

49
src/eDFT/lda_exchange_individual_energy.f90
View File

@ -1,49 +0,0 @@
subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
rho,Cx_choice,doNcentered,LZx,Ex)
! Compute LDA exchange energy for individual states
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Output variables
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional
select case (DFA)
case (1)
call S51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex)
case (2)
call CC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho, &
Cx_choice,doNcentered,LZx,Ex)
case default
call print_warning('!!! LDA exchange individual energy not available !!!')
stop
end select
end subroutine lda_exchange_individual_energy

49
src/eDFT/lda_exchange_potential.f90
View File

@ -1,49 +0,0 @@
subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho, &
Cx_choice,doNcentered,Fx)
! Select LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nCC
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Select exchange functional
select case (DFA)
case (1)
call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
case (2)
call CC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx)
case default
call print_warning('!!! LDA exchange potential not available !!!')
stop
end select
end subroutine lda_exchange_potential

34
src/eDFT/mgga_correlation_derivative_discontinuity.f90
View File

@ -1,34 +0,0 @@
subroutine mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! Compute the correlation MGGA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
! Local variables
! Output variables
double precision,intent(out) :: Ec(nsp,nEns)
! Select correlation functional
select case (DFA)
case default
call print_warning('!!! MGGA correlation functional not available !!!')
stop
end select
end subroutine mgga_correlation_derivative_discontinuity

36
src/eDFT/mgga_correlation_energy.f90
View File

@ -1,36 +0,0 @@
subroutine mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
! Compute unrestricted MGGA correlation energy
implicit none
include 'parameters.h'
! Input variables
integer,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
integer :: iG
double precision :: ra,rb,ga,gb
! Output variables
double precision :: Ec(nsp)
select case (DFA)
case default
call print_warning('!!! MGGA correlation energy not available !!!')
stop
end select
end subroutine mgga_correlation_energy

38
src/eDFT/mgga_correlation_potential.f90
View File

@ -1,38 +0,0 @@
subroutine mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute unrestricted MGGA correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(3,nGrid,nspin)
! Local variables
! Output variables
double precision,intent(out) :: Fc(nBas,nBas,nspin)
! Select MGGA exchange functional
select case (DFA)
case default
call print_warning('!!! MGGA correlation potential not available !!!')
stop
end select
end subroutine mgga_correlation_potential

36
src/eDFT/mgga_exchange_derivative_discontinuity.f90
View File

@ -1,36 +0,0 @@
subroutine mgga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD)
! Compute the exchange MGGA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,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 exchange functional
select case (DFA)
case default
call print_warning('!!! MGGA exchange derivative discontinuity not available !!!')
stop
end select
end subroutine mgga_exchange_derivative_discontinuity

32
src/eDFT/mgga_exchange_energy.f90
View File

@ -1,32 +0,0 @@
subroutine mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
! Select MGGA exchange functional for energy calculation
implicit none
include 'parameters.h'
! Input variables
integer,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)
! Output variables
double precision :: Ex
select case (DFA)
case default
call print_warning('!!! MGGA exchange energy not available !!!')
stop
end select
end subroutine mgga_exchange_energy

36
src/eDFT/mgga_exchange_individual_energy.f90
View File

@ -1,36 +0,0 @@
subroutine mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! Compute MGGA exchange energy for individual states
implicit none
include 'parameters.h'
! Input variables
integer,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,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns)
! Output variables
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional
select case (DFA)
case default
call print_warning('!!! MGGA exchange individual energy not available !!!')
stop
end select
end subroutine mgga_exchange_individual_energy

36
src/eDFT/mgga_exchange_potential.f90
View File

@ -1,36 +0,0 @@
subroutine mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Select MGGA exchange functional for potential calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Select MGGA exchange functional
select case (DFA)
case default
call print_warning('!!! MGGA exchange potential not available !!!')
stop
end select
end subroutine mgga_exchange_potential

1
src/eDFT/obj/.gitignore vendored
View File

@ -1 +0,0 @@
*.o

47
src/eDFT/one_electron_density.f90
View File

@ -1,47 +0,0 @@
subroutine one_electron_density(nGrid,nBas,P,AO,dAO,rho,drho)
! Calculate one-electron density
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
! Local variables
integer :: ixyz,iG,mu,nu
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: rho(nGrid)
double precision,intent(out) :: drho(3,nGrid)
rho(:) = 0d0
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
enddo
enddo
enddo
drho(:,:) = 0d0
do ixyz=1,3
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
drho(ixyz,iG) = drho(ixyz,iG) &
+ P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
enddo
enddo
enddo
enddo
end subroutine one_electron_density

167
src/eDFT/print_UKS.f90
View File

@ -1,167 +0,0 @@
subroutine print_UKS(nBas,nEns,occnum,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
! Print one- and two-electron energies and other stuff for KS calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
double precision,intent(in) :: occnum(nBas,nspin,nEns)
double precision,intent(in) :: Ov(nBas,nBas)
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: eps(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET(nspin)
double precision,intent(in) :: EV(nspin)
double precision,intent(in) :: EH(nsp)
double precision,intent(in) :: Ex(nspin)
double precision,intent(in) :: Ec(nsp)
double precision,intent(in) :: Ew
double precision,intent(in) :: dipole(ncart)
! Local variables
integer :: ixyz
integer :: ispin
integer :: iEns
integer :: iBas
integer :: iHOMOa,iHOMOb
integer :: iLUMOa,iLUMOb
double precision :: HOMOa,HOMOb,HOMO
double precision :: LUMOa,LUMOb,LUMO
double precision :: Gapa,Gapb,Gap
! double precision :: S_exact,S2_exact
! double precision :: S,S2
double precision :: nO(nspin)
! Compute the number of spin-up and spin-down electrons
nO(:) = 0d0
do ispin=1,nspin
do iEns=1,nEns
do iBas=1,nBas
nO(ispin) = nO(ispin) + wEns(iEns)*occnum(iBas,ispin,iEns)
end do
end do
end do
! HOMO and LUMO
iHOMOa = ceiling(nO(1))
iLUMOa = iHOMOa + 1
iHOMOb = ceiling(nO(2))
iLUMOb = iHOMOb + 1
HOMOa = -huge(0d0)
if(iHOMOa > 0) HOMOa = eps(iHOMOa,1)
LUMOa = +huge(0d0)
if(iLUMOa <= nBas) LUMOa = eps(iLUMOa,1)
HOMOb = -huge(0d0)
if(iHOMOb > 0) HOMOb = eps(iHOMOb,2)
LUMOb = +huge(0d0)
if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,2)
HOMO = max(HOMOa,HOMOb)
LUMO = min(LUMOa,LUMOb)
Gapa = LUMOa - HOMOa
Gapb = LUMOb - HOMOb
Gap = LUMO - HOMO
! Spin comtamination
! S2_exact = (nO(1) - nO(2))/2d0*(nO(1) - nO(2))/2d0 + 1d0
! S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
! S_exact = 0.5d0*dble(nO(1) - nO(2))
! S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2)
! Dump results
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40)') ' Summary '
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',sum(ET(:)) + sum(EV(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' One-electron a energy: ',ET(1) + EV(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' One-electron b energy: ',ET(2) + EV(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',sum(ET(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Kinetic a energy: ',ET(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Kinetic b energy: ',ET(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',sum(EV(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron a energy: ',sum(EH(:)) + sum(Ex(:)) + sum(Ec(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EH(1) + Ex(1) + Ec(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EH(2) + Ec(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EH(3) + Ex(2) + Ec(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EH(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EH(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree ab energy: ',EH(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Hartree bb energy: ',EH(3),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au'
write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',Ec(3),' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',Ew,' au'
write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' Kohn-Sham energy: ',Ew + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' KS HOMO a energy:',HOMOa*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS LUMO a energy:',LUMOa*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS HOMOa-LUMOa gap:',Gapa*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' KS HOMO b energy:',HOMOb*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS LUMO b energy:',LUMOb*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS HOMOb-LUMOb gap :',Gapb*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,F13.6,A3)') ' KS HOMO energy:',HOMO*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS LUMO energy:',LUMO*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS HOMO -LUMO gap :',Gap*HatoeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A40,1X,F16.6)') ' S (exact) :',2d0*S_exact + 1d0
! write(*,'(A40,1X,F16.6)') ' S :',2d0*S + 1d0
! write(*,'(A40,1X,F16.6)') ' <S**2> (exact) :',S2_exact
! write(*,'(A40,1X,F16.6)') ' <S**2> :',S2
! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A45)') ' Dipole moment (Debye) '
write(*,'(19X,4A10)') 'X','Y','Z','Tot.'
write(*,'(19X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
! Print results
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'Kohn-Sham spin-up orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c(:,:,1))
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'Kohn-Sham spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c(:,:,2))
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Kohn-Sham spin-up orbital energies '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eps(:,1))
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Kohn-Sham spin-down orbital energies '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eps(:,2))
write(*,*)
end subroutine print_UKS

246
src/eDFT/print_individual_energy.f90
View File

@ -1,246 +0,0 @@
subroutine print_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, &
Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD)
! Print individual energies for eDFT calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: ENuc
double precision,intent(in) :: Ew
double precision,intent(in) :: ET(nspin,nEns)
double precision,intent(in) :: EV(nspin,nEns)
double precision,intent(in) :: EH(nsp,nEns)
double precision,intent(in) :: Ex(nspin,nEns)
double precision,intent(in) :: Ec(nsp,nEns)
double precision,intent(in) :: Eaux(nspin,nEns)
double precision :: LZH(nsp)
double precision :: LZx(nspin)
double precision :: LZc(nsp)
double precision,intent(in) :: ExDD(nspin,nEns)
double precision,intent(in) :: EcDD(nsp,nEns)
double precision,intent(in) :: E(nEns)
double precision,intent(in) :: OmH(nEns)
double precision,intent(in) :: Omx(nEns)
double precision,intent(in) :: Omc(nEns)
double precision,intent(in) :: Omaux(nEns)
double precision,intent(in) :: OmxDD(nEns)
double precision,intent(in) :: OmcDD(nEns)
double precision,intent(in) :: Om(nEns)
! Local variables
integer :: iEns
!------------------------------------------------------------------------
! Ensemble energies
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' ENSEMBLE ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Individual energies
!------------------------------------------------------------------------
! write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A60)') ' INDIVIDUAL TOTAL ENERGIES'
! write(*,'(A60)') '-------------------------------------------------'
! do iEns=1,nEns
! write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
! end do
! write(*,'(A60)') '-------------------------------------------------'
! write(*,*)
!------------------------------------------------------------------------
! Kinetic energy
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Potential energy
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Hartree energy
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EH(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Exchange energy
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Correlation energy
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Auxiliary energies
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' AUXILIARY KS ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') 'Auxiliary KS energy state ',iEns,': ',sum(Eaux(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Print Levy-Zahariev shifts
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' LEVY-ZAHARIEV SHIFTS CONTRIBUTIONS'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:)),' au'
write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:)),' au'
write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:)),' au'
write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZc(:)),' au'
write(*,*)
write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZc(:)))*HaToeV,' eV'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Compute derivative discontinuities
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' ENSEMBLE DERIVATIVE CONTRIBUTIONS'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,*)
write(*,'(A40,I2,A2,F16.10,A3)') ' x ensemble derivative state ',iEns,': ',sum(ExDD(:,iEns)), ' au'
write(*,'(A40,I2,A2,F16.10,A3)') ' c ensemble derivative state ',iEns,': ',sum(EcDD(:,iEns)), ' au'
write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',sum(ExDD(:,iEns))+sum(EcDD(:,iEns)),' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Total Energy and IP and EA
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' ENERGY DIFFERENCES FROM AUXILIARY ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns
write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Omaux(iEns)+OmxDD(iEns)+OmcDD(iEns),' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxDD(iEns)+OmcDD(iEns), ' au'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',(Omaux(iEns)+OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV'
write(*,*)
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
end do
write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns
write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' H energy contribution : ',OmH(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' Hxc energy contribution : ',OmH(iEns)+Omx(iEns)+Omc(iEns), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxDD(iEns)+OmcDD(iEns),' au'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' H energy contribution : ',OmH(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' Hxc energy contribution : ',(OmH(iEns)+Omx(iEns)+Omc(iEns))*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxDD(iEns)+OmcDD(iEns))*HaToeV,' eV'
write(*,*)
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
end subroutine print_individual_energy

49
src/eDFT/read_grid.f90
View File

@ -1,49 +0,0 @@
subroutine read_grid(SGn,radial_precision,nRad,nAng)
! Read grid type
implicit none
! Input variables
integer,intent(in) :: SGn
! Output variables
double precision,intent(out) :: radial_precision
integer,intent(out) :: nRad
integer,intent(out) :: nAng
write(*,*)'----------------------------------------------------------'
write(*,'(A22,I1)')' Quadrature grid: SG-',SGn
write(*,*)'----------------------------------------------------------'
select case (SGn)
case(0)
radial_precision = 1d-5
nRad = 23
nAng = 170
case(1)
radial_precision = 1d-7
nRad = 50
nAng = 194
case(2)
radial_precision = 1d-9
nRad = 75
nAng = 302
case(3)
radial_precision = 1d-11
nRad = 99
nAng = 590
case default
call print_warning('!!! Quadrature grid not available !!!')
stop
end select
end subroutine read_grid

431
src/eDFT/read_options_dft.f90
View File

@ -1,431 +0,0 @@
subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,nCC,aCC, &
doNcentered,occnum,Cx_choice)
! Read DFT options
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
! Local variables
integer :: iBas
integer :: iEns
integer :: iCC
character(len=1) :: answer
double precision,allocatable :: nEl(:)
character(len=12) :: x_func
character(len=12) :: c_func
! Output variables
character(len=8),intent(out) :: method
integer,intent(out) :: x_rung,c_rung
integer,intent(out) :: x_DFA,c_DFA
integer,intent(out) :: SGn
integer,intent(out) :: nEns
logical,intent(out) :: doNcentered
double precision,intent(out) :: wEns(maxEns)
integer,intent(out) :: nCC
double precision,intent(out) :: aCC(maxCC,maxEns-1)
double precision,intent(out) :: occnum(nBas,nspin,maxEns)
integer,intent(out) :: Cx_choice
! Open file with method specification
open(unit=1,file='input/dft')
! Default values
method = 'eDFT-UKS'
x_rung = 1
c_rung = 1
x_DFA = 1
c_DFA = 1
SGn = 0
wEns(:) = 0d0
! Restricted or unrestricted calculation
read(1,*)
read(1,*) method
!---------------------------------------!
! EXCHANGE: read rung of Jacob's ladder !
!---------------------------------------!
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*) x_rung,x_func
select case (x_rung) ! exchange functionals
case (0) ! Hartree
select case (x_func)
case ('H')
x_DFA = 1
case default
call print_warning('!!! Hartree exchange functional not available !!!')
stop
end select
case (1) ! LDA
select case (x_func)
case ('S51')
x_DFA = 1
case ('CC-S51')
x_DFA = 2
case default
call print_warning('!!! LDA exchange functional not available !!!')
stop
end select
case (2) ! GGA
select case (x_func)
case ('G96')
x_DFA = 1
case ('B88')
x_DFA = 2
case ('PBE')
x_DFA = 3
case ('CC-B88')
x_DFA = 4
case default
call print_warning('!!! GGA exchange functional not available !!!')
stop
end select
case (3) ! MGGA
select case (x_func)
case default
call print_warning('!!! MGGA exchange functional not available !!!')
stop
end select
case (4) ! Hybrid
select case (x_func)
case ('HF')
x_DFA = 1
case ('B3LYP')
x_DFA = 2
case ('BHHLYP')
x_DFA = 3
case ('PBE')
x_DFA = 4
case default
call print_warning('!!! Hybrid exchange functional not available !!!')
stop
end select
case default
call print_warning('!!! Exchange rung not available !!!')
stop
end select
! Select rung for exchange
write(*,*)
write(*,*) '*******************************************************************'
write(*,*) '* Exchange rung *'
write(*,*) '*******************************************************************'
call select_rung(x_rung,x_func)
!------------------------------------------!
! CORRELATION: read rung of Jacob's ladder !
!------------------------------------------!
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*) c_rung,c_func
select case (c_rung) ! correlation functionals
case (0) ! Hartree
select case (c_func)
case ('H')
c_DFA = 1
case default
call print_warning('!!! Hartree correlation functional not available !!!')
stop
end select
case (1) ! LDA
select case (c_func)
case ('W38')
c_DFA = 1
case ('PW92')
c_DFA = 2
case ('VWN3')
c_DFA = 3
case ('VWN5')
c_DFA = 4
case ('eVWN5')
c_DFA = 5
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
case (2) ! GGA
select case (c_func)
case ('LYP')
c_DFA = 1
case ('PBE')
c_DFA = 2
case default
call print_warning('!!! GGA correlation functional not available !!!')
stop
end select
case (3) ! MGGA
select case (c_func)
case default
call print_warning('!!! MGGA correlation functional not available !!!')
stop
end select
case (4) ! Hybrid
select case (c_func)
case ('HF')
c_DFA = 1
case ('B3LYP')
c_DFA = 2
case ('BHHLYP')
c_DFA = 3
case ('PBE')
c_DFA = 4
case default
call print_warning('!!! Hybrid correlation functional not available !!!')
stop
end select
case default
call print_warning('!!! Correlation rung not available !!!')
stop
end select
! Select rung for correlation
write(*,*)
write(*,*) '*******************************************************************'
write(*,*) '* Correlation rung *'
write(*,*) '*******************************************************************'
call select_rung(c_rung,c_func)
! Read SG-n grid
read(1,*)
read(1,*) SGn
! Read number of states in ensemble
read(1,*)
read(1,*) nEns
if(nEns.gt.maxEns) then
write(*,*) ' Number of states in ensemble too big!! '
stop
endif
write(*,*)'----------------------------------------------------------'
write(*,'(A33,I3)')' Number of states in ensemble = ',nEns
write(*,*)'----------------------------------------------------------'
write(*,*)
! Read occupation numbers for orbitals nO and nO+1
occnum(:,:,:) = 0d0
do iEns=1,maxEns
read(1,*)
read(1,*) (occnum(iBas,1,iEns),iBas=1,nBas)
read(1,*) (occnum(iBas,2,iEns),iBas=1,nBas)
end do
do iEns=1,nEns
write(*,*)
write(*,*) '==============='
write(*,*) 'State n.',iEns
write(*,*) '==============='
write(*,*)
write(*,*) 'Spin-up occupation numbers'
write(*,*) (int(occnum(iBas,1,iEns)),iBas=1,nBas)
write(*,*) 'Spin-down occupation numbers'
write(*,*) (int(occnum(iBas,2,iEns)),iBas=1,nBas)
write(*,*)
end do
! Read ensemble weights for real physical (fractional number of electrons) ensemble (w1,w2)
allocate(nEl(maxEns))
nEl(:) = 0d0
do iEns=1,maxEns
do iBas=1,nBas
nEl(iEns) = nEl(iEns) + occnum(iBas,1,iEns) + occnum(iBas,2,iEns)
end do
end do
doNcentered = .false.
read(1,*)
read(1,*) (wEns(iEns),iEns=2,nEns)
read(1,*)
read(1,*) answer
if(answer == 'T') doNcentered = .true.
wEns(1) = 1d0
do iEns=2,nEns
wEns(1) = wEns(1) - wEns(iEns)
end do
if (doNcentered) then
do iEns=2,nEns
if(nEl(iEns) > 0d0) then
wEns(iEns) = (nEl(1)/nEl(iEns))*wEns(iEns)
else
wEns(iENs) = 0d0
end if
end do
end if
write(*,*)'----------------------------------------------------------'
write(*,*)' Ensemble weights '
write(*,*)'----------------------------------------------------------'
call matout(nEns,1,wEns)
write(*,*)
! Read parameters for weight-dependent functional
read(1,*)
read(1,*) nCC
do iEns=2,maxEns
read(1,*) (aCC(iCC,iEns-1),iCC=1,nCC)
end do
! Read choice of exchange coefficient
read(1,*)
read(1,*) Cx_choice
write(*,*)'----------------------------------------------------------'
write(*,*)' Parameters for weight-dependent exchange functional '
do iEns=2,maxEns
write(*,*)'----------------------------------------------------------'
write(*,'(A10,I2,A2)') ' State ',iEns,':'
write(*,*)'----------------------------------------------------------'
write(*,*)
call matout(nCC,1,acc(:,iEns-1))
write(*,*)
end do
write(*,*)
! Close file with options
close(unit=1)
end subroutine read_options_dft

49
src/eDFT/select_rung.f90
View File

@ -1,49 +0,0 @@
subroutine select_rung(rung,DFA)
! Select rung of Jacob's ladder
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
select case (rung)
! Hartree calculation
case(0)
write(*,*) "* 0th rung of Jacob's ladder: Hartree calculation *"
! LDA functionals
case(1)
write(*,*) "* 1st rung of Jacob's ladder: local-density approximation (LDA) *"
! GGA functionals
case(2)
write(*,*) "* 2nd rung of Jacob's ladder: generalized gradient approximation (GGA) *"
! meta-GGA functionals
case(3)
write(*,*) "* 3rd rung of Jacob's ladder: meta-GGA functional (MGGA) *"
! Hybrid functionals
case(4)
write(*,*) "* 4th rung of Jacob's ladder: hybrid functional *"
! Default
case default
write(*,*) "!!! rung not available !!!"
stop
end select
! Print selected functional
write(*,*) '* You have selected the following functional: ',DFA,' *'
write(*,*) '*******************************************************************'
write(*,*)
end subroutine select_rung

Some files were not shown because too many files have changed in this diff Show More