From 2c2f02d701ef4e8a47535ba971a9fb3aabeed4d8 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 6 Jul 2020 15:42:21 +0200 Subject: [PATCH] clean up --- input/basis | 30 +-- input/weight | 30 +-- .../RB88_gga_exchange_individual_energy.f90 | 2 +- .../RCC_lda_exchange_individual_energy.f90 | 2 +- .../RMFL20_lda_exchange_individual_energy.f90 | 2 +- ...VWN5_lda_correlation_individual_energy.f90 | 2 +- src/eDFT/US51_lda_exchange_energy.f90 | 5 +- .../US51_lda_exchange_individual_energy.f90 | 10 +- src/eDFT/UVWN5_lda_correlation_energy.f90 | 2 +- src/eDFT/UW38_lda_correlation_energy.f90 | 2 +- ...UW38_lda_correlation_individual_energy.f90 | 2 +- ...UWN5_lda_correlation_individual_energy.f90 | 204 ------------------ src/eDFT/eDFT_UKS.f90 | 2 +- .../elda_correlation_individual_energy.f90 | 2 +- src/eDFT/lda_exchange_energy.f90 | 2 +- .../print_unrestricted_individual_energy.f90 | 38 ++-- ...ted_elda_correlation_individual_energy.f90 | 2 +- src/eDFT/unrestricted_correlation_energy.f90 | 8 +- src/eDFT/unrestricted_individual_energy.f90 | 20 +- 19 files changed, 61 insertions(+), 306 deletions(-) delete mode 100644 src/eDFT/UWN5_lda_correlation_individual_energy.f90 diff --git a/input/basis b/input/basis index cd6cb22..6796e3b 100644 --- a/input/basis +++ b/input/basis @@ -1,25 +1,9 @@ -1 10 -S 4 - 1 528.5000000 0.0009400 - 2 79.3100000 0.0072140 - 3 18.0500000 0.0359750 - 4 5.0850000 0.1277820 +1 3 +S 3 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 1.6090000 1.0000000 -S 1 - 1 0.5363000 1.0000000 -S 1 - 1 0.1833000 1.0000000 + 1 0.2976000 1.0000000 P 1 - 1 5.9940000 1.0000000 -P 1 - 1 1.7450000 1.0000000 -P 1 - 1 0.5600000 1.0000000 -D 1 - 1 4.2990000 1.0000000 -D 1 - 1 1.2230000 1.0000000 -F 1 - 1 2.6800000 1.0000000 - + 1 1.2750000 1.0000000 diff --git a/input/weight b/input/weight index cd6cb22..6796e3b 100644 --- a/input/weight +++ b/input/weight @@ -1,25 +1,9 @@ -1 10 -S 4 - 1 528.5000000 0.0009400 - 2 79.3100000 0.0072140 - 3 18.0500000 0.0359750 - 4 5.0850000 0.1277820 +1 3 +S 3 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 1.6090000 1.0000000 -S 1 - 1 0.5363000 1.0000000 -S 1 - 1 0.1833000 1.0000000 + 1 0.2976000 1.0000000 P 1 - 1 5.9940000 1.0000000 -P 1 - 1 1.7450000 1.0000000 -P 1 - 1 0.5600000 1.0000000 -D 1 - 1 4.2990000 1.0000000 -D 1 - 1 1.2230000 1.0000000 -F 1 - 1 2.6800000 1.0000000 - + 1 1.2750000 1.0000000 diff --git a/src/eDFT/RB88_gga_exchange_individual_energy.f90 b/src/eDFT/RB88_gga_exchange_individual_energy.f90 index 761c13e..23c05bb 100644 --- a/src/eDFT/RB88_gga_exchange_individual_energy.f90 +++ b/src/eDFT/RB88_gga_exchange_individual_energy.f90 @@ -40,7 +40,7 @@ subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho, r = max(0d0,0.5d0*rhow(iG)) rI = max(0d0,0.5d0*rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) x = sqrt(g)/r**(4d0/3d0) diff --git a/src/eDFT/RCC_lda_exchange_individual_energy.f90 b/src/eDFT/RCC_lda_exchange_individual_energy.f90 index 641a4c6..66a272e 100644 --- a/src/eDFT/RCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RCC_lda_exchange_individual_energy.f90 @@ -83,7 +83,7 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then e_p = Cx*r**(1d0/3d0) dedr = 1d0/3d0*Cx*r**(-2d0/3d0) diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 index e582ac1..75176b0 100644 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -42,7 +42,7 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then e_p = Cxw*r**(1d0/3d0) dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) diff --git a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 index 8ff4a95..d6fb0a1 100644 --- a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 @@ -42,7 +42,7 @@ subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) x = sqrt(rs) diff --git a/src/eDFT/US51_lda_exchange_energy.f90 b/src/eDFT/US51_lda_exchange_energy.f90 index 920e059..a71cf58 100644 --- a/src/eDFT/US51_lda_exchange_energy.f90 +++ b/src/eDFT/US51_lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex,wEns,nEns) +subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex) ! Compute Slater's LDA exchange energy @@ -11,9 +11,6 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex,wEns,nEns) double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) - ! Local variables integer :: iG diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index 5afc204..139be0f 100644 --- a/src/eDFT/US51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/US51_lda_exchange_individual_energy.f90 @@ -32,12 +32,18 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .or. rI > threshold) then + if(r > threshold) then e = alpha*r**(1d0/3d0) dedr = 1d0/3d0*alpha*r**(-2d0/3d0) - Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) + Ex = Ex - weight(iG)*dedr*r*r + + if(rI > threshold) then + + Ex = Ex + weight(iG)*(e*rI + dedr*r*rI) + + endif endif diff --git a/src/eDFT/UVWN5_lda_correlation_energy.f90 b/src/eDFT/UVWN5_lda_correlation_energy.f90 index 07e8d24..8dadbc4 100644 --- a/src/eDFT/UVWN5_lda_correlation_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_energy.f90 @@ -73,7 +73,7 @@ subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) ! alpha-beta contribution - if(ra > threshold .and. rb > threshold) then + if(ra > threshold .or. rb > threshold) then r = ra + rb rs = (4d0*pi*r/3d0)**(-1d0/3d0) diff --git a/src/eDFT/UW38_lda_correlation_energy.f90 b/src/eDFT/UW38_lda_correlation_energy.f90 index c6f2905..391c8c7 100644 --- a/src/eDFT/UW38_lda_correlation_energy.f90 +++ b/src/eDFT/UW38_lda_correlation_energy.f90 @@ -35,7 +35,7 @@ subroutine UW38_lda_correlation_energy(nGrid,weight,rho,Ec) ra = max(0d0,rho(iG,1)) rb = max(0d0,rho(iG,2)) - if(ra > threshold .and. rb > threshold) then + if(ra > threshold .or. rb > threshold) then r = ra + rb diff --git a/src/eDFT/UW38_lda_correlation_individual_energy.f90 b/src/eDFT/UW38_lda_correlation_individual_energy.f90 index dc00815..6c859c2 100644 --- a/src/eDFT/UW38_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UW38_lda_correlation_individual_energy.f90 @@ -45,7 +45,7 @@ subroutine UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) r = ra + rb rI = raI + rbI - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then epsc = ra*rb/(r + d*r**(2d0/3d0)) dFcdra = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/ra) diff --git a/src/eDFT/UWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UWN5_lda_correlation_individual_energy.f90 deleted file mode 100644 index 03b185f..0000000 --- a/src/eDFT/UWN5_lda_correlation_individual_energy.f90 +++ /dev/null @@ -1,204 +0,0 @@ -subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) - -! Compute VWN5 LDA correlation potential - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nGrid - double precision,intent(in) :: weight(nGrid) - double precision,intent(in) :: rhow(nGrid,nspin) - double precision,intent(in) :: rho(nGrid,nspin) - -! Local variables - - integer :: iG - double precision :: ra,rb,r,raI,rbI,rI,rs,x,z - double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p - double precision :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f - double precision :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a - double precision :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a - double precision :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra - double precision :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb - double precision :: ec_z,ec_p,ec_f,ec_a - double precision :: fz,d2fz - -! Output variables - - double precision :: Ec(nspin) - -! Parameters of the functional - - a_p = +0.0621814D0/2D0 - x0_p = -0.10498d0 - b_p = +3.72744d0 - c_p = +12.9352d0 - - a_f = +0.0621814D0/4D0 - x0_f = -0.325d0 - b_f = +7.06042d0 - c_f = +18.0578d0 - - a_a = -1d0/(6d0*pi**2) - x0_a = -0.0047584D0 - b_a = +1.13107d0 - c_a = +13.0045d0 - -! Initialization - - Ec(:) = 0d0 - - do iG=1,nGrid - - ra = max(0d0,rhow(iG,1)) - rb = max(0d0,rhow(iG,2)) - - raI = max(0d0,rho(iG,1)) - rbI = max(0d0,rho(iG,2)) - -! spin-up contribution - - if(ra > threshold .and. raI > threshold) then - - r = ra + rb - rI = raI + rbI - - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - z = (ra - rb)/r - x = sqrt(rs) - - fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0 - fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0)) - - d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0)) - - x_p = x*x + b_p*x + c_p - x_f = x*x + b_f*x + c_f - x_a = x*x + b_a*x + c_a - - xx0_p = x0_p*x0_p + b_p*x0_p + c_p - xx0_f = x0_f*x0_f + b_f*x0_f + c_f - xx0_a = x0_a*x0_a + b_a*x0_a + c_a - - q_p = sqrt(4d0*c_p - b_p*b_p) - q_f = sqrt(4d0*c_f - b_f*b_f) - q_a = sqrt(4d0*c_a - b_a*b_a) - - ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - - ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & - - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) ) - - ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & - - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) ) - - ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4 - - dzdra = (1d0 - z)/r - dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0)) - dfzdra = dzdra*dfzdz - - drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) - - dxdx_p = 2d0*x + b_p - dxdx_f = 2d0*x + b_f - dxdx_a = 2d0*x + b_a - - decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) - - decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f & - - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) ) - - decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a & - - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) ) - - decdra_p = drsdra*dxdrs*decdx_p - decdra_f = drsdra*dxdrs*decdx_f - decdra_a = drsdra*dxdrs*decdx_a - - decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 & - + (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3 - - Ec(2) = Ec(2) + weight(iG)*(ec_z + decdra*r)*rI - - end if - -! spin-down contribution - - if(rb > threshold .and. rbI > threshold) then - - r = ra + rb - rI = raI + rbI - - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - z = (ra - rb)/r - x = sqrt(rs) - - fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0 - fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0)) - - d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0)) - - x_p = x*x + b_p*x + c_p - x_f = x*x + b_f*x + c_f - x_a = x*x + b_a*x + c_a - - xx0_p = x0_p*x0_p + b_p*x0_p + c_p - xx0_f = x0_f*x0_f + b_f*x0_f + c_f - xx0_a = x0_a*x0_a + b_a*x0_a + c_a - - q_p = sqrt(4d0*c_p - b_p*b_p) - q_f = sqrt(4d0*c_f - b_f*b_f) - q_a = sqrt(4d0*c_a - b_a*b_a) - - ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - - ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & - - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) ) - - ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & - - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) ) - - ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4 - - dzdrb = -(1d0 + z)/r - dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0)) - dfzdrb = dzdrb*dfzdz - - drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) - - dxdx_p = 2d0*x + b_p - dxdx_f = 2d0*x + b_f - dxdx_a = 2d0*x + b_a - - decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) - - decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f & - - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) ) - - decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a & - - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) ) - - decdrb_p = drsdrb*dxdrs*decdx_p - decdrb_f = drsdrb*dxdrs*decdx_f - decdrb_a = drsdrb*dxdrs*decdx_a - - decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 & - + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3 - - Ec(2) = Ec(2) + weight(iG)*(ec_z + decdrb*r)*rI - - end if - - end do - -end subroutine UVWN5_lda_correlation_individual_energy diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index 315359e..e3dd564 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -311,7 +311,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) EJ(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) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) ! Exchange energy diff --git a/src/eDFT/elda_correlation_individual_energy.f90 b/src/eDFT/elda_correlation_individual_energy.f90 index 77efc7e..79000fe 100644 --- a/src/eDFT/elda_correlation_individual_energy.f90 +++ b/src/eDFT/elda_correlation_individual_energy.f90 @@ -39,7 +39,7 @@ subroutine elda_correlation_individual_energy(aLF,nGrid,weight,rhow,rho,Ec) r = ra + rb rI = raI + rbI - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0)) diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 index a1bdfb2..b7141c6 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -27,7 +27,7 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,we case ('US51') - call US51_lda_exchange_energy(nGrid,weight,rho,Ex,wEns,nEns) + call US51_lda_exchange_energy(nGrid,weight,rho,Ex) case ('RS51') diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index 0702d86..c2abf33 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -136,18 +136,18 @@ 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(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',Omaux(2)+OmxcDD(2),' au' + 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(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3 :',Omaux(3)+OmxcDD(3),' au' + 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' @@ -155,17 +155,17 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3),' au' write(*,*) - write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' write(*,*) - write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',(Omaux(2)+OmxcDD(2))*HaToeV,' eV' + 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(*,'(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' @@ -173,18 +173,18 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3)*HaToeV,' eV' write(*,*) - write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' write(*,*) - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' IP and EA FROM INDIVIDUAL ENERGIES ' - 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(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',Om(2), ' au' + 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' @@ -194,7 +194,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au' write(*,*) - write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3 :',Om(3), ' au' + 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' @@ -205,9 +205,9 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3),' au' write(*,*) - write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',Om(2)*HaToeV, ' eV' + 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' @@ -217,7 +217,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, 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 :',Om(3)*HaToeV, ' eV' + 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' @@ -228,7 +228,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3)*HaToeV,' eV' write(*,*) - write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' write(*,*) diff --git a/src/eDFT/restricted_elda_correlation_individual_energy.f90 b/src/eDFT/restricted_elda_correlation_individual_energy.f90 index d2f4d81..cd3212b 100644 --- a/src/eDFT/restricted_elda_correlation_individual_energy.f90 +++ b/src/eDFT/restricted_elda_correlation_individual_energy.f90 @@ -33,7 +33,7 @@ subroutine restricted_elda_correlation_individual_energy(aMFL,nGrid,weight,rhow, r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) - if(r > threshold .and. rI > threshold) then + if(r > threshold .or. rI > threshold) then ec_p = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0)) diff --git a/src/eDFT/unrestricted_correlation_energy.f90 b/src/eDFT/unrestricted_correlation_energy.f90 index 9a49cde..578991f 100644 --- a/src/eDFT/unrestricted_correlation_energy.f90 +++ b/src/eDFT/unrestricted_correlation_energy.f90 @@ -38,13 +38,13 @@ subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,d case(1) - call unrestricted_lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),Ec(:)) + call unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) ! GGA functionals case(2) - call unrestricted_gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),Ec(:)) + call unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Hybrid functionals @@ -52,8 +52,8 @@ subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,d aC = 0.81d0 - call unrestricted_lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),EcLDA(:)) - call unrestricted_gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),EcGGA(:)) + call unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,EcLDA) + call unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,EcGGA) Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index 173fbd9..1a5a34f 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -42,7 +42,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision,intent(in) :: Fc(nBas,nBas,nspin) double precision :: Ew - ! Local variables double precision :: ET(nspin,nEns) @@ -105,26 +104,16 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & - + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & - ! 2.0d0*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))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) - ! if (iEns.ne.2) then EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - ! end if + end do - ! if (nO(2) > 1) then - ! EJ(3,2) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & - ! - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - ! else - ! EJ(3,2) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) - ! end if - - !------------------------------------------------------------------------ ! Checking Hartree contributions for each individual states !------------------------------------------------------------------------ @@ -143,7 +132,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered !------------------------------------------------------------------------ ! Individual exchange energy !------------------------------------------------------------------------ - print*,'old Ex(2,2)=',Ex(2,2) + do iEns=1,nEns do ispin=1,nspin call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, & @@ -151,7 +140,6 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns)) end do end do - print*,'new Ex(2,2)=',Ex(2,2) !------------------------------------------------------------------------ ! Checking exchange contributions for each individual states