From 52c224416f3bf668db1e9538549189fe4dc4f0db Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 2 Nov 2021 11:07:34 +0100 Subject: [PATCH] fix problem with threshold --- input/dft | 18 +++++----- ...VWN3_lda_correlation_individual_energy.f90 | 33 ++++++++++++++----- ...VWN5_lda_correlation_individual_energy.f90 | 33 ++++++++++++++----- 3 files changed, 57 insertions(+), 27 deletions(-) diff --git a/input/dft b/input/dft index ecbc4ae..415fecc 100644 --- a/input/dft +++ b/input/dft @@ -19,21 +19,21 @@ # Number of states in ensemble (nEns) 4 # 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 0 0 0 0 0 0 0 0 0 0 0 0 0 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 0 0 0 0 0 0 0 0 0 0 0 0 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.25 0.25 0.25 + 1.00 0.00 0.00 # N-centered? - T + F # Parameters for CC weight-dependent exchange functional 0.0 0.0 0.0 0.0 0.0 0.0 diff --git a/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 index dcdbc00..5aff535 100644 --- a/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 @@ -72,7 +72,7 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent r = ra rI = raI - if(r > threshold .or. rI > threshold) then + if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) x = sqrt(rs) @@ -95,9 +95,14 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent decdr_f = drsdr*dxdrs*decdx_f 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 + + 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 @@ -105,7 +110,7 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent r = ra + rb rI = raI + rbI - if(r > threshold .or. rI > threshold) then + if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) z = (ra - rb)/r @@ -167,9 +172,14 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent + (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 - EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI - EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI + + if(rI > threshold) then + + EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI + EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI + end if + end if ! spin-down contribution @@ -177,7 +187,7 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent r = rb rI = rbI - if(r > threshold .or. rI > threshold) then + if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) x = sqrt(rs) @@ -200,8 +210,13 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent decdr_f = drsdr*dxdrs*decdx_f 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 + + 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 diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 index 3861750..d6407ac 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 @@ -72,7 +72,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent r = ra rI = raI - if(r > threshold .or. rI > threshold) then + if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) x = sqrt(rs) @@ -95,8 +95,13 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent decdr_f = drsdr*dxdrs*decdx_f 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 + + 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 @@ -105,7 +110,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent r = ra + rb rI = raI + rbI - if(r > threshold .or. rI > threshold) then + if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) z = (ra - rb)/r @@ -167,8 +172,13 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent + (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 - EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI - EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI + + if(rI > threshold) then + + EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI + EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI + + end if end if @@ -177,7 +187,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent r = rb rI = rbI - if(r > threshold .or. rI > threshold) then + if(r > threshold) then rs = (4d0*pi*r/3d0)**(-1d0/3d0) x = sqrt(rs) @@ -200,8 +210,13 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent decdr_f = drsdr*dxdrs*decdx_f 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 + + 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