LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-26 15:12:17 +02:00
Code Releases Activity

individual energies are working

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-12-01 10:54:51 +01:00
parent 1b6e7e9e0f
commit 13185c8b64
27 changed files with 401 additions and 392 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
input/dft
View File

@ -6,7 +6,7 @@
# GGA = 2: B88,G96,PBE
# MGGA = 3:
# Hybrid = 4 HF,B3LYP,PBE
1 S51
1 S51
# correlation rung:
# Hartree = 0: H
# LDA = 1: PW92,VWN3,VWN5,eVWN5
@ -19,11 +19,11 @@
# Number of states in ensemble (nEns)
2
# occupation numbers
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
@ -31,7 +31,7 @@
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00 0.0 0.0
0.5 0.0 0.0
# Ncentered ?
F
# Parameters for CC weight-dependent exchange functional

2
src/eDFT/UCC_lda_exchange_energy.f90
View File

@ -103,7 +103,7 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,
r = max(0d0,rho(iG))
if(r > threshold) Ex = Ex + weight(iG)*Cx*r**(4d0/3d0)
if(r > threshold) Ex = Ex + weight(iG)*Cx*r**(1d0/3d0)*r
enddo

2
src/eDFT/UCC_lda_exchange_potential.f90
View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx,Cx_choice,doNcentered)
subroutine UCC_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

8
src/eDFT/US51_lda_exchange_energy.f90
View File

@ -25,13 +25,9 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex)
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
r = max(0d0,rho(iG))
if(r > threshold) then
Ex = Ex + weight(iG)*CxLSDA*r**(4d0/3d0)
endif
if(r > threshold) Ex = Ex + weight(iG)*CxLSDA*r**(1d0/3d0)*r
enddo

46
src/eDFT/US51_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,Ex)
subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex)
! Compute the restricted version of Slater's LDA exchange individual energy
@ -7,34 +7,54 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,Ex)
! Input variables
integer,intent(in) :: nEns
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(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) :: Ex
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
! Compute LDA exchange matrix in the AO basis
LZx(:) = 0d0
Ex(:,:) = 0d0
Ex = 0d0
do ispin=1,nspin
do iG=1,nGrid
do iG=1,nGrid
r = max(0d0,rhow(iG,ispin))
r = max(0d0,rhow(iG))
if(r > threshold) then
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
dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0)
Ex = Ex - weight(iG)*dedr*r*r
endif
enddo
enddo

249
src/eDFT/UVWN3_lda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,Ec)
subroutine UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
! Compute VWN3 LDA correlation potential
@ -8,30 +8,30 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered
! 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 :: dzdr ,dfzdr ,drsdr ,decdr_p ,decdr_f ,decdr_a, decdr
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
double precision :: Ecrr(nsp)
double precision :: EcrI(nsp)
double precision :: EcrrI(nsp)
! Output variables
double precision :: Ec(nsp)
double precision,intent(out) :: LZc(nsp)
double precision,intent(out) :: Ec(nsp,nEns)
! Parameters of the functional
@ -52,187 +52,130 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered
! Initialization
! Ec(:) = 0d0
! Ecrr(:) = 0d0
! EcrI(:) = 0d0
! EcrrI(:) = 0d0
LZc(:) = 0d0
Ec(:,:) = 0d0
! do iG=1,nGrid
do iG=1,nGrid
! ra = max(0d0,rhow(iG,1))
! rb = max(0d0,rhow(iG,2))
ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2))
! raI = max(0d0,rho(iG,1))
! rbI = max(0d0,rho(iG,2))
!
! spin-up contribution
!
! r = ra
! rI = raI
! if(r > threshold) then
! rs = (4d0*pi*r/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)) ) )
!
! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
! dxdrs = 0.5d0/sqrt(rs)
! dxdx_f = 2d0*x + b_f
! 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 ) )
! decdr_f = drsdr*dxdrs*decdx_f
! Ecrr(1) = Ecrr(1) - weight(iG)*decdr_f*r*r
! if(rI > threshold) then
! EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI
! EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI
!
! end if
! end if
r = ra + rb
! up-down contribution
!
! r = ra + rb
! rI = raI + rbI
! if(r > threshold) then
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
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
! dzdr = (1d0 - z)/r
! dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
! dfzdr = dzdr*dfzdz
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)
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
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_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*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
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*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
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
decdr_p = drsdr*dxdrs*decdx_p
decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a
! decdr = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdr/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdr*z**3 &
! + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdr*z**4 + 4d0*(ec_f - ec_p)*fz*dzdr*z**3
dzdra = + (1d0 - z)/r
dfzdra = dzdra*dfzdz
! Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r
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
! if(rI > threshold) then
dzdrb = - (1d0 + z)/r
dfzdrb = dzdrb*dfzdz
! EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
! EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
!
! end if
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
! end if
if(ra > threshold) then
LZc(1) = LZc(1) - weight(iG)*decdra*ra*ra
do iEns=1,nEns
! spin-down contribution
!
! r = rb
! rI = rbI
!
! if(r > threshold) then
raI = max(0d0,rho(iG,1,iEns))
! rs = (4d0*pi*r/3d0)**(-1d0/3d0)
! x = sqrt(rs)
if(raI > threshold) then
! 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(1,iEns) = Ec(1,iEns) + weight(iG)*(ec_z + decdra*ra)*raI
Ec(2,iEns) = Ec(2,iEns) + weight(iG)*decdra*rb*raI
! 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)) ) )
end if
! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
! dxdrs = 0.5d0/sqrt(rs)
end do
! dxdx_f = 2d0*x + b_f
end if
! 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 ) )
! spin-down contribution
! decdr_f = drsdr*dxdrs*decdx_f
if(rb > threshold) then
! Ecrr(3) = Ecrr(3) - weight(iG)*decdr_f*r*r
LZc(3) = LZc(3) - weight(iG)*decdrb*rb*rb
do iEns=1,nEns
! if(rI > threshold) then
rbI = max(0d0,rho(iG,2,iEns))
! EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
! EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
if(rbI > threshold) then
! end if
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 if
! end do
end do
! Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3)
! EcrI(2) = EcrI(2) - EcrI(1) - EcrI(3)
! EcrrI(2) = EcrrI(2) - EcrrI(1) - EcrrI(3)
end if
! De-scaling for N-centered ensemble
end if
! if(doNcentered) then
! Ecrr(:) = kappa*Ecrr(:)
! EcrI(:) = kappa*EcrI(:)
! endif
! Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:)
end do
end subroutine UVWN3_lda_correlation_individual_energy

69
src/eDFT/UVWN5_lda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,LZc)
subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
! Compute VWN5 LDA correlation potential
@ -8,14 +8,17 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered
! 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
@ -27,7 +30,8 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered
! Output variables
double precision :: LZc(nspin)
double precision,intent(out) :: LZc(nsp)
double precision,intent(out) :: Ec(nsp,nEns)
! Parameters of the functional
@ -48,7 +52,8 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered
! Initialization
LZc(:) = 0d0
LZc(:) = 0d0
Ec(:,:) = 0d0
do iG=1,nGrid
@ -115,27 +120,57 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered
decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a
if(ra > threshold) then
dzdra = + (1d0 - z)/r
dfzdra = dzdra*dfzdz
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
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
LZc(1) = LZc(1) - weight(iG)*decdra*ra*r
end if
if(rb > threshold) then
! spin-down contribution
dzdrb = - (1d0 + z)/r
dfzdrb = dzdrb*dfzdz
if(rb > threshold) then
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
LZc(2) = LZc(2) - weight(iG)*decdrb*rb*r
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

50
src/eDFT/eDFT_UKS.f90
View File

@ -1,5 +1,5 @@
subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw,Vxc)
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
@ -16,7 +16,9 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
double precision,intent(in) :: aCC(nCC,nEns-1)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: maxSCF,max_diis,guess_type
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
@ -81,7 +83,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Output variables
double precision,intent(out) :: Ew
double precision,intent(out) :: eps(nBas,nspin)
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)
@ -103,26 +105,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Rung of Jacob's ladder
!------------------------------------------------------------------------
! Select rung for exchange
! write(*,*)
! write(*,*) '*******************************************************************'
! write(*,*) '* Exchange rung *'
! write(*,*) '*******************************************************************'
! call select_rung(x_rung,x_DFA)
! Select rung for correlation
! write(*,*)
! write(*,*) '*******************************************************************'
! write(*,*) '* Correlation rung *'
! write(*,*) '*******************************************************************'
! call select_rung(c_rung,c_DFA)
! Overall rung
xc_rung = max(x_rung,c_rung)
! Memory allocation
@ -144,7 +126,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
do ispin=1,nspin
cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin))
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
@ -254,15 +236,15 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Compute Hartree potential
do ispin=1,nspin
call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin))
call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin))
end do
! Compute exchange potential
do ispin=1,nspin
call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nCC,aCC,nGrid,weight(:),nBas, &
Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
Fx(:,:,ispin),FxHF(:,:,ispin),Cx_choice,doNcentered)
call unrestricted_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
@ -305,7 +287,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
@ -337,9 +319,9 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Exchange energy
do ispin=1,nspin
call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin), &
Cx_choice,doNcentered)
call unrestricted_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
@ -385,7 +367,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Compute final KS energy
call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_UKS(nBas,nEns,nO,S,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
call print_UKS(nBas,nEns,nO,S,wEns,eKS,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipole)
! Compute Vxc for post-HF calculations
@ -396,6 +378,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
!------------------------------------------------------------------------
call unrestricted_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,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,occnum,Cx_choice,doNcentered)
AO,dAO,T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew)
end subroutine eDFT_UKS

3
src/eDFT/print_unrestricted_individual_energy.f90
View File

@ -1,5 +1,4 @@
subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E, &
Om,Omx,Omc,Omaux,OmxDD,OmcDD)
subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,ExDD,EcDD,E,Om,Omx,Omc,Omaux,OmxDD,OmcDD)
! Print individual energies for eDFT calculation

17
src/eDFT/unrestricted_correlation_individual_energy.f90
View File

@ -1,4 +1,5 @@
subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc)
subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight, &
rhow,drhow,rho,drho,doNcentered,LZc,Ec)
! Compute the correlation energy of individual states
@ -16,16 +17,14 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
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)
logical,intent(in) :: doNcentered
! Local variables
double precision :: LZcLDA(nspin)
double precision :: LZcGGA(nspin)
! Output variables
double precision,intent(out) :: LZc(nspin)
double precision,intent(out) :: LZc(nsp)
double precision,intent(out) :: Ec(nsp,nEns)
select case (rung)
@ -39,7 +38,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
case(1)
call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,doNcentered,LZc)
call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
! GGA functionals
@ -57,7 +56,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
case(4)
call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,LZc)
call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,doNcentered,LZc,Ec)
end select

7
src/eDFT/unrestricted_exchange_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
rho,drho,Ex,Cx_choice,doNcentered)
rho,drho,Cx_choice,doNcentered,Ex)
! Compute the exchange energy
@ -43,8 +43,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,
case(1)
call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,&
rho,Ex,Cx_choice,doNcentered)
call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
! GGA functionals
@ -63,7 +62,7 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,
case(4)
call unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
rho,drho,Ex,Cx_choice)
rho,drho,Cx_choice,doNcentered,Ex)
end select

22
src/eDFT/unrestricted_exchange_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
ERI,Pw,rhow,drhow,Cx_choice,doNcentered,Ex)
ERI,Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex)
! Compute the exchange individual energy
@ -19,15 +19,19 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
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)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
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) :: Ex
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
select case (rung)
@ -42,25 +46,25 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
case(1)
call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,&
rhow,Cx_choice,doNcentered,Ex)
rhow,rho,Cx_choice,doNcentered,LZx,Ex)
! GGA functionals
case(2)
call unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ex)
call unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! MGGA functionals
case(3)
call unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ex)
call unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! Hybrid functionals
case(4)
call unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,Ex)
call unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,P,rho,drho,LZx,Ex)
end select

14
src/eDFT/unrestricted_exchange_potential.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice,doNcentered)
ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF)
! Compute the exchange potential
@ -29,12 +29,14 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,a
! Local variables
double precision,allocatable :: FxLDA(:,:),FxGGA(:,:)
double precision,allocatable :: FxLDA(:,:)
double precision,allocatable :: FxGGA(:,:)
double precision :: cX,aX
! Output variables
double precision,intent(out) :: Fx(nBas,nBas),FxHF(nBas,nBas)
double precision,intent(out) :: Fx(nBas,nBas)
double precision,intent(out) :: FxHF(nBas,nBas)
! Memory allocation
@ -50,8 +52,8 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,a
case(1)
call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx,&
Cx_choice,doNcentered)
call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,&
Cx_choice,doNcentered,Fx)
! GGA functionals
@ -70,7 +72,7 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,a
case(4)
call unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice)
ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF)
end select

30
src/eDFT/unrestricted_fock_exchange_individual_energy.f90
View File

@ -1,30 +1,46 @@
subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,ERI,Ex)
subroutine unrestricted_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
double precision,intent(in) :: Pw(nBas,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,allocatable :: Fx(:,:,:)
double precision,external :: trace_matrix
integer :: iEns
integer :: ispin
! Output variables
double precision,intent(out) :: Ex
double precision,intent(out) :: LZx(nspin)
double precision,intent(out) :: Ex(nspin,nEns)
! Compute HF exchange matrix
allocate(Fx(nBas,nBas))
allocate(Fx(nBas,nBas,nspin))
call unrestricted_fock_exchange_potential(nBas,Pw,ERI,Fx)
do ispin=1,nspin
call unrestricted_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
Ex = - 0.5d0*trace_matrix(nBas,matmul(Pw,Fx))
end subroutine unrestricted_fock_exchange_individual_energy

11
src/eDFT/unrestricted_gga_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ex)
subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! Compute GGA exchange energy for individual states
@ -12,12 +12,15 @@ subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weigh
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) :: 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 :: Ex
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional

51
src/eDFT/unrestricted_hartree_individual_energy.f90 Normal file
View File

@ -0,0 +1,51 @@
subroutine unrestricted_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))
do ispin=1,nspin
call unrestricted_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 unrestricted_hartree_individual_energy

12
src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90
View File

@ -1,4 +1,5 @@
subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight, &
rhow,drhow,rho,drho,doNcentered,LZc,Ec)
! Compute the hybrid correlation energy for individual states
@ -14,10 +15,14 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,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)
logical,intent(in) :: doNcentered
! Output variables
double precision :: Ec(nsp)
double precision :: LZc(nsp)
double precision :: Ec(nsp,nEns)
! Select correlation functional
@ -25,7 +30,8 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid
case (1)
Ec(:) = 0d0
LZc(:) = 0d0
Ec(:,:) = 0d0
case default

5
src/eDFT/unrestricted_hybrid_exchange_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, &
rho,drho,Ex,Cx_choice)
rho,drho,Cx_choice,doNcentered,Ex)
! Compute the exchange energy for hybrid functionals
@ -22,6 +22,7 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aC
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Local variables
@ -44,7 +45,7 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aC
aX = 0.72d0
call unrestricted_lda_exchange_energy(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,&
rho,ExLDA,Cx_choice)
rho,Cx_choice,doNcentered,ExLDA)
call unrestricted_gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF)

15
src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90
View File

@ -1,4 +1,5 @@
subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,Ex)
subroutine unrestricted_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
@ -12,16 +13,20 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we
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) :: 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 :: Ex
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional
@ -29,7 +34,7 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we
case (1)
call unrestricted_fock_exchange_individual_energy(nBas,Pw,ERI,Ex)
call unrestricted_fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex)
case default

14
src/eDFT/unrestricted_hybrid_exchange_potential.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice)
ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF)
! Compute the exchange potential for hybrid functionals
@ -24,15 +24,19 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC
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(:,:),FxGGA(:,:)
double precision :: a0,aX
double precision,allocatable :: FxLDA(:,:)
double precision,allocatable :: FxGGA(:,:)
double precision :: a0
double precision :: aX
! Output variables
double precision,intent(out) :: Fx(nBas,nBas),FxHF(nBas,nBas)
double precision,intent(out) :: Fx(nBas,nBas)
double precision,intent(out) :: FxHF(nBas,nBas)
! Memory allocation
@ -51,7 +55,7 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC
aX = 0.72d0
call unrestricted_lda_exchange_potential(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight, &
nBas,AO,rho,FxLDA,Cx_choice)
nBas,AO,rho,Cx_choice,doNcentered,FxLDA)
call unrestricted_gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF)

109
src/eDFT/unrestricted_individual_energy.f90
View File

@ -1,6 +1,5 @@
subroutine unrestricted_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,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,occnum,&
Cx_choice,doNcentered)
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
@ -27,7 +26,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: eps(nBas,nspin)
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)
@ -117,94 +116,33 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Individual Hartree energy
!------------------------------------------------------------------------
do ispin=1,nspin
call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin))
end do
do iEns=1,nEns
if(doNcentered) then
EH(1,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1)))
EH(2,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
+ kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1)))
EH(3,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2)))
else
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 if
end do
! Levy-Zahariev shif for Hartree
if(doNcentered) then
! Do I need to scale this term?
! LZH(:) = 0d0
else
! 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)))
LZH(1) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
LZH(2) = 0d0
! print*,- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1)+Pw(:,:,2),J(:,:,2)))
print*,- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2)))
print*,- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
LZH(3) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
print*,LZH(3)
end if
LZH(:) = 0d0
EH(:,:) = 0d0
call unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH)
!------------------------------------------------------------------------
! Individual exchange energy
!------------------------------------------------------------------------
do iEns=1,nEns
do ispin=1,nspin
call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
P(:,:,ispin,iEns),FxHF(:,:,ispin),rho(:,ispin,iEns),drho(:,:,ispin,iEns), &
Ex(ispin,iEns),Cx_choice,doNcentered)
end do
end do
! Levy-Zahariev shif for exchange
do ispin=1,nspin
call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, &
Pw(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered, &
LZx(ispin))
end do
LZx(:) = 0d0
Ex(:,:) = 0d0
call unrestricted_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
!------------------------------------------------------------------------
do iEns=1,nEns
call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns))
end do
! Levy-Zahariev shif for correlation
call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc)
LZc(:) = 0d0
Ec(:,:) = 0d0
call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, &
rhow,drhow,rho,drho,doNcentered,LZc,Ec)
!------------------------------------------------------------------------
! Compute auxiliary energies
!------------------------------------------------------------------------
call unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux)
call unrestricted_auxiliary_energy(nBas,nEns,eKS,occnum,doNcentered,Eaux)
!------------------------------------------------------------------------
! Compute derivative discontinuities
@ -221,22 +159,21 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Total energy
!------------------------------------------------------------------------
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
print*,E
do iEns=1,nEns
E(iEns) = sum(Eaux(:,iEns)) &
+ sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) &
+ sum(ExDD(:,iEns)) + sum(EcDD(:,iEns))
end do
print*,E
! Alternative way of calculating individual energies
! 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

10
src/eDFT/unrestricted_lda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,doNcentered,LZc)
subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
! Compute LDA correlation energy for individual states
@ -14,11 +14,13 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,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
! Output variables
double precision :: LZc(nspin)
double precision :: LZc(nsp)
double precision :: Ec(nsp,nEns)
! Select correlation functional
@ -34,11 +36,11 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,
case (3)
call UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,LZc)
! call UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
case (4)
call UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,LZc)
call UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec)
case default

2
src/eDFT/unrestricted_lda_exchange_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Ex,Cx_choice,doNcentered)
subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex)
! Select LDA exchange functional

14
src/eDFT/unrestricted_lda_exchange_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
Cx_choice,doNcentered,Ex)
rho,Cx_choice,doNcentered,LZx,Ex)
! Compute LDA exchange energy for individual states
@ -16,14 +16,16 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
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) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin,nEns)
integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered
! Output variables
double precision :: Ex
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional
@ -31,12 +33,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn
case (1)
call US51_lda_exchange_individual_energy(nGrid,weight,rhow,Ex)
call US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex)
case (2)
call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow, &
Cx_choice,doNcentered,Ex)
call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho, &
Cx_choice,doNcentered,LZx,Ex)
case default

6
src/eDFT/unrestricted_lda_exchange_potential.f90
View File

@ -1,5 +1,5 @@
subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx &
,Cx_choice,doNcentered)
subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho, &
Cx_choice,doNcentered,Fx)
! Select LDA correlation potential
@ -37,7 +37,7 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aC
case (2)
call UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Fx,Cx_choice,doNcentered)
call UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx)
case default

2
src/eDFT/unrestricted_mgga_exchange_energy.f90
View File

@ -14,7 +14,7 @@ subroutine unrestricted_mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Output variables

11
src/eDFT/unrestricted_mgga_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ex)
subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex)
! Compute MGGA exchange energy for individual states
@ -12,12 +12,15 @@ subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weig
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) :: 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 :: Ex
double precision :: LZx(nspin)
double precision :: Ex(nspin,nEns)
! Select correlation functional