From 48b0952b557b69a6a55040a98d0f7d1633f12229 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 6 Jun 2020 18:07:26 +0200 Subject: [PATCH] removed floating points exceptions in DFT --- src/ao_basis/aos_in_r.irp.f | 2 + src/dft_utils_func/on_top_from_ueg.irp.f | 30 +++-- src/dft_utils_func/routines_exc_sr_lda.irp.f | 123 ++++++++++++++----- 3 files changed, 117 insertions(+), 38 deletions(-) diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 29e52169..ecba6027 100644 --- a/src/ao_basis/aos_in_r.irp.f +++ b/src/ao_basis/aos_in_r.irp.f @@ -97,6 +97,7 @@ subroutine give_all_aos_at_r(r,aos_array) dz2 = dz**power_ao(3) do l = 1,ao_prim_num(k) beta = ao_expo_ordered_transp_per_nucl(l,j,i) + if(dabs(beta*r2).gt.40.d0)cycle aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) enddo aos_array(k) = aos_array(k) * dx2 * dy2 * dz2 @@ -162,6 +163,7 @@ subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) accu_2 = 0.d0 do l = 1,ao_prim_num(k) beta = ao_expo_ordered_transp_per_nucl(l,j,i) + if(beta*r2.gt.40.d0)cycle contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) accu_1 += contrib accu_2 += contrib * beta diff --git a/src/dft_utils_func/on_top_from_ueg.irp.f b/src/dft_utils_func/on_top_from_ueg.irp.f index 70560a7a..8912d954 100644 --- a/src/dft_utils_func/on_top_from_ueg.irp.f +++ b/src/dft_utils_func/on_top_from_ueg.irp.f @@ -32,10 +32,10 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b) C = 0.08193d0 D = -0.01277d0 E = 0.001859d0 - if (dabs(rho) > 1.d-12) then + x = -d2*rs + if (dabs(rho) > 1.d-12.and.dabs(x).lt.20.d0) then rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 - x = -d2*rs - g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*exp(x) + 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 @@ -67,7 +67,11 @@ double precision function g0_UEG_mu(mu,rho_a,rho_b) kf = (alpha*rs)**(-1d0) zeta = mu / kf 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.40.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 @@ -111,11 +115,23 @@ end D1 = -0.0127713d0 E1 = 0.00185898d0 B1 = 0.7317d0 - F1 - rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) + if(dabs(rho).gt.1.d-12)then + rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) + else + rs = (3.d0 / (4.d0*pi*1.d-12))**(1.d0/3.d0) + endif 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) - dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs + if(dabs(F1*rs).lt.40.d0)then + 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(rho**2.gt.1.d-12)then + dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs + else + dg0drho = -((6.d0*dsqrt(pi)*1.d-12)**(-2.d0/3.d0))*dg0drs + endif end subroutine g0_dg0 diff --git a/src/dft_utils_func/routines_exc_sr_lda.irp.f b/src/dft_utils_func/routines_exc_sr_lda.irp.f index ea1dcd69..676ac5d5 100644 --- a/src/dft_utils_func/routines_exc_sr_lda.irp.f +++ b/src/dft_utils_func/routines_exc_sr_lda.irp.f @@ -153,8 +153,13 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b) !Intermediate values of a elseif (a.le.100d0) then - ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3)) - vx_a = -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(dexp(-f14/a2)-z1)+mu/sqpi * derf(f12/a) + if(dabs(f14/a2).lt.40.d0)then + ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3)) + vx_a = -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(dexp(-f14/a2)-z1)+mu/sqpi * derf(f12/a) + else + ex_a = 0.d0 + vx_a = 0.d0 + endif !Expansion for large a @@ -185,8 +190,13 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b) !Intermediate values of a elseif (a.le.100d0) then - ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3)) - vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(dexp(-f14/a2)-z1)+mu/sqpi* derf(f12/a) + if(dabs(f14/a2).lt.40.d0)then + ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3)) + vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(dexp(-f14/a2)-z1)+mu/sqpi* derf(f12/a) + else + ex_b = 0.d0 + vx_b = 0.d0 + endif !Expansion for large a elseif (a.lt.1.d+9) then @@ -254,7 +264,11 @@ end double precision derf eta=19.0d0 - fak=2.540118935556d0*dexp(-eta*a*a) + if(dabs(eta*a*a).lt.40.d0)then + fak=2.540118935556d0*dexp(-eta*a*a) + else + fak = 0.d0 + endif if(a .lt. 0.075d0) then ! expansion for small mu to avoid numerical problems @@ -301,7 +315,11 @@ end double precision t1,t2,tdexp,t3,t4,t5 eta=19.0d0 - fak=2.540118935556d0*dexp(-eta*a*a) + if(dabs(eta*a*a).lt.40.d0)then + fak=2.540118935556d0*dexp(-eta*a*a) + else + fak=0.d0 + endif dfakda=-2.0d0*eta*a*fak if(a .lt. 0.075d0) then @@ -372,9 +390,16 @@ subroutine ecorrlr(rs,z,mu,eclr) t3a = 0.31d0 b0=adib*rs - - d2anti=(q1a*rs+q2a*rs**2)*exp(-abs(q3a)*rs)/rs**2 - d3anti=(t1a*rs+t2a*rs**2)*exp(-abs(t3a)*rs)/rs**3 + if(dabs(q3a*rs).lt.40.d0)then + d2anti=(q1a*rs+q2a*rs**2)*dexp(-dabs(q3a)*rs)/rs**2 + else + d2anti=0.d0 + endif + if(dabs(t3a*rs).lt.40.d0)then + d3anti=(t1a*rs+t2a*rs**2)*dexp(-dabs(t3a)*rs)/rs**3 + else + d3anti=0.d0 + endif coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0f(rs)-0.5d0) @@ -468,20 +493,44 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) b0=adib*rs - d2anti=(q1a+q2a*rs)*exp(-q3a*rs)/rs - d3anti=(t1a+t2a*rs)*exp(-t3a*rs)/rs**2 + if(dabs(q3a*rs).lt.40.d0)then + d2anti=(q1a+q2a*rs)*dexp(-q3a*rs)/rs + else + d2anti=0.d0 + endif + if(dabs(t3a*rs).lt.40.d0)then + d3anti=(t1a+t2a*rs)*dexp(-t3a*rs)/rs**2 + else + d3anti=0.d0 + endif - d2antid=-((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*exp(-q3a*rs) - d3antid=-((rs*t2a*(1d0 + rs*t3a) + t1a*(2d0 + rs*t3a))/rs**3)*exp(-rs*t3a) + if(dabs(q3a*rs).lt.40.d0)then + d2antid=-((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*dexp(-q3a*rs) + else + d2antid=0.d0 + endif + if(dabs(t3a*rs).lt.40.d0)then + d3antid=-((rs*t2a*(1d0 + rs*t3a) + t1a*(2d0 + rs*t3a))/rs**3)*dexp(-rs*t3a) + else + d3antid=0.d0 + endif !SCD - d2antidd = exp(-q3a*rs)/rs**3*( & - q3a**2*q1a*rs**2+q2a*q3a**2*rs**3 & - +2.d0*q3a*q1a*rs+2.d0*q1a) - d3antidd = exp(-t3a*rs)/rs**4* & + if(dabs(q3a*rs).lt.40.d0)then + d2antidd = dexp(-q3a*rs)/rs**3*( & + q3a**2*q1a*rs**2+q2a*q3a**2*rs**3 & + +2.d0*q3a*q1a*rs+2.d0*q1a) + else + d2antidd = 0.d0 + endif + if(dabs(t3a*rs).lt.40.d0)then + d3antidd = dexp(-t3a*rs)/rs**4* & (2.d0*t3a*t2a*rs**2 + 2.d0*t2a*rs & + t1a*t3a**2*rs**2 + t2a*t3a**2*rs**3 & + 4.d0*t1a*t3a*rs + 6.d0*t1a) + else + d3antidd=0.d0 + endif !SCF coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0f(rs)-0.5d0) coe2rs=-3.d0/8.d0/rs**3*(1.d0-z**2)*g0d(rs)+ & @@ -759,8 +808,12 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) D0f = 0.752411d0 E0f = -0.0127713d0 F0f = 0.00185898d0 - g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+ & - F0f*x**4)*exp(-abs(D0f)*x)/2.d0 + if(dabs(D0f*x).lt.40.d0)then + g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+ & + F0f*x**4)*dexp(-dabs(D0f)*x)/2.d0 + else + g0f = 0.d0 + endif return end @@ -774,11 +827,15 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) Dg0 = -0.0127713d0 Eg0 = 0.00185898d0 Bg0 =0.7317d0-Fg0 - expsum=exp(-Fg0*rs) - g0d=(-Bg0+2d0*Cg0*rs+3d0*Dg0*rs**2+4d0*Eg0*rs**3)/2.d0 & - *expsum & - - (Fg0*(1d0 - Bg0*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4))/ & - 2.d0*expsum + if(dabs(Fg0*rs).lt.40.d0)then + expsum=dexp(-Fg0*rs) + g0d=(-Bg0+2d0*Cg0*rs+3d0*Dg0*rs**2+4d0*Eg0*rs**3)/2.d0 & + *expsum & + - (Fg0*(1d0 - Bg0*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4))/ & + 2.d0*expsum + else + g0d = 0.d0 + endif return end !SCD @@ -791,13 +848,17 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) Dg0 = -0.0127713d0 Eg0 = 0.00185898d0 Bg0 = 0.7317d0-Fg0 - expsum=exp(-Fg0*rs) - g0dd = (2.d0*Cg0+6.d0*Dg0*rs+12.d0*Eg0*rs**2)/2.d0* & - expsum & - - (-Bg0+2.d0*Cg0*rs+3.d0*Dg0*rs**2+4.d0*Eg0*rs**3)*Fg0* & - expsum & - + (1.d0-Bg0*rs+Cg0*rs**2+Dg0*rs**3+Eg0*rs**4)*Fg0**2* & - expsum/(2.d0) + if(dabs(Fg0*rs).lt.40.d0)then + expsum=dexp(-Fg0*rs) + g0dd = (2.d0*Cg0+6.d0*Dg0*rs+12.d0*Eg0*rs**2)/2.d0* & + expsum & + - (-Bg0+2.d0*Cg0*rs+3.d0*Dg0*rs**2+4.d0*Eg0*rs**3)*Fg0* & + expsum & + + (1.d0-Bg0*rs+Cg0*rs**2+Dg0*rs**3+Eg0*rs**4)*Fg0**2* & + expsum/(2.d0) + else + g0dd = 0.d0 + endif return end !SCF