diff --git a/input/dft b/input/dft index 8c87da5..0a32c89 100644 --- a/input/dft +++ b/input/dft @@ -31,7 +31,7 @@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.75 0.0 0.0 + 1.00 0.0 0.0 # Ncentered ? F # Parameters for CC weight-dependent exchange functional diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90 index 9977288..01e0dfb 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -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 @@ -14,17 +14,14 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) - double precision,intent(in) :: rho(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! Local variables integer :: iG - double precision :: r,rI - double precision :: e_p,dedr - double precision :: Exrr,ExrI,ExrrI + double precision :: r + double precision :: dedr double precision :: a1,b1,c1,d1,w1 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 -! External variable - - double precision,external :: electron_number - - ! Defining enhancements factor for weight-dependent functionals 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 Ex = 0d0 - Exrr = 0d0 - ExrI = 0d0 - ExrrI = 0d0 - do iG=1,nGrid r = max(0d0,rhow(iG)) - rI = max(0d0,rho(iG)) if(r > threshold) then - e_p = Cx*r**(1d0/3d0) dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - Exrr = Exrr - weight(iG)*dedr*r*r - - if(rI > threshold) then - - ExrI = ExrI + weight(iG)*e_p*rI - ExrrI = ExrrI + weight(iG)*dedr*r*rI - - endif + Ex = Ex - weight(iG)*dedr*r*r endif 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 diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 index f62c173..1012f81 100644 --- a/src/eDFT/US51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/US51_lda_exchange_individual_energy.f90 @@ -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 @@ -10,16 +10,12 @@ subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,doNcentered integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) - double precision,intent(in) :: rho(nGrid) - logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! Local variables integer :: iG - double precision :: r,rI - double precision :: e,dedr - double precision :: Exrr,ExrI,ExrrI + double precision :: r + double precision :: dedr ! 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 - Exrr = 0d0 - ExrI = 0d0 - ExrrI = 0d0 + Ex = 0d0 do iG=1,nGrid r = max(0d0,rhow(iG)) - rI = max(0d0,rho(iG)) - if(r > threshold) then - e = CxLSDA*r**(1d0/3d0) dedr = 1d0/3d0*CxLSDA*r**(-2d0/3d0) - Exrr = Exrr - weight(iG)*dedr*r*r - - if(rI > threshold) then - - ExrI = ExrI + weight(iG)*e*rI - ExrrI = ExrrI + weight(iG)*dedr*r*rI - - endif + Ex = Ex - weight(iG)*dedr*r*r endif 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 diff --git a/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 index 5aff535..3ecd4f1 100644 --- a/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 @@ -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 @@ -11,9 +11,7 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) - double precision,intent(in) :: rho(nGrid,nspin) logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! Local variables @@ -54,187 +52,187 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent ! Initialization - Ec(:) = 0d0 - Ecrr(:) = 0d0 - EcrI(:) = 0d0 - EcrrI(:) = 0d0 +! Ec(:) = 0d0 +! Ecrr(:) = 0d0 +! EcrI(:) = 0d0 +! EcrrI(:) = 0d0 - do iG=1,nGrid +! do iG=1,nGrid - ra = max(0d0,rhow(iG,1)) - rb = max(0d0,rhow(iG,2)) +! ra = max(0d0,rhow(iG,1)) +! rb = max(0d0,rhow(iG,2)) - raI = max(0d0,rho(iG,1)) - rbI = max(0d0,rho(iG,2)) - +! raI = max(0d0,rho(iG,1)) +! rbI = max(0d0,rho(iG,2)) +! ! spin-up contribution - - r = ra - rI = raI +! +! r = ra +! rI = raI - if(r > threshold) then +! 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) +! 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 +! 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 ) ) +! 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 +! 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 - EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI - - end if +! EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI +! EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI +! +! end if - end if +! end if ! up-down contribution - - r = ra + rb - rI = raI + rbI +! +! r = ra + rb +! rI = raI + rbI - if(r > threshold) then +! if(r > threshold) then - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - z = (ra - rb)/r - x = sqrt(rs) - - fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0 - fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0)) - - d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0)) - - x_p = x*x + b_p*x + c_p - x_f = x*x + b_f*x + c_f - x_a = x*x + b_a*x + c_a - - xx0_p = x0_p*x0_p + b_p*x0_p + c_p - xx0_f = x0_f*x0_f + b_f*x0_f + c_f - xx0_a = x0_a*x0_a + b_a*x0_a + c_a - - q_p = sqrt(4d0*c_p - b_p*b_p) - q_f = sqrt(4d0*c_f - b_f*b_f) - q_a = sqrt(4d0*c_a - b_a*b_a) - - ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - - ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & - - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) ) - - ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & - - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) ) - - ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4 +! rs = (4d0*pi*r/3d0)**(-1d0/3d0) +! z = (ra - rb)/r +! x = sqrt(rs) +! +! fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0 +! fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0)) +! +! d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0)) +! +! x_p = x*x + b_p*x + c_p +! x_f = x*x + b_f*x + c_f +! x_a = x*x + b_a*x + c_a +! +! xx0_p = x0_p*x0_p + b_p*x0_p + c_p +! xx0_f = x0_f*x0_f + b_f*x0_f + c_f +! xx0_a = x0_a*x0_a + b_a*x0_a + c_a +! +! q_p = sqrt(4d0*c_p - b_p*b_p) +! q_f = sqrt(4d0*c_f - b_f*b_f) +! q_a = sqrt(4d0*c_a - b_a*b_a) +! +! ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & +! - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) +! +! ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & +! - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) ) +! +! ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & +! - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) ) +! +! ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4 - dzdr = (1d0 - z)/r - dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0)) - dfzdr = dzdr*dfzdz +! dzdr = (1d0 - z)/r +! dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0)) +! dfzdr = dzdr*dfzdz - drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) +! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) +! dxdrs = 0.5d0/sqrt(rs) - dxdx_p = 2d0*x + b_p - dxdx_f = 2d0*x + b_f - dxdx_a = 2d0*x + b_a +! dxdx_p = 2d0*x + b_p +! dxdx_f = 2d0*x + b_f +! dxdx_a = 2d0*x + b_a - decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) +! decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & +! - b_p*x0_p/xx0_p*( 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_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 ) ) +! 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 ) ) - decdr_p = drsdr*dxdrs*decdx_p - decdr_f = drsdr*dxdrs*decdx_f - decdr_a = drsdr*dxdrs*decdx_a +! decdr_p = drsdr*dxdrs*decdx_p +! decdr_f = drsdr*dxdrs*decdx_f +! decdr_a = drsdr*dxdrs*decdx_a - decdr = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdr/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdr*z**3 & - + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdr*z**4 + 4d0*(ec_f - ec_p)*fz*dzdr*z**3 +! 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 - 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 - EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI - - end if +! EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI +! EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI +! +! end if - end if +! end if ! spin-down contribution - - r = rb - rI = rbI - - if(r > threshold) then +! +! r = rb +! rI = rbI +! +! if(r > threshold) then - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - x = sqrt(rs) +! 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) +! x_f = x*x + b_f*x + c_f +! xx0_f = x0_f*x0_f + b_f*x0_f + c_f +! q_f = sqrt(4d0*c_f - b_f*b_f) - ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & - - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) ) +! ec_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) +! drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) +! 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 & - - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) ) +! decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f & +! - b_f*x0_f/xx0_f*( 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 - EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI +! EcrI(3) = EcrI(3) + weight(iG)*ec_f*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) - EcrI(2) = EcrI(2) - EcrI(1) - EcrI(3) - EcrrI(2) = EcrrI(2) - EcrrI(1) - EcrrI(3) +! 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 +! if(doNcentered) then - Ecrr(:) = kappa*Ecrr(:) - EcrI(:) = kappa*EcrI(:) +! Ecrr(:) = kappa*Ecrr(:) +! EcrI(:) = kappa*EcrI(:) - endif +! endif - Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:) +! Ec(:) = Ecrr(:) + EcrI(:) + EcrrI(:) end subroutine UVWN3_lda_correlation_individual_energy diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 index 1fb0e8b..1f327c9 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 @@ -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 @@ -11,9 +11,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) - double precision,intent(in) :: rho(nGrid,nspin) logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! 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_a,x0_a,xx0_a,b_a,c_a,x_a,q_a double precision :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a - double precision :: dzdra,dzdrb,dfzdra,dfzdrb,drsdr,decdr_p,decdr_f,decdr_a,decdra,decdrb,decdr + double precision :: dzdra,dzdrb,dfzdra,dfzdrb,drsdr,decdr_p,decdr_f,decdr_a,decdra,decdrb double precision :: ec_z,ec_p,ec_f,ec_a double precision :: fz,d2fz - double precision :: Ecrr(nsp) - double precision :: EcrI(nsp) - double precision :: EcrrI(nsp) - ! Output variables - double precision :: Ec(nsp) + double precision :: LZc(nspin) ! Parameters of the functional @@ -54,56 +48,14 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent ! Initialization - Ec(:) = 0d0 - Ecrr(:) = 0d0 - EcrI(:) = 0d0 - EcrrI(:) = 0d0 + LZc(:) = 0d0 do iG=1,nGrid ra = max(0d0,rhow(iG,1)) rb = max(0d0,rhow(iG,2)) - raI = max(0d0,rho(iG,1)) - rbI = max(0d0,rho(iG,2)) - r = ra + rb - rI = raI + rbI - -! 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 @@ -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)) - dzdra = + (1d0 - z)/r - dfzdra = dzdra*dfzdz - - dzdrb = - (1d0 + z)/r - dfzdrb = dzdrb*dfzdz - drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) 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_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 & - + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3 + if(ra > threshold) then - decdrb = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 & - + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3 + dzdra = + (1d0 - z)/r + dfzdra = dzdra*dfzdz + + decdra = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 & + + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3 + + 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 - 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 + LZc(2) = LZc(2) - weight(iG)*decdrb*rb*r 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 - 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 diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index 01f9217..aa4c46e 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -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 :: drho(:,:,:,:) - double precision :: E(nEns) - double precision :: Om(nEns) - integer :: ispin,iEns,iBas ! 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 -!!!!! -! 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 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 !------------------------------------------------------------------------ - 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, & - Cx_choice,doNcentered) + call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & + AO,dAO,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,occnum,Cx_choice,doNcentered) end subroutine eDFT_UKS diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index bd8e7d8..5ce1641 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -34,7 +34,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' ENSEMBLE ENERGIES' 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(*,*) @@ -55,66 +55,66 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, ! Kinetic energy !------------------------------------------------------------------------ - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES' - write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au' - end do - write(*,'(A60)') '-------------------------------------------------' - write(*,*) +! write(*,'(A60)') '-------------------------------------------------' +! write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES' +! write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' +! write(*,*) !------------------------------------------------------------------------ ! Potential energy !------------------------------------------------------------------------ - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES' - write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au' - end do - write(*,'(A60)') '-------------------------------------------------' - write(*,*) +! write(*,'(A60)') '-------------------------------------------------' +! write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES' +! write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' +! write(*,*) !------------------------------------------------------------------------ ! Hartree energy !------------------------------------------------------------------------ - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES' - write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EJ(:,iEns)),' au' - end do - write(*,'(A60)') '-------------------------------------------------' - write(*,*) +! write(*,'(A60)') '-------------------------------------------------' +! write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES' +! write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EJ(:,iEns)),' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' +! write(*,*) !------------------------------------------------------------------------ ! Exchange energy !------------------------------------------------------------------------ - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES' - write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au' - end do - write(*,'(A60)') '-------------------------------------------------' - write(*,*) +! write(*,'(A60)') '-------------------------------------------------' +! write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES' +! write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' +! write(*,*) !------------------------------------------------------------------------ ! Correlation energy !------------------------------------------------------------------------ - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES' - write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au' - end do - write(*,'(A60)') '-------------------------------------------------' - write(*,*) +! write(*,'(A60)') '-------------------------------------------------' +! write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES' +! write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' +! write(*,*) !------------------------------------------------------------------------ ! Auxiliary energies @@ -179,18 +179,18 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL ENERGIES ' write(*,'(A60)') '-------------------------------------------------' - do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' - end do - write(*,'(A60)') '-------------------------------------------------' +! do iEns=1,nEns +! write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' +! end do +! write(*,'(A60)') '-------------------------------------------------' do iEns=2,nEns write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au' write(*,*) - write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(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(*,*) +! 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)') ' xc energy contribution : ',Omxc(iEns), ' au' +! write(*,*) write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',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(*,*) - 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)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV' - write(*,*) +! 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)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV' +! write(*,*) write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' diff --git a/src/eDFT/unrestricted_correlation_individual_energy.f90 b/src/eDFT/unrestricted_correlation_individual_energy.f90 index 4bcc8af..0265d4d 100644 --- a/src/eDFT/unrestricted_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_correlation_individual_energy.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho, & - doNcentered,kappa,Ec) +subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc) ! 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) :: rhow(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 - double precision,intent(in) :: kappa ! Local variables - double precision :: EcLDA(nsp) - double precision :: EcGGA(nsp) + double precision :: LZcLDA(nspin) + double precision :: LZcGGA(nspin) ! Output variables - double precision,intent(out) :: Ec(nsp) + double precision,intent(out) :: LZc(nspin) select case (rung) @@ -37,14 +33,13 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns case(0) - Ec(:) = 0d0 + LZc(:) = 0d0 ! LDA functionals case(1) - call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho, & - doNcentered,kappa,Ec) + call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,doNcentered,LZc) ! GGA functionals @@ -62,7 +57,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns case(4) - call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) + call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,LZc) end select diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90 index 721425e..3144892 100644 --- a/src/eDFT/unrestricted_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_exchange_individual_energy.f90 @@ -1,5 +1,5 @@ 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 @@ -20,14 +20,10 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE integer,intent(in) :: nBas double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: Pw(nBas,nBas) - double precision,intent(in) :: P(nBas,nBas) double precision,intent(in) :: rhow(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 logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! Output variables @@ -46,25 +42,25 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE case(1) call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,& - rhow,rho,Cx_choice,doNcentered,kappa,Ex) + rhow,Cx_choice,doNcentered,Ex) ! GGA functionals 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 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 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 diff --git a/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 b/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 index bec1c8a..f60dfdb 100644 --- a/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 @@ -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 @@ -8,7 +8,6 @@ subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) integer,intent(in) :: nBas double precision,intent(in) :: Pw(nBas,nBas) - double precision,intent(in) :: P(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! 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) - Ex = trace_matrix(nBas,matmul(P ,Fx)) & - - 0.5d0*trace_matrix(nBas,matmul(Pw,Fx)) + Ex = - 0.5d0*trace_matrix(nBas,matmul(Pw,Fx)) end subroutine unrestricted_fock_exchange_individual_energy diff --git a/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 b/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 index 25b3834..ba9645b 100644 --- a/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 @@ -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 @@ -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) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) - double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(ncart,nGrid) ! Output variables diff --git a/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 b/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 index 13de347..76ab137 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 @@ -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 @@ -14,8 +14,6 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) - double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(ncart,nGrid) ! Output variables diff --git a/src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 b/src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 index d542993..eb95447 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 @@ -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 @@ -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) :: rhow(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 double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: Pw(nBas,nBas) - double precision,intent(in) :: P(nBas,nBas) ! Output variables @@ -32,7 +29,7 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we 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 diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index aeae00a..34b867a 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -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, & - 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) ! 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) :: FxHF(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) integer,intent(in) :: Cx_choice 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 :: Ec(nsp,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 :: 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 :: kappa(:) - + double precision :: E(nEns) + double precision :: Om(nEns) + 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)) nEl(:) = 0d0 @@ -95,83 +97,91 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Kinetic energy !------------------------------------------------------------------------ - do ispin=1,nspin - do iEns=1,nEns - ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) - end do - end do +! do ispin=1,nspin +! do iEns=1,nEns +! ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) +! end do +! end do !------------------------------------------------------------------------ ! Potential energy !------------------------------------------------------------------------ - do iEns=1,nEns - do ispin=1,nspin - EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) - end do - end do +! do iEns=1,nEns +! do ispin=1,nspin +! EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) +! end do +! end do !------------------------------------------------------------------------ ! Individual Hartree energy !------------------------------------------------------------------------ - do iEns=1,nEns +! do iEns=1,nEns - do ispin=1,nspin - call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) - end do +! do ispin=1,nspin +! call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) +! end do - if(doNcentered) then - - 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))) - - 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))) & - - 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))) - - 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))) - else +! if(doNcentered) then +! +! 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))) +! +! 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))) & +! - 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))) +! +! 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))) +! else - EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) +! 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))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) +! 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))) - EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & - - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - end if +! 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 + +!------------------------------------------------------------------------ +! Individual Hartree energy +!------------------------------------------------------------------------ + + do ispin=1,nspin + LZH(ispin) = -0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1)+Pw(:,:,2),J(:,:,ispin))) end do !------------------------------------------------------------------------ ! Individual exchange energy !------------------------------------------------------------------------ - do iEns=1,nEns - do ispin=1,nspin - call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, & - Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & - rho(:,ispin,iEns),drho(:,:,ispin,iEns),Cx_choice,doNcentered,kappa(iEns), & - Ex(ispin,iEns)) - - end do + do ispin=1,nspin + call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, & + Pw(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered, & + LZx(ispin)) end do !------------------------------------------------------------------------ ! Individual correlation energy !------------------------------------------------------------------------ - do iEns=1,nEns - 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 + call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,doNcentered,LZc) + +!------------------------------------------------------------------------ +! Individual exchange-correlation energy +!------------------------------------------------------------------------ + + print*,LZc + LZHxc(:) = LZH(:) + LZx(:) + LZc(:) !------------------------------------------------------------------------ ! Compute auxiliary energies @@ -199,10 +209,14 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! 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 - 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)) + E(iEns) = sum(Eaux(:,iEns)) + sum(LZHxc(:)) + sum(ExcDD(:,iEns)) end do !------------------------------------------------------------------------ @@ -212,9 +226,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns Om(iEns) = E(iEns) - E(1) - Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1)) - Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1)) - Omxc(iEns) = Exc(iEns) - Exc(1) +! Omx(iEns) = sum(Ex(:,iEns)) - sum(Ex(:,1)) +! Omc(iEns) = sum(Ec(:,iEns)) - sum(Ec(:,1)) +! Omxc(iEns) = Exc(iEns) - Exc(1) Omaux(iEns) = sum(Eaux(:,iEns)) - sum(Eaux(:,1)) diff --git a/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 b/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 index d654414..eef2395 100644 --- a/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 @@ -1,5 +1,4 @@ -subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho, & - doNcentered,kappa,Ec) +subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,doNcentered,LZc) ! 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 double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid,nspin) - double precision,intent(in) :: rho(nGrid,nspin) logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! Output variables - double precision :: Ec(nsp) + double precision :: LZc(nspin) ! Select correlation functional @@ -37,11 +34,11 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns, 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) - 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 diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 index ca75792..2564086 100644 --- a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,& - Cx_choice,doNcentered,kappa,Ex) +subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& + Cx_choice,doNcentered,Ex) ! 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 double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rhow(nGrid) - double precision,intent(in) :: rho(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered - double precision,intent(in) :: kappa ! Output variables @@ -33,12 +31,12 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn 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) - call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,& - Cx_choice,doNcentered,kappa,Ex) + call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow, & + Cx_choice,doNcentered,Ex) case default diff --git a/src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 b/src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 index 5f49807..b34154f 100644 --- a/src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 +++ b/src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 @@ -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 @@ -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) :: rhow(nGrid) double precision,intent(in) :: drhow(ncart,nGrid) - double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(ncart,nGrid) ! Output variables