10
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00

fix bug in individual energies for VWN

This commit is contained in:
Pierre-Francois Loos 2021-11-01 21:51:02 +01:00
parent 95ae0073ee
commit dc2f8b054c
9 changed files with 274 additions and 254 deletions

View File

@ -17,19 +17,19 @@
# quadrature grid SG-n
0
# Number of states in ensemble (nEns)
4
1
# 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
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
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 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
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
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
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
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
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
1 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
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
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
1 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
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00 0.00 0.00
# N-centered?

View File

@ -1,5 +1,5 @@
# RHF UHF KS MOM
T F F F
F F T F
# MP2* MP3 MP2-F12
F F F
# CCD DCD CCSD CCSD(T)
@ -13,7 +13,7 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW
T F F F F
F F F F F
# G0T0 evGT qsGT
F F F
# MCMP2

View File

@ -54,6 +54,8 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered
enddo
! De-scaling for N-centered ensemble
if(doNcentered) then
Exrr = kappa*Exrr

View File

@ -23,11 +23,14 @@ subroutine UVWN3_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_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
double precision :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
double precision :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
double precision :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
double precision :: dzdr ,dfzdr ,drsdr ,decdr_p ,decdr_f ,decdr_a, decdr
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)
@ -63,24 +66,22 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
! spin-up contribution
if(ra > threshold .or. raI > threshold) then
r = ra
rI = raI
if(r > threshold .or. rI > 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)) ) )
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_f = 2d0*x + b_f
@ -88,19 +89,21 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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 ) )
decdra_f = drsdra*dxdrs*decdx_f
decdr_f = drsdr*dxdrs*decdx_f
Ec(1) = Ec(1) + weight(iG)*(ec_f + decdra_f*r)*rI
Ecrr(1) = Ecrr(1) - weight(iG)*decdr_f*r*r
EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI
EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI
end if
! up-down contribution
if(ra > threshold .or. raI > threshold) then
r = ra + rb
rI = raI + rbI
if(r > threshold .or. rI > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r
x = sqrt(rs)
@ -133,11 +136,11 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
dzdra = (1d0 - z)/r
dzdr = (1d0 - z)/r
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
dfzdra = dzdra*dfzdz
dfzdr = dzdr*dfzdz
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_p = 2d0*x + b_p
@ -153,37 +156,37 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
- b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
decdra_p = drsdra*dxdrs*decdx_p
decdra_f = drsdra*dxdrs*decdx_f
decdra_a = drsdra*dxdrs*decdx_a
decdr_p = drsdr*dxdrs*decdx_p
decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a
decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+ (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
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
Ec(2) = Ec(2) + weight(iG)*(ec_z + decdra*r)*rI
Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r
EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
end if
! spin-down contribution
if(rb > threshold .or. rbI > threshold) then
r = rb
rI = rbI
if(r > threshold .or. rI > 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)) ) )
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_f = 2d0*x + b_f
@ -191,18 +194,29 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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 ) )
decdra_f = drsdra*dxdrs*decdx_f
decdr_f = drsdr*dxdrs*decdx_f
Ec(3) = Ec(3) + weight(iG)*(ec_f + decdra_f*r)*rI
Ecrr(3) = Ecrr(3) - weight(iG)*decdr_f*r*r
EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
end if
end do
! De-scaling for N-centered
Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3)
EcrI(2) = EcrI(2) - EcrI(1) - EcrI(3)
EcrrI(2) = EcrrI(2) - EcrrI(1) - EcrrI(3)
if(doNcentered) Ec(:) = kappa*Ec(:)
! De-scaling for N-centered ensemble
Ec(2) = Ec(2) - Ec(1) - Ec(3)
if(doNcentered) then
Ecrr(:) = kappa*Ecrr(:)
EcrI(:) = kappa*EcrI(:)
endif
Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:)
end subroutine UVWN3_lda_correlation_individual_energy

View File

@ -23,29 +23,32 @@ 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_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
double precision :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
double precision :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
double precision :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
double precision :: dzdr ,dfzdr ,drsdr ,decdr_p ,decdr_f ,decdr_a, decdr
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)
! Parameters of the functional
a_p = +0.0621814D0/2D0
a_p = +0.0621814d0/2d0
x0_p = -0.10498d0
b_p = +3.72744d0
c_p = +12.9352d0
a_f = +0.0621814D0/4D0
a_f = +0.0621814d0/4d0
x0_f = -0.325d0
b_f = +7.06042d0
c_f = +18.0578d0
a_a = -1d0/(6d0*pi**2)
x0_a = -0.0047584D0
x0_a = -0.0047584d0
b_a = +1.13107d0
c_a = +13.0045d0
@ -63,24 +66,22 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
! spin-up contribution
if(ra > threshold .or. raI > threshold) then
r = ra
rI = raI
if(r > threshold .or. rI > 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)) ) )
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_f = 2d0*x + b_f
@ -88,19 +89,21 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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 ) )
decdra_f = drsdra*dxdrs*decdx_f
decdr_f = drsdr*dxdrs*decdx_f
Ec(1) = Ec(1) + weight(iG)*(ec_f + decdra_f*r)*rI
Ecrr(1) = Ecrr(1) - weight(iG)*decdr_f*r*r
EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI
EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI
end if
! up-down contribution
if(ra > threshold .or. raI > threshold) then
r = ra + rb
rI = raI + rbI
if(r > threshold .or. rI > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r
x = sqrt(rs)
@ -133,11 +136,11 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
dzdra = (1d0 - z)/r
dzdr = (1d0 - z)/r
dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
dfzdra = dzdra*dfzdz
dfzdr = dzdr*dfzdz
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_p = 2d0*x + b_p
@ -153,37 +156,37 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
- b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
decdra_p = drsdra*dxdrs*decdx_p
decdra_f = drsdra*dxdrs*decdx_f
decdra_a = drsdra*dxdrs*decdx_a
decdr_p = drsdr*dxdrs*decdx_p
decdr_f = drsdr*dxdrs*decdx_f
decdr_a = drsdr*dxdrs*decdx_a
decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+ (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
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
Ec(2) = Ec(2) + weight(iG)*(ec_z + decdra*r)*rI
Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r
EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
end if
! spin-down contribution
if(rb > threshold .or. rbI > threshold) then
r = rb
rI = rbI
if(r > threshold .or. rI > 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)) ) )
drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
dxdrs = 0.5d0/sqrt(rs)
dxdx_f = 2d0*x + b_f
@ -191,18 +194,29 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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 ) )
decdra_f = drsdra*dxdrs*decdx_f
decdr_f = drsdr*dxdrs*decdx_f
Ec(3) = Ec(3) + weight(iG)*(ec_f + decdra_f*r)*rI
Ecrr(3) = Ecrr(3) - weight(iG)*decdr_f*r*r
EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
end if
end do
! De-scaling for N-centered
Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3)
EcrI(2) = EcrI(2) - EcrI(1) - EcrI(3)
EcrrI(2) = EcrrI(2) - EcrrI(1) - EcrrI(3)
if(doNcentered) Ec(:) = kappa*Ec(:)
! De-scaling for N-centered ensemble
Ec(2) = Ec(2) - Ec(1) - Ec(3)
if(doNcentered) then
Ecrr(:) = kappa*Ecrr(:)
EcrI(:) = kappa*EcrI(:)
endif
Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:)
end subroutine UVWN5_lda_correlation_individual_energy

View File

@ -402,7 +402,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
! Compute individual energies from ensemble energy
!------------------------------------------------------------------------
if(nEns > 1) &
call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om,occnum, &
Cx_choice,doNcentered)

View File

@ -38,6 +38,19 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Individual energies
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' INDIVIDUAL TOTAL ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
end do
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! Kinetic energy
!------------------------------------------------------------------------
@ -136,135 +149,134 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
! Total Energy and IP and EA
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' IP AND EA FROM AUXILIARY ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A60)') ' IP AND EA FROM AUXILIARY ENERGIES '
! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Omaux(2)+OmxcDD(2),' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-(Omaux(3)+OmxcDD(3)),' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',-Omaux(3), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3),' au'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',Omaux(2)+OmxcDD(2)+(Omaux(3)+OmxcDD(3)),' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2)+Omaux(3), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)+OmxDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)+OmcDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)+OmxcDD(3),' au'
write(*,*)
! write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Omaux(2)+OmxcDD(2),' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2), ' au'
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2), ' au'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-(Omaux(3)+OmxcDD(3)),' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',-Omaux(3), ' au'
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3),' au'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',Omaux(2)+OmxcDD(2)+(Omaux(3)+OmxcDD(3)),' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2)+Omaux(3), ' au'
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)+OmxDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)+OmcDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)+OmxcDD(3),' au'
! write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
! write(*,'(A60)') '-------------------------------------------------'
! write(*,*)
write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2:',(Omaux(2)+OmxcDD(2))*HaToeV,' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)*HaToeV,' eV'
write(*,*)
write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3:',-(Omaux(3)+OmxcDD(3))*HaToeV,' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',-Omaux(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3)*HaToeV,' eV'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',(Omaux(2)+OmxcDD(2)+(Omaux(3)+OmxcDD(3)))*HaToeV,' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',(Omaux(2)+Omaux(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',(OmxDD(2)+OmxDD(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',(OmcDD(2)+OmcDD(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxcDD(2)+OmxcDD(3))*HaToeV,' eV'
write(*,*)
! write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2:',(Omaux(2)+OmxcDD(2))*HaToeV,' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)*HaToeV,' eV'
! write(*,*)
! write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3:',-(Omaux(3)+OmxcDD(3))*HaToeV,' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',-Omaux(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3)*HaToeV,' eV'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',(Omaux(2)+OmxcDD(2)+(Omaux(3)+OmxcDD(3)))*HaToeV,' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',(Omaux(2)+Omaux(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',(OmxDD(2)+OmxDD(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',(OmcDD(2)+OmcDD(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxcDD(2)+OmxcDD(3))*HaToeV,' eV'
! write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
! write(*,'(A60)') '-------------------------------------------------'
! write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' IP and EA FROM INDIVIDUAL ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
end do
write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A60)') ' IP and EA FROM INDIVIDUAL ENERGIES '
! write(*,'(A60)') '-------------------------------------------------'
! do iEns=1,nEns
! write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
! end do
! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Om(2), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-Om(3), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',-Omx(3), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',-Omc(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',-Omxc(3), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3),' au'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',Om(2)+Om(3), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2)+Omx(3), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2)+Omc(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2)+Omxc(3), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)+OmxDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)+OmcDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)+OmxcDD(3),' au'
write(*,*)
! write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Om(2), ' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2), ' au'
! write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2), ' au'
! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2), ' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2), ' au'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-Om(3), ' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x energy contribution : ',-Omx(3), ' au'
! write(*,'(A44, F16.10,A3)') ' c energy contribution : ',-Omc(3), ' au'
! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',-Omxc(3), ' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3),' au'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',Om(2)+Om(3), ' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2)+Omx(3), ' au'
! write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2)+Omc(3), ' au'
! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2)+Omxc(3), ' au'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)+OmxDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)+OmcDD(3), ' au'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)+OmxcDD(3),' au'
! write(*,*)
write(*,'(A60)') '-------------------------------------------------'
! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Om(2)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)*HaToeV,' eV'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-Om(3)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',-Omx(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',-Omc(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',-Omxc(3)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3)*HaToeV,' eV'
write(*,*)
write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',(Om(2)+Om(3))*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',(Omx(2)+Omx(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',(Omc(2)+Omc(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',(Omxc(2)+Omxc(3))*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',(OmxDD(2)+OmxDD(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',(OmcDD(2)+OmcDD(3))*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxcDD(2)+OmxcDD(3))*HaToeV,' eV'
write(*,*)
! write(*,'(A43,F16.10,A4)') ' Ionization Potential 1 -> 2:',Om(2)*HaToeV, ' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2)*HaToeV, ' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)*HaToeV,' eV'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Electronic Affinity 1 -> 3:',-Om(3)*HaToeV, ' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x energy contribution : ',-Omx(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c energy contribution : ',-Omc(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',-Omxc(3)*HaToeV, ' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',-OmxDD(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',-OmcDD(3)*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',-OmxcDD(3)*HaToeV,' eV'
! write(*,*)
! write(*,'(A43,F16.10,A4)') ' Fundamental Gap :',(Om(2)+Om(3))*HaToeV, ' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x energy contribution : ',(Omx(2)+Omx(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c energy contribution : ',(Omc(2)+Omc(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',(Omxc(2)+Omxc(3))*HaToeV, ' eV'
! write(*,*)
! write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',(OmxDD(2)+OmxDD(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',(OmcDD(2)+OmcDD(3))*HaToeV, ' eV'
! write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',(OmxcDD(2)+OmxcDD(3))*HaToeV,' eV'
! write(*,*)
!
! write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
! write(*,*)
end subroutine print_unrestricted_individual_energy

View File

@ -102,7 +102,7 @@ rule build_lib
description = Linking $out
"""
LIBS="$LDIR/libxcf90.a $LDIR/libxc.a $LDIR/libnumgrid.a "
LIBS="$LDIR/libnumgrid.a "
rule_build_exe = """
LIBS = {0} $LAPACK $STDCXX
@ -136,25 +136,6 @@ build $LDIR/libnumgrid.a: make_numgrid
generator = true
"""
build_libxc = """
rule make_libxc
command = cd $QUACK_ROOT/libxc-tools ; LIBXC_VERSION="$LIBXC_VERSION" QUACK_ROOT="$QUACK_ROOT" CC="$CC" CXX="$CXX" FC="$FC" ./install_libxc.sh
description = Building libxc
pool = console
rule install_libxc
command = cd $QUACK_ROOT/libxc-tools ; LIBXC_VERSION="$LIBXC_VERSION" QUACK_ROOT="$QUACK_ROOT" CC="$CC" CXX="$CXX" FC="$FC" ./install_libxc.sh install
description = Installing libxc
pool = console
build $QUACK_ROOT/libxc-tools/libxc-$LIBXC_VERSION/src/.libs/libxcf90.a: make_libxc
generator = true
build $LDIR/libxc.a $LDIR/libxcf90.a $IDIR/xc.h $IDIR/xc_funcs_removed.h $IDIR/xc_f90_lib_m.mod $IDIR/xc_funcs_worker.h $IDIR/xc_funcs.h $IDIR/xc_version.h: install_libxc $QUACK_ROOT/libxc-tools/libxc-$LIBXC_VERSION/src/.libs/libxcf90.a
generator = true
"""
build_qcaml = """
rule install_qcaml
command = cd $QUACK_ROOT/qcaml-tools ; ./install_qcaml.sh
@ -189,7 +170,6 @@ build_in_exe_dir = "\n".join([
compiler,
rule_fortran,
rule_build_exe,
build_libxc,
])
build_main = "\n".join([
@ -197,7 +177,6 @@ build_main = "\n".join([
compiler,
rule_git_clone,
build_numgrid,
build_libxc,
build_qcaml,
build_GoDuck,
])