LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-07 06:33:55 +01:00
Code Releases Activity

individual energies seems to work

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-11-29 23:32:49 +01:00
parent 851ac5eb46
commit b9cf1fe6d5
17 changed files with 325 additions and 509 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/dft
View File

@ -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
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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) # Ensemble weights: wEns(1),...,wEns(nEns-1)
0.75 0.0 0.0 1.00 0.0 0.0
# Ncentered ? # Ncentered ?
F F
# Parameters for CC weight-dependent exchange functional # Parameters for CC weight-dependent exchange functional

41
src/eDFT/UCC_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,kappa,Ex) subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,Ex)
! Compute the unrestricted version of the curvature-corrected exchange functional ! Compute the unrestricted version of the curvature-corrected exchange functional
@ -14,17 +14,14 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho
integer,intent(in) :: nGrid integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables ! Local variables
integer :: iG integer :: iG
double precision :: r,rI double precision :: r
double precision :: e_p,dedr double precision :: dedr
double precision :: Exrr,ExrI,ExrrI
double precision :: a1,b1,c1,d1,w1 double precision :: a1,b1,c1,d1,w1
double precision :: a2,b2,c2,d2,w2 double precision :: a2,b2,c2,d2,w2
@ -34,11 +31,6 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho
double precision,intent(out) :: Ex double precision,intent(out) :: Ex
! External variable
double precision,external :: electron_number
! Defining enhancements factor for weight-dependent functionals ! Defining enhancements factor for weight-dependent functionals
if(doNcentered) then if(doNcentered) then
@ -105,44 +97,19 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho
! Compute LDA exchange matrix in the AO basis ! Compute LDA exchange matrix in the AO basis
Ex = 0d0 Ex = 0d0
Exrr = 0d0
ExrI = 0d0
ExrrI = 0d0
do iG=1,nGrid do iG=1,nGrid
r = max(0d0,rhow(iG)) r = max(0d0,rhow(iG))
rI = max(0d0,rho(iG))
if(r > threshold) then if(r > threshold) then
e_p = Cx*r**(1d0/3d0)
dedr = 1d0/3d0*Cx*r**(-2d0/3d0) dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
Exrr = Exrr - weight(iG)*dedr*r*r Ex = Ex - weight(iG)*dedr*r*r
if(rI > threshold) then
ExrI = ExrI + weight(iG)*e_p*rI
ExrrI = ExrrI + weight(iG)*dedr*r*rI
endif
endif endif
enddo enddo
! De-scaling for N-centered ensemble
if(doNcentered) then
Exrr = kappa*Exrr
ExrI = kappa*ExrI
endif
Ex = Exrr + ExrI + ExrrI
end subroutine UCC_lda_exchange_individual_energy end subroutine UCC_lda_exchange_individual_energy

37
src/eDFT/US51_lda_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ex) subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,Ex)
! Compute the restricted version of Slater's LDA exchange individual energy ! Compute the restricted version of Slater's LDA exchange individual energy
@ -10,16 +10,12 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
integer,intent(in) :: nGrid integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables ! Local variables
integer :: iG integer :: iG
double precision :: r,rI double precision :: r
double precision :: e,dedr double precision :: dedr
double precision :: Exrr,ExrI,ExrrI
! Output variables ! Output variables
@ -27,42 +23,19 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
! Compute LDA exchange matrix in the AO basis ! Compute LDA exchange matrix in the AO basis
Exrr = 0d0 Ex = 0d0
ExrI = 0d0
ExrrI = 0d0
do iG=1,nGrid do iG=1,nGrid
r = max(0d0,rhow(iG)) r = max(0d0,rhow(iG))
rI = max(0d0,rho(iG))
if(r > threshold) then if(r > threshold) then
e = CxLSDA*r**(1d0/3d0)
dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0) dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0)
Exrr = Exrr - weight(iG)*dedr*r*r Ex = Ex - weight(iG)*dedr*r*r
if(rI > threshold) then
ExrI = ExrI + weight(iG)*e*rI
ExrrI = ExrrI + weight(iG)*dedr*r*rI
endif
endif endif
enddo enddo
! De-scaling for N-centered ensemble
if(doNcentered) then
Exrr = kappa*Exrr
ExrI = kappa*ExrI
endif
Ex = Exrr + ExrI + ExrrI
end subroutine US51_lda_exchange_individual_energy end subroutine US51_lda_exchange_individual_energy

266
src/eDFT/UVWN3_lda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ec) subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,Ec)
! Compute VWN3 LDA correlation potential ! Compute VWN3 LDA correlation potential
@ -11,9 +11,7 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
integer,intent(in) :: nGrid integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin)
logical,intent(in) :: doNcentered logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables ! Local variables
@ -54,187 +52,187 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
! Initialization ! Initialization
Ec(:) = 0d0 ! Ec(:) = 0d0
Ecrr(:) = 0d0 ! Ecrr(:) = 0d0
EcrI(:) = 0d0 ! EcrI(:) = 0d0
EcrrI(:) = 0d0 ! EcrrI(:) = 0d0
do iG=1,nGrid ! do iG=1,nGrid
ra = max(0d0,rhow(iG,1)) ! ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2)) ! rb = max(0d0,rhow(iG,2))
raI = max(0d0,rho(iG,1)) ! raI = max(0d0,rho(iG,1))
rbI = max(0d0,rho(iG,2)) ! rbI = max(0d0,rho(iG,2))
!
! spin-up contribution ! spin-up contribution
!
r = ra ! r = ra
rI = raI ! rI = raI
if(r > threshold) then ! if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) ! rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs) ! x = sqrt(rs)
!
x_f = x*x + b_f*x + c_f ! x_f = x*x + b_f*x + c_f
xx0_f = x0_f*x0_f + b_f*x0_f + c_f ! xx0_f = x0_f*x0_f + b_f*x0_f + c_f
q_f = sqrt(4d0*c_f - b_f*b_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)) & ! 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)) ) ) ! - 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) ! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs) ! dxdrs = 0.5d0/sqrt(rs)
dxdx_f = 2d0*x + b_f ! 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 & ! 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 ) ) ! - 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 ! decdr_f = drsdr*dxdrs*decdx_f
Ecrr(1) = Ecrr(1) - weight(iG)*decdr_f*r*r ! Ecrr(1) = Ecrr(1) - weight(iG)*decdr_f*r*r
if(rI > threshold) then ! if(rI > threshold) then
EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI ! EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI
EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI ! EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI
!
end if ! end if
end if ! end if
! up-down contribution ! up-down contribution
!
r = ra + rb ! r = ra + rb
rI = raI + rbI ! rI = raI + rbI
if(r > threshold) then ! if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) ! rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r ! z = (ra - rb)/r
x = sqrt(rs) ! x = sqrt(rs)
!
fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0 ! fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0)) ! fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
!
d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0)) ! d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
!
x_p = x*x + b_p*x + c_p ! x_p = x*x + b_p*x + c_p
x_f = x*x + b_f*x + c_f ! x_f = x*x + b_f*x + c_f
x_a = x*x + b_a*x + c_a ! x_a = x*x + b_a*x + c_a
!
xx0_p = x0_p*x0_p + b_p*x0_p + c_p ! 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_f = x0_f*x0_f + b_f*x0_f + c_f
xx0_a = x0_a*x0_a + b_a*x0_a + c_a ! xx0_a = x0_a*x0_a + b_a*x0_a + c_a
!
q_p = sqrt(4d0*c_p - b_p*b_p) ! q_p = sqrt(4d0*c_p - b_p*b_p)
q_f = sqrt(4d0*c_f - b_f*b_f) ! q_f = sqrt(4d0*c_f - b_f*b_f)
q_a = sqrt(4d0*c_a - b_a*b_a) ! 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)) & ! 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)) ) ) ! - 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)) & ! 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)) ) ) ! - 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)) & ! 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)) ) ) ! - 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_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
dzdr = (1d0 - z)/r ! dzdr = (1d0 - z)/r
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0)) ! dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
dfzdr = dzdr*dfzdz ! dfzdr = dzdr*dfzdz
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) ! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs) ! dxdrs = 0.5d0/sqrt(rs)
dxdx_p = 2d0*x + b_p ! dxdx_p = 2d0*x + b_p
dxdx_f = 2d0*x + b_f ! dxdx_f = 2d0*x + b_f
dxdx_a = 2d0*x + b_a ! 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 & ! 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 ) ) ! - 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 & ! 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 ) ) ! - 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 & ! 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 ) ) ! - 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 ) )
decdr_p = drsdr*dxdrs*decdx_p ! decdr_p = drsdr*dxdrs*decdx_p
decdr_f = drsdr*dxdrs*decdx_f ! decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a ! 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 = 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 ! + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdr*z**4 + 4d0*(ec_f - ec_p)*fz*dzdr*z**3
Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r ! Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r
if(rI > threshold) then ! if(rI > threshold) then
EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI ! EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI ! EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
!
end if ! end if
end if ! end if
! spin-down contribution ! spin-down contribution
!
r = rb ! r = rb
rI = rbI ! rI = rbI
!
if(r > threshold) then ! if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) ! rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs) ! x = sqrt(rs)
x_f = x*x + b_f*x + c_f ! x_f = x*x + b_f*x + c_f
xx0_f = x0_f*x0_f + b_f*x0_f + c_f ! xx0_f = x0_f*x0_f + b_f*x0_f + c_f
q_f = sqrt(4d0*c_f - b_f*b_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)) & ! 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)) ) ) ! - 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) ! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs) ! dxdrs = 0.5d0/sqrt(rs)
dxdx_f = 2d0*x + b_f ! 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 & ! 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 ) ) ! - 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 ! decdr_f = drsdr*dxdrs*decdx_f
Ecrr(3) = Ecrr(3) - weight(iG)*decdr_f*r*r ! Ecrr(3) = Ecrr(3) - weight(iG)*decdr_f*r*r
if(rI > threshold) then ! if(rI > threshold) then
EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI ! EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI ! EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
end if ! end if
end if ! end if
end do ! end do
Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3) ! Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3)
EcrI(2) = EcrI(2) - EcrI(1) - EcrI(3) ! EcrI(2) = EcrI(2) - EcrI(1) - EcrI(3)
EcrrI(2) = EcrrI(2) - EcrrI(1) - EcrrI(3) ! EcrrI(2) = EcrrI(2) - EcrrI(1) - EcrrI(3)
! De-scaling for N-centered ensemble ! De-scaling for N-centered ensemble
if(doNcentered) then ! if(doNcentered) then
Ecrr(:) = kappa*Ecrr(:) ! Ecrr(:) = kappa*Ecrr(:)
EcrI(:) = kappa*EcrI(:) ! EcrI(:) = kappa*EcrI(:)
endif ! endif
Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:) ! Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:)
end subroutine UVWN3_lda_correlation_individual_energy end subroutine UVWN3_lda_correlation_individual_energy

147
src/eDFT/UVWN5_lda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ec) subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,LZc)
! Compute VWN5 LDA correlation potential ! Compute VWN5 LDA correlation potential
@ -11,9 +11,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
integer,intent(in) :: nGrid integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin)
logical,intent(in) :: doNcentered logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables ! Local variables
@ -23,17 +21,13 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
double precision :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f 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 :: 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 :: 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,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 :: ec_z,ec_p,ec_f,ec_a
double precision :: fz,d2fz double precision :: fz,d2fz
double precision :: Ecrr(nsp)
double precision :: EcrI(nsp)
double precision :: EcrrI(nsp)
! Output variables ! Output variables
double precision :: Ec(nsp) double precision :: LZc(nspin)
! Parameters of the functional ! Parameters of the functional
@ -54,56 +48,14 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
! Initialization ! Initialization
Ec(:) = 0d0 LZc(:) = 0d0
Ecrr(:) = 0d0
EcrI(:) = 0d0
EcrrI(:) = 0d0
do iG=1,nGrid do iG=1,nGrid
ra = max(0d0,rhow(iG,1)) ra = max(0d0,rhow(iG,1))
rb = max(0d0,rhow(iG,2)) rb = max(0d0,rhow(iG,2))
raI = max(0d0,rho(iG,1))
rbI = max(0d0,rho(iG,2))
r = ra + rb r = ra + rb
rI = raI + rbI
! spin-up contribution
! 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
! up-down contribution ! up-down contribution
@ -143,12 +95,6 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0)) dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
dzdra = + (1d0 - z)/r
dfzdra = dzdra*dfzdz
dzdrb = - (1d0 + z)/r
dfzdrb = dzdrb*dfzdz
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs) dxdrs = 0.5d0/sqrt(rs)
@ -169,77 +115,32 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
decdr_f = drsdr*dxdrs*decdx_f decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a decdr_a = drsdr*dxdrs*decdx_a
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 & if(ra > threshold) then
+ (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
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 & dzdra = + (1d0 - z)/r
+ (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3 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
LZc(1) = LZc(1) - weight(iG)*decdra*ra*r
end if
if(rb > threshold) then
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
decdr = 0d0 LZc(2) = LZc(2) - weight(iG)*decdrb*rb*r
if(ra > threshold) decdr = decdr + decdra
if(rb > threshold) decdr = decdr + decdrb
Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r
if(rI > threshold) then
EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
end if end if
end if end if
! spin-down contribution
! 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(3) = Ecrr(3) - weight(iG)*decdr_f*r*r
! if(rI > threshold) then
! EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
! EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
! 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)
! De-scaling for N-centered ensemble
if(doNcentered) then
Ecrr(:) = kappa*Ecrr(:)
EcrI(:) = kappa*EcrI(:)
endif
Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:)
end subroutine UVWN5_lda_correlation_individual_energy end subroutine UVWN5_lda_correlation_individual_energy

16
src/eDFT/eDFT_UKS.f90
View File

@ -76,9 +76,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
double precision,allocatable :: rho(:,:,:) double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: drho(:,:,:,:) double precision,allocatable :: drho(:,:,:,:)
double precision :: E(nEns)
double precision :: Om(nEns)
integer :: ispin,iEns,iBas integer :: ispin,iEns,iBas
! Output variables ! Output variables
@ -382,14 +379,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
end if end if
!!!!!
! do iEns=1,nEns
! print*,'occnum=',occnum(1,1,iEns),occnum(2,1,iEns),occnum(1,2,iEns),occnum(2,2,iEns)
! print*,'nel up and down and total=', electron_number(nGrid,weight,&
! rho(:,1,iEns)),electron_number(nGrid,weight,rho(:,2,iEns)),sum(nEl(:))
! end do
!!!!!
! Compute final KS energy ! Compute final KS energy
call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
@ -403,8 +392,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Compute individual energies from ensemble energy ! Compute individual energies from ensemble energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & 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,E,Om,occnum, & AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,occnum,Cx_choice,doNcentered)
Cx_choice,doNcentered)
end subroutine eDFT_UKS end subroutine eDFT_UKS

106
src/eDFT/print_unrestricted_individual_energy.f90
View File

@ -34,7 +34,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' ENSEMBLE ENERGIES' write(*,'(A60)') ' ENSEMBLE ENERGIES'
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au' write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
write(*,*) write(*,*)
@ -55,66 +55,66 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
! Kinetic energy ! Kinetic energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES' ! write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES'
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns ! do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au' ! write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au'
end do ! end do
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,*) ! write(*,*)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Potential energy ! Potential energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES' ! write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES'
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns ! do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au' ! write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au'
end do ! end do
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,*) ! write(*,*)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Hartree energy ! Hartree energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES' ! write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES'
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns ! do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EJ(:,iEns)),' au' ! write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EJ(:,iEns)),' au'
end do ! end do
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,*) ! write(*,*)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Exchange energy ! Exchange energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES' ! write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES'
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns ! do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au' ! write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au'
end do ! end do
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,*) ! write(*,*)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Correlation energy ! Correlation energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES' ! write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES'
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns ! do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au' ! write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au'
end do ! end do
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
write(*,*) ! write(*,*)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Auxiliary energies ! Auxiliary energies
@ -179,18 +179,18 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL ENERGIES ' write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL ENERGIES '
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns ! do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' ! write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
end do ! end do
write(*,'(A60)') '-------------------------------------------------' ! write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns do iEns=2,nEns
write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au' write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au'
write(*,*) write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(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)') ' c energy contribution : ',Omc(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns), ' au' ! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns), ' au'
write(*,*) ! write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(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)') ' c ensemble derivative : ',OmcDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iENs),' au' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iENs),' au'
@ -200,10 +200,10 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns)*HaToeV, ' eV' write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns)*HaToeV, ' eV'
write(*,*) write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(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)') ' c energy contribution : ',Omc(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV' ! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV'
write(*,*) ! write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(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)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV'

19
src/eDFT/unrestricted_correlation_individual_energy.f90
View File

@ -1,5 +1,4 @@
subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho, & subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc)
doNcentered,kappa,Ec)
! Compute the correlation energy of individual states ! Compute the correlation energy of individual states
@ -17,19 +16,16 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: drhow(ncart,nGrid,nspin) double precision,intent(in) :: drhow(ncart,nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin)
double precision,intent(in) :: drho(ncart,nGrid,nspin)
logical,intent(in) :: doNcentered logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Local variables ! Local variables
double precision :: EcLDA(nsp) double precision :: LZcLDA(nspin)
double precision :: EcGGA(nsp) double precision :: LZcGGA(nspin)
! Output variables ! Output variables
double precision,intent(out) :: Ec(nsp) double precision,intent(out) :: LZc(nspin)
select case (rung) select case (rung)
@ -37,14 +33,13 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
case(0) case(0)
Ec(:) = 0d0 LZc(:) = 0d0
! LDA functionals ! LDA functionals
case(1) case(1)
call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho, & call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,doNcentered,LZc)
doNcentered,kappa,Ec)
! GGA functionals ! GGA functionals
@ -62,7 +57,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
case(4) case(4)
call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,LZc)
end select end select

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

6
src/eDFT/unrestricted_fock_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,ERI,Ex)
! Compute the HF individual energy in the unrestricted formalism ! Compute the HF individual energy in the unrestricted formalism
@ -8,7 +8,6 @@ subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex)
integer,intent(in) :: nBas integer,intent(in) :: nBas
double precision,intent(in) :: Pw(nBas,nBas) double precision,intent(in) :: Pw(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
@ -26,7 +25,6 @@ subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex)
call unrestricted_fock_exchange_potential(nBas,Pw,ERI,Fx) call unrestricted_fock_exchange_potential(nBas,Pw,ERI,Fx)
Ex = trace_matrix(nBas,matmul(P ,Fx)) & Ex = - 0.5d0*trace_matrix(nBas,matmul(Pw,Fx))
- 0.5d0*trace_matrix(nBas,matmul(Pw,Fx))
end subroutine unrestricted_fock_exchange_individual_energy end subroutine unrestricted_fock_exchange_individual_energy

4
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,rho,drho,Ex) subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ex)
! Compute GGA exchange energy for individual states ! Compute GGA exchange energy for individual states
@ -14,8 +14,6 @@ subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weigh
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid) double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Output variables ! Output variables

4
src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
! Compute the hybrid correlation energy for individual states ! Compute the hybrid correlation energy for individual states
@ -14,8 +14,6 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid) double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Output variables ! Output variables

7
src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,P,rhow,drhow,rho,drho,Ex) subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,Ex)
! Compute the hybrid exchange energy for individual states ! Compute the hybrid exchange energy for individual states
@ -14,13 +14,10 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid) double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: nBas integer,intent(in) :: nBas
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Pw(nBas,nBas) double precision,intent(in) :: Pw(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
! Output variables ! Output variables
@ -32,7 +29,7 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we
case (1) case (1)
call unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) call unrestricted_fock_exchange_individual_energy(nBas,Pw,ERI,Ex)
case default case default

138
src/eDFT/unrestricted_individual_energy.f90
View File

@ -1,5 +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, & 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,E,Om,occnum,& T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,occnum,&
Cx_choice,doNcentered) Cx_choice,doNcentered)
! Compute unrestricted individual energies as well as excitation energies ! Compute unrestricted individual energies as well as excitation energies
@ -40,7 +40,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,intent(in) :: Fx(nBas,nBas,nspin) double precision,intent(in) :: Fx(nBas,nBas,nspin)
double precision,intent(in) :: FxHF(nBas,nBas,nspin) double precision,intent(in) :: FxHF(nBas,nBas,nspin)
double precision,intent(in) :: Fc(nBas,nBas,nspin) double precision,intent(in) :: Fc(nBas,nBas,nspin)
double precision :: Ew double precision,intent(in) :: Ew
double precision,intent(in) :: occnum(nBas,nspin,nEns) double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice integer,intent(in) :: Cx_choice
logical,intent(in) :: doNcentered logical,intent(in) :: doNcentered
@ -54,6 +54,10 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision :: Ex(nspin,nEns) double precision :: Ex(nspin,nEns)
double precision :: Ec(nsp,nEns) double precision :: Ec(nsp,nEns)
double precision :: Exc(nEns) double precision :: Exc(nEns)
double precision :: LZH(nspin)
double precision :: LZx(nspin)
double precision :: LZc(nspin)
double precision :: LZHxc(nspin)
double precision :: Eaux(nspin,nEns) double precision :: Eaux(nspin,nEns)
double precision :: ExDD(nspin,nEns) double precision :: ExDD(nspin,nEns)
@ -70,15 +74,13 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
double precision,allocatable :: nEl(:) double precision,allocatable :: nEl(:)
double precision,allocatable :: kappa(:) double precision,allocatable :: kappa(:)
double precision :: E(nEns)
double precision :: Om(nEns)
double precision,external :: electron_number double precision,external :: electron_number
! Output variables ! Compute scaling factor for N-centered ensembles
double precision,intent(out) :: E(nEns)
double precision,intent(out) :: Om(nEns)
allocate(nEl(nEns),kappa(nEns)) allocate(nEl(nEns),kappa(nEns))
nEl(:) = 0d0 nEl(:) = 0d0
@ -95,83 +97,91 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Kinetic energy ! Kinetic energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do ispin=1,nspin ! do ispin=1,nspin
do iEns=1,nEns ! do iEns=1,nEns
ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) ! ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
end do ! end do
end do ! end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Potential energy ! Potential energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do iEns=1,nEns ! do iEns=1,nEns
do ispin=1,nspin ! do ispin=1,nspin
EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) ! EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
end do ! end do
end do ! end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Individual Hartree energy ! Individual Hartree energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do iEns=1,nEns ! do iEns=1,nEns
do ispin=1,nspin ! do ispin=1,nspin
call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) ! call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin))
end do ! end do
if(doNcentered) then ! if(doNcentered) then
!
EJ(1,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & ! EJ(1,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) &
- 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) ! - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
!
EJ(2,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & ! EJ(2,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
+ kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & ! + kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) &
- 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & ! - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) &
- 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) ! - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
!
EJ(3,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & ! EJ(3,iEns) = kappa(iEns)*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
- 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) ! - 0.5d0*kappa(iEns)*kappa(iEns)*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
else ! else
EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & ! EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) ! - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & ! EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) &
+ trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & ! + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & ! - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) ! - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1)))
EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & ! EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) ! - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
end if ! end if
! end do
!------------------------------------------------------------------------
! Individual Hartree energy
!------------------------------------------------------------------------
do ispin=1,nspin
LZH(ispin) = -0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1)+Pw(:,:,2),J(:,:,ispin)))
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Individual exchange energy ! Individual exchange energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do iEns=1,nEns do ispin=1,nspin
do ispin=1,nspin call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, &
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, &
Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & LZx(ispin))
rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,kappa(iEns), &
Ex(ispin,iEns))
end do
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Individual correlation energy ! Individual correlation energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
do iEns=1,nEns call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc)
call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, &
rhow,drhow,rho(:,:,iEns),drho(:,:,:,iEns),doNcentered,kappa(iEns),Ec(:,iEns)) !------------------------------------------------------------------------
end do ! Individual exchange-correlation energy
!------------------------------------------------------------------------
print*,LZc
LZHxc(:) = LZH(:) + LZx(:) + LZc(:)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Compute auxiliary energies ! Compute auxiliary energies
@ -199,10 +209,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
! Total energy ! Total energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! do iEns=1,nEns
! Exc(iEns) = sum(Ex(:,iEns)) + sum(Ec(:,iEns))
! E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) &
! + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(ExcDD(:,iEns))
! end do
do iEns=1,nEns do iEns=1,nEns
Exc(iEns) = sum(Ex(:,iEns)) + sum(Ec(:,iEns)) E(iEns) = sum(Eaux(:,iEns)) + sum(LZHxc(:)) + sum(ExcDD(:,iEns))
E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) &
+ sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(ExcDD(:,iEns))
end do end do
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -212,9 +226,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
do iEns=1,nEns do iEns=1,nEns
Om(iEns) = E(iEns) - E(1) Om(iEns) = E(iEns) - E(1)
Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1)) ! Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1))
Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1)) ! Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1))
Omxc(iEns) = Exc(iEns) - Exc(1) ! Omxc(iEns) = Exc(iEns) - Exc(1)
Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1)) Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1))

11
src/eDFT/unrestricted_lda_correlation_individual_energy.f90
View File

@ -1,5 +1,4 @@
subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho, & subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,doNcentered,LZc)
doNcentered,kappa,Ec)
! Compute LDA correlation energy for individual states ! Compute LDA correlation energy for individual states
@ -15,13 +14,11 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,
integer,intent(in) :: nGrid integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin) double precision,intent(in) :: rhow(nGrid,nspin)
double precision,intent(in) :: rho(nGrid,nspin)
logical,intent(in) :: doNcentered logical,intent(in) :: doNcentered
double precision,intent(in) :: kappa
! Output variables ! Output variables
double precision :: Ec(nsp) double precision :: LZc(nspin)
! Select correlation functional ! Select correlation functional
@ -37,11 +34,11 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,
case (3) case (3)
call UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ec) call UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,LZc)
case (4) case (4)
call UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcentered,kappa,Ec) call UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,doNcentered,LZc)
case default case default

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

4
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,rho,drho,Ex) subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ex)
! Compute MGGA exchange energy for individual states ! Compute MGGA exchange energy for individual states
@ -14,8 +14,6 @@ subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weig
double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid) double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Output variables ! Output variables