LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-30 15:15:38 +01:00
Code Releases Activity

removed a lot of floating point exceptions in DFT

Browse Source
This commit is contained in:
Emmanuel Giner 2020-06-08 12:52:08 +02:00
parent ff15a50895
commit 3987b9794d
1 changed files with 27 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
src/dft_utils_func/on_top_from_ueg.irp.f
View File

@ -35,7 +35,11 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
if (dabs(rho) > 1.d-12) then if (dabs(rho) > 1.d-12) then
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
x = -d2*rs x = -d2*rs
g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*exp(x) if(dabs(x).lt.50.d0)then
g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*dexp(x)
else
g0_UEG_mu_inf= 0.d0
endif
else else
g0_UEG_mu_inf= 0.d0 g0_UEG_mu_inf= 0.d0
endif endif
@ -67,7 +71,11 @@ double precision function g0_UEG_mu(mu,rho_a,rho_b)
kf = (alpha*rs)**(-1d0) kf = (alpha*rs)**(-1d0)
zeta = mu / kf zeta = mu / kf
x = -d2*rs*h_func(zeta)/ahd x = -d2*rs*h_func(zeta)/ahd
g0_UEG_mu = (exp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) ) if(dabs(x).lt.50.d0)then
g0_UEG_mu = (dexp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
else
g0_UEG_mu = 0.d0
endif
end end
@ -81,11 +89,11 @@ double precision function h_func(zeta)
pi = 4d0 * datan(1d0) pi = 4d0 * datan(1d0)
ahd = -0.36583d0 ahd = -0.36583d0
alpha = (4d0/(9d0*pi))**(1d0/3d0) alpha = (4d0/(9d0*pi))**(1d0/3d0)
a1 = -(6d0*alpha/pi)*(1d0-log(2d0)) a1 = -(6d0*alpha/pi)*(1d0-dlog(2d0))
b1 = 1.4919d0 b1 = 1.4919d0
b3 = 1.91528d0 b3 = 1.91528d0
a2 = ahd * b3 a2 = ahd * b3
b2 = (a1 - (b3*alpha/sqrt(pi)))/ahd b2 = (a1 - (b3*alpha/dsqrt(pi)))/ahd
h_func = (a1*zeta**2d0 + a2*zeta**3d0) / (1d0 + b1*zeta + b2*zeta**2d0 + b3*zeta**3d0) h_func = (a1*zeta**2d0 + a2*zeta**3d0) / (1d0 + b1*zeta + b2*zeta**2d0 + b3*zeta**3d0)
end end
@ -111,11 +119,23 @@ end
D1 = -0.0127713d0 D1 = -0.0127713d0
E1 = 0.00185898d0 E1 = 0.00185898d0
B1 = 0.7317d0 - F1 B1 = 0.7317d0 - F1
rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) if(dabs(rho).gt.1.d-20)then
rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0)
else
rs = (3.d0 / (4.d0*pi*1.d-20))**(1.d0/3.d0)
endif
g0 = g0_UEG_mu_inf(rho_a, rho_b) g0 = g0_UEG_mu_inf(rho_a, rho_b)
dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*exp(-F1*rs) if(dabs(F1*rs).lt.50.d0)then
dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*dexp(-F1*rs)
else
dg0drs = 0.d0
endif
if(dabs(rho).gt.1.d-20)then
dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
else
dg0drho = -((6.d0*dsqrt(pi)*1.d-40)**(-2.d0/3.d0))*dg0drs
endif
end subroutine g0_dg0 end subroutine g0_dg0