LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

fix problem with threshold

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-11-02 11:07:34 +01:00
parent f95e0d4dfc
commit 52c224416f
3 changed files with 57 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
input/dft
View File

@ -19,21 +19,21 @@
# Number of states in ensemble (nEns) # Number of states in ensemble (nEns)
4 4
# occupation numbers # occupation numbers
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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 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 0 0 0 0 0 0 0 0 0 0 0 0 0 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 0 0 0 0 0 0 0 0 0 0 0 0 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) # Ensemble weights: wEns(1),...,wEns(nEns-1)
0.25 0.25 0.25 1.00 0.00 0.00
# N-centered? # N-centered?
T F
# Parameters for CC weight-dependent exchange functional # Parameters for CC weight-dependent exchange functional
0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0

21
src/eDFT/UVWN3_lda_correlation_individual_energy.f90
View File

@ -72,7 +72,7 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
r = ra r = ra
rI = raI rI = raI
if(r > threshold .or. rI > threshold) then if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs) x = sqrt(rs)
@ -95,17 +95,22 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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
EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI
EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI
end if end if
end if
! up-down contribution ! up-down contribution
r = ra + rb r = ra + rb
rI = raI + rbI rI = raI + rbI
if(r > threshold .or. rI > threshold) then if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r z = (ra - rb)/r
@ -167,17 +172,22 @@ 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 + (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
EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
end if end if
end if
! spin-down contribution ! spin-down contribution
r = rb r = rb
rI = rbI rI = rbI
if(r > threshold .or. rI > threshold) then if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs) x = sqrt(rs)
@ -200,11 +210,16 @@ subroutine UVWN3_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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
EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
end if end if
end if
end do end do
Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3) Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3)

21
src/eDFT/UVWN5_lda_correlation_individual_energy.f90
View File

@ -72,7 +72,7 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
r = ra r = ra
rI = raI rI = raI
if(r > threshold .or. rI > threshold) then if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs) x = sqrt(rs)
@ -95,17 +95,22 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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
EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI EcrI(1) = EcrI(1) + weight(iG)*ec_f*rI
EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI EcrrI(1) = EcrrI(1) + weight(iG)*decdr_f*r*rI
end if end if
end if
! up-down contribution ! up-down contribution
r = ra + rb r = ra + rb
rI = raI + rbI rI = raI + rbI
if(r > threshold .or. rI > threshold) then if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) rs = (4d0*pi*r/3d0)**(-1d0/3d0)
z = (ra - rb)/r z = (ra - rb)/r
@ -167,17 +172,22 @@ 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 + (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
EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI EcrI(2) = EcrI(2) + weight(iG)*ec_z*rI
EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI EcrrI(2) = EcrrI(2) + weight(iG)*decdr*r*rI
end if end if
end if
! spin-down contribution ! spin-down contribution
r = rb r = rb
rI = rbI rI = rbI
if(r > threshold .or. rI > threshold) then if(r > threshold) then
rs = (4d0*pi*r/3d0)**(-1d0/3d0) rs = (4d0*pi*r/3d0)**(-1d0/3d0)
x = sqrt(rs) x = sqrt(rs)
@ -200,11 +210,16 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent
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
EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI EcrI(3) = EcrI(3) + weight(iG)*ec_f*rI
EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI EcrrI(3) = EcrrI(3) + weight(iG)*decdr_f*r*rI
end if end if
end if
end do end do
Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3) Ecrr(2) = Ecrr(2) - Ecrr(1) - Ecrr(3)