From 346ea702c44a1a2fd99265c504830df84bc7bf2f Mon Sep 17 00:00:00 2001 From: Emmanuel Giner LCT Date: Tue, 19 May 2020 11:46:29 +0200 Subject: [PATCH 01/11] minor modifs on basis correction --- src/basis_correction/basis_correction.irp.f | 1 + 1 file changed, 1 insertion(+) diff --git a/src/basis_correction/basis_correction.irp.f b/src/basis_correction/basis_correction.irp.f index 11e53fcd..8d3264bc 100644 --- a/src/basis_correction/basis_correction.irp.f +++ b/src/basis_correction/basis_correction.irp.f @@ -7,6 +7,7 @@ program basis_correction touch read_wf no_core_density = .True. touch no_core_density + provide ao_two_e_integrals_in_map provide mo_two_e_integrals_in_map call print_basis_correction ! call print_e_b From 48b0952b557b69a6a55040a98d0f7d1633f12229 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 6 Jun 2020 18:07:26 +0200 Subject: [PATCH 02/11] 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 From 47528188b595e99abe597d82baf5ac946cfdfdb1 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 11:50:41 +0200 Subject: [PATCH 03/11] fixed floating point exception in AOs --- src/ao_basis/aos_in_r.irp.f | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 29e52169..7d400222 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 From 764801033157abdf459b96a0aa8cbca800a125ae Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 11:55:44 +0200 Subject: [PATCH 04/11] fixed another floating point exception in aos_in_r.irp.f --- src/ao_basis/aos_in_r.irp.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 7d400222..7fcb980a 100644 --- a/src/ao_basis/aos_in_r.irp.f +++ b/src/ao_basis/aos_in_r.irp.f @@ -163,6 +163,8 @@ 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) + contrib = 0.d0 + if(beta*r2.gt.50.d0)cycle contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) accu_1 += contrib accu_2 += contrib * beta From 408257dbfd594e4534b6b70801eaeef78fc15442 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 12:00:25 +0200 Subject: [PATCH 05/11] removed a floating point exception in routines_exc_sr_lda.irp.f --- src/dft_utils_func/routines_exc_sr_lda.irp.f | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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..129185ed 100644 --- a/src/dft_utils_func/routines_exc_sr_lda.irp.f +++ b/src/dft_utils_func/routines_exc_sr_lda.irp.f @@ -254,7 +254,11 @@ end double precision derf eta=19.0d0 - fak=2.540118935556d0*dexp(-eta*a*a) + if(dabs(eta*a*a).lt.50.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 From ff15a508954dd2a2078bab77d91e4ec33777294f Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 12:24:35 +0200 Subject: [PATCH 06/11] removed a floating point exception in routines_exc_sr_lda.irp.f --- src/dft_utils_func/routines_exc_sr_lda.irp.f | 121 +++++++++++++------ 1 file changed, 82 insertions(+), 39 deletions(-) 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 129185ed..21fde301 100644 --- a/src/dft_utils_func/routines_exc_sr_lda.irp.f +++ b/src/dft_utils_func/routines_exc_sr_lda.irp.f @@ -114,6 +114,7 @@ subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b) double precision :: f12,f13,f14,f32,f23,f43,f16 double precision :: ckf double precision :: a, akf,a2, a3 + double precision :: exp_f14a2 z0 = 0.D0 z1 = 1.D0 @@ -153,8 +154,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.50.d0)then + exp_f14a2 = dexp(-f14/a2) + else + exp_f14a2 = 0.d0 + endif + ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)* exp_f14a2 -z3*a+z4*a3)) + vx_a = -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(exp_f14a2 - z1)+mu/sqpi * derf(f12/a) !Expansion for large a @@ -185,8 +191,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.50.d0)then + exp_f14a2 = dexp(-f14/a2) + else + exp_f14a2 = 0.d0 + endif + ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*exp_f14a2-z3*a+z4*a3)) + vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(exp_f14a2-z1)+mu/sqpi* derf(f12/a) !Expansion for large a elseif (a.lt.1.d+9) then @@ -305,7 +316,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.50.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 @@ -377,17 +392,29 @@ subroutine ecorrlr(rs,z,mu,eclr) 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 + double precision :: exp_q3a_rs + if(dabs(q3a*rs).lt.50.d0)then + exp_q3a_rs = dexp(-dabs(q3a)*rs) + else + exp_q3a_rs = 0.d0 + endif + d2anti=(q1a*rs+q2a*rs**2)*exp_q3a_rs/rs**2 + double precision :: exp_t3a_rs + if(dabs(t3a*rs).lt.50.d0)then + exp_t3a_rs = dexp(-dabs(t3a)*rs) + else + exp_t3a_rs = 0.d0 + endif + d3anti=(t1a*rs+t2a*rs**2)*exp_t3a_rs/rs**3 coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0f(rs)-0.5d0) - coe3=-(1.d0-z**2)*g0f(rs)/(sqrt(2.d0*pi)*rs**3) + coe3=-(1.d0-z**2)*g0f(rs)/(dsqrt(2.d0*pi)*rs**3) - if(abs(z).eq.1.d0) then + if(dabs(z).eq.1.d0) then coe4=-9.d0/64.d0/rs**3*(dpol(rs) -cf**2*2d0**(5.d0/3.d0)/5.d0/rs**2) - coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*dpol(rs) + coe5=-9.d0/40.d0/(dsqrt(2.d0*pi)*rs**3)*dpol(rs) else @@ -397,7 +424,7 @@ subroutine ecorrlr(rs,z,mu,eclr) (1.-z**2)*d2anti-cf**2/10.d0*((1.d0+z)**(8.d0/3.d0) & +(1.-z)**(8.d0/3.d0))/rs**2) - coe5=-9.d0/40.d0/(sqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2 & + coe5=-9.d0/40.d0/(dsqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2 & *dpol(rs*(2.d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2 & *dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)* & d3anti) @@ -413,13 +440,13 @@ subroutine ecorrlr(rs,z,mu,eclr) a3=b0**8*coe3 a4=b0**6*(b0**2*coe2+4.d0*ec) - if(mu*sqrt(rs)/phi.lt.0.d0)then + if(mu*dsqrt(rs)/phi.lt.0.d0)then print*,'phi',phi print*,'mu ',mu print*,'rs ',rs stop -1 endif - eclr=(phi**3*Qrpa(mu*sqrt(rs)/phi)+a1*mu**3+a2*mu**4+a3*mu**5+ & + eclr=(phi**3*Qrpa(mu*dsqrt(rs)/phi)+a1*mu**3+a2*mu**4+a3*mu**5+ & a4*mu**6+b0**8*mu**8*ec)/((1.d0+b0**2*mu**2)**4) return @@ -471,18 +498,29 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) !SCF b0=adib*rs + double precision :: exp_q3a_rs,exp_t3a_rs + if(dabs(q3a*rs).lt.50.d0)then + exp_q3a_rs = dexp(-q3a*rs) + else + exp_q3a_rs = 0.d0 + endif + if(dabs(t3a*rs).lt.50.d0)then + exp_t3a_rs = dexp(-t3a*rs) + else + exp_t3a_rs = 0.d0 + endif - d2anti=(q1a+q2a*rs)*exp(-q3a*rs)/rs - d3anti=(t1a+t2a*rs)*exp(-t3a*rs)/rs**2 + d2anti=(q1a+q2a*rs)*exp_q3a_rs/rs + d3anti=(t1a+t2a*rs)*exp_t3a_rs/rs**2 - 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) + 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_t3a_rs !SCD - d2antidd = exp(-q3a*rs)/rs**3*( & + 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* & + d3antidd = exp_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) @@ -530,7 +568,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) +dpoldd(rs)*rs**4) coe4zd = 0.d0 - coe5rsd = -9.d0/40.d0/sqrt(2.d0/pi)/rs**5* & + coe5rsd = -9.d0/40.d0/dsqrt(2.d0/pi)/rs**5* & (12.d0*dpol(rs)-6.d0*rs*dpold(rs) & +rs**2*dpoldd(rs)) coe5zd = 0.d0 @@ -674,7 +712,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) a5zd= 0.d0 !SCF - x=mu*sqrt(rs)/phi + x=mu*dsqrt(rs)/phi eclr=(phi**3*Qrpa(x)+a1*mu**3+a2*mu**4+a3*mu**5+ & a4*mu**6+a5*mu**8)/((1.d0+b0**2*mu**2)**4) @@ -763,8 +801,14 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) D0f = 0.752411d0 E0f = -0.0127713d0 F0f = 0.00185898d0 + double precision :: exp_d0fx + if(dabs(D0f*x).lt.50.d0)then + exp_d0fx = dexp(-dabs(D0f)*x) + else + exp_d0fx = 0.d0 + endif g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+ & - F0f*x**4)*exp(-abs(D0f)*x)/2.d0 + F0f*x**4)*exp_d0fx/2.d0 return end @@ -778,7 +822,11 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) Dg0 = -0.0127713d0 Eg0 = 0.00185898d0 Bg0 =0.7317d0-Fg0 - expsum=exp(-Fg0*rs) + if(dabs(Fg0*rs).lt.50.d0)then + expsum=dexp(-Fg0*rs) + else + expsum = 0.d0 + endif 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))/ & @@ -795,7 +843,11 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) Dg0 = -0.0127713d0 Eg0 = 0.00185898d0 Bg0 = 0.7317d0-Fg0 - expsum=exp(-Fg0*rs) + if(dabs(Fg0*rs).lt.50.d0)then + expsum=dexp(-Fg0*rs) + else + expsum=0.d0 + endif 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* & @@ -860,19 +912,12 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) implicit none double precision pi,a2,b2,c2,d2,x,Acoul pi=dacos(-1.d0) - Acoul=2.d0*(log(2.d0)-1.d0)/pi**2 + Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2 a2 = 5.84605d0 c2 = 3.91744d0 d2 = 3.44851d0 b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0) - !if(((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)).le.0.d0)then - ! print*,(1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2) - ! print*,(1.d0+a2*x+b2*x**2+c2*x**3),(1.d0+a2*x+d2*x**2) - ! print*,x - ! pause - !endif - !Qrpa=Acoul*log(dabs((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2))) - Qrpa=Acoul*log((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)) + Qrpa=Acoul*dlog((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)) return end @@ -880,7 +925,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) implicit none double precision pi,a2,b2,c2,d2,x,Acoul pi=dacos(-1.d0) - Acoul=2.d0*(log(2.d0)-1.d0)/pi**2 + Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2 a2 = 5.84605d0 c2 = 3.91744d0 d2 = 3.44851d0 @@ -898,7 +943,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) double precision pi,a2,b2,c2,d2,x,Acoul double precision uQ,duQ,dduQ,vQ,dvQ,ddvQ pi=dacos(-1.d0) - Acoul=2.d0*(log(2.d0)-1.d0)/pi**2 + Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2 a2 = 5.84605d0 c2 = 3.91744d0 d2 = 3.44851d0 @@ -938,7 +983,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) ff=((1.d0+y)**(4.d0/3.d0)+(1.d0-y)**(4.d0/3.d0)- & 2.d0)/(2.d0**(4.d0/3.d0)-2.d0) - aaa=(1.d0-log(2.d0))/pi**2 + aaa=(1.d0-dlog(2.d0))/pi**2 call GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0, & 1.6382d0,0.49294d0,G,Gd,Gdd) ec0=G @@ -973,8 +1018,6 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) return end -! subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G,Gd) -!SCD subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G,Gd,Gdd) !SCF !cc Gd is d/drs G @@ -986,7 +1029,7 @@ subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd) double precision A,dA,ddA,B !SCF double precision sqrtx - sqrtx=sqrt(x) + sqrtx=dsqrt(x) G=-2.d0*Ac*(1.d0+alfa1*x)*dlog(1.d0+1.d0/(2.d0* & Ac*(beta1*x**0.5d0+ & beta2*x+beta3*x**1.5d0+beta4*x**2))) From 3987b9794d6e17d00860d0543d12646d3e2d63f0 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 12:52:08 +0200 Subject: [PATCH 07/11] removed a lot of floating point exceptions in DFT --- src/dft_utils_func/on_top_from_ueg.irp.f | 34 +++++++++++++++++++----- 1 file changed, 27 insertions(+), 7 deletions(-) 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..1965f37f 100644 --- a/src/dft_utils_func/on_top_from_ueg.irp.f +++ b/src/dft_utils_func/on_top_from_ueg.irp.f @@ -35,7 +35,11 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b) if (dabs(rho) > 1.d-12) 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) + 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 g0_UEG_mu_inf= 0.d0 endif @@ -67,7 +71,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.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 @@ -81,11 +89,11 @@ double precision function h_func(zeta) pi = 4d0 * datan(1d0) ahd = -0.36583d0 alpha = (4d0/(9d0*pi))**(1d0/3d0) - a1 = -(6d0*alpha/pi)*(1d0-log(2d0)) + a1 = -(6d0*alpha/pi)*(1d0-dlog(2d0)) b1 = 1.4919d0 b3 = 1.91528d0 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) end @@ -111,11 +119,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-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) - 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.50.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(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 From 4244e0b27e5d88cadd8824ef66ca469d34168ec4 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 14:52:23 +0200 Subject: [PATCH 08/11] removed FPE in g0_UEG_mu_inf, src/dft_utils_func/on_top_from_ueg.irp.f --- src/dft_utils_func/on_top_from_ueg.irp.f | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) 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 446b2846..24e40353 100644 --- a/src/dft_utils_func/on_top_from_ueg.irp.f +++ b/src/dft_utils_func/on_top_from_ueg.irp.f @@ -33,7 +33,7 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b) D = -0.01277d0 E = 0.001859d0 x = -d2*rs - if (dabs(rho) > 1.d-12.and.dabs(x).lt.20.d0) then + if (dabs(rho) > 1.d-20.and.dabs(x).lt.20.d0) then rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 x = -d2*rs if(dabs(x).lt.50.d0)then @@ -68,7 +68,11 @@ double precision function g0_UEG_mu(mu,rho_a,rho_b) C = 0.08193d0 D = -0.01277d0 E = 0.001859d0 - rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 + if(rho.gt.1.d-20)then + rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 + else + rs = (3d0 / (4d0*pi*1.d-20))**(1d0/3d0) + endif kf = (alpha*rs)**(-1d0) zeta = mu / kf x = -d2*rs*h_func(zeta)/ahd From 0850fa6f723755a0d5a725cd9a476940e14e2e3c Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 8 Jun 2020 15:17:53 +0200 Subject: [PATCH 09/11] fixed bug in g0_UEG_mu_inf, src/dft_utils_func/on_top_from_ueg.irp.f --- src/dft_utils_func/on_top_from_ueg.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 24e40353..717081a7 100644 --- a/src/dft_utils_func/on_top_from_ueg.irp.f +++ b/src/dft_utils_func/on_top_from_ueg.irp.f @@ -33,7 +33,7 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b) D = -0.01277d0 E = 0.001859d0 x = -d2*rs - if (dabs(rho) > 1.d-20.and.dabs(x).lt.20.d0) then + if (dabs(rho) > 1.d-20) then rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 x = -d2*rs if(dabs(x).lt.50.d0)then From 28605c76b9f7e582358132bd68fb287a61a68c4b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 17 Jun 2020 13:00:09 +0200 Subject: [PATCH 10/11] added all possible angular grid points as possible choices --- src/becke_numerical_grid/grid_becke.irp.f | 69 ----------------------- 1 file changed, 69 deletions(-) diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index bedf0a84..79f15c9a 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/src/becke_numerical_grid/grid_becke.irp.f @@ -39,75 +39,6 @@ BEGIN_PROVIDER [integer, n_points_grid_per_atom] END_DOC n_points_grid_per_atom = n_points_integration_angular * n_points_radial_grid -END_PROVIDER - - BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular,3) ] -&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)] - implicit none - BEGIN_DOC - ! weights and grid points for the integration on the angular variables on - ! the unit sphere centered on (0,0,0) - ! According to the LEBEDEV scheme - END_DOC - - include 'constants.include.F' - integer :: i - double precision :: accu - double precision :: degre_rad - double precision :: x(n_points_integration_angular) - double precision :: y(n_points_integration_angular) - double precision :: z(n_points_integration_angular) - double precision :: w(n_points_integration_angular) - - degre_rad = pi/180.d0 - accu = 0.d0 - - select case (n_points_integration_angular) - - case (5810) - call LD5810(X,Y,Z,W,n_points_integration_angular) - - case (2030) - call LD2030(X,Y,Z,W,n_points_integration_angular) - - case (1202) - call LD1202(X,Y,Z,W,n_points_integration_angular) - - case (0590) - call LD0590(X,Y,Z,W,n_points_integration_angular) - - case (302) - call LD0302(X,Y,Z,W,n_points_integration_angular) - - case (266) - call LD0266(X,Y,Z,W,n_points_integration_angular) - - case (194) - call LD0194(X,Y,Z,W,n_points_integration_angular) - - case (170) - call LD0170(X,Y,Z,W,n_points_integration_angular) - - case (74) - call LD0074(X,Y,Z,W,n_points_integration_angular) - - case (50) - call LD0050(X,Y,Z,W,n_points_integration_angular) - - case default - print *, irp_here//': wrong n_points_integration_angular. Expected:' - print *, '[ 50 | 74 | 170 | 194 | 266 | 302 | 590 | 1202 | 2030 | 5810 ]' - stop -1 - end select - - do i = 1, n_points_integration_angular - angular_quadrature_points(i,1) = x(i) - angular_quadrature_points(i,2) = y(i) - angular_quadrature_points(i,3) = z(i) - weights_angular_points(i) = w(i) * 4.d0 * pi - accu += w(i) - enddo - END_PROVIDER BEGIN_PROVIDER [integer , m_knowles] From 2a4497d0674a57795201ef19b1d5b82b490739c6 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 27 Jun 2020 13:31:29 +0200 Subject: [PATCH 11/11] added all the angular integration grid --- .../angular_grid_pts.irp.f | 104 ++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 src/becke_numerical_grid/angular_grid_pts.irp.f diff --git a/src/becke_numerical_grid/angular_grid_pts.irp.f b/src/becke_numerical_grid/angular_grid_pts.irp.f new file mode 100644 index 00000000..65a45977 --- /dev/null +++ b/src/becke_numerical_grid/angular_grid_pts.irp.f @@ -0,0 +1,104 @@ + + BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular,3) ] +&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)] + implicit none + BEGIN_DOC + ! weights and grid points for the integration on the angular variables on + ! the unit sphere centered on (0,0,0) + ! According to the LEBEDEV scheme + END_DOC + + include 'constants.include.F' + integer :: i + double precision :: accu + double precision :: degre_rad + double precision :: x(n_points_integration_angular) + double precision :: y(n_points_integration_angular) + double precision :: z(n_points_integration_angular) + double precision :: w(n_points_integration_angular) + + degre_rad = pi/180.d0 + accu = 0.d0 + + select case (n_points_integration_angular) + + + case (0006) + call LD0006(X,Y,Z,W,n_points_integration_angular) + case (0014) + call LD0014(X,Y,Z,W,n_points_integration_angular) + case (0026) + call LD0026(X,Y,Z,W,n_points_integration_angular) + case (0038) + call LD0038(X,Y,Z,W,n_points_integration_angular) + case (0050) + call LD0050(X,Y,Z,W,n_points_integration_angular) + case (0074) + call LD0074(X,Y,Z,W,n_points_integration_angular) + case (0086) + call LD0086(X,Y,Z,W,n_points_integration_angular) + case (0110) + call LD0110(X,Y,Z,W,n_points_integration_angular) + case (0146) + call LD0146(X,Y,Z,W,n_points_integration_angular) + case (0170) + call LD0170(X,Y,Z,W,n_points_integration_angular) + case (0194) + call LD0194(X,Y,Z,W,n_points_integration_angular) + case (0230) + call LD0230(X,Y,Z,W,n_points_integration_angular) + case (0266) + call LD0266(X,Y,Z,W,n_points_integration_angular) + case (0302) + call LD0302(X,Y,Z,W,n_points_integration_angular) + case (0350) + call LD0350(X,Y,Z,W,n_points_integration_angular) + case (0434) + call LD0434(X,Y,Z,W,n_points_integration_angular) + case (0590) + call LD0590(X,Y,Z,W,n_points_integration_angular) + case (0770) + call LD0770(X,Y,Z,W,n_points_integration_angular) + case (0974) + call LD0974(X,Y,Z,W,n_points_integration_angular) + case (1202) + call LD1202(X,Y,Z,W,n_points_integration_angular) + case (1454) + call LD1454(X,Y,Z,W,n_points_integration_angular) + case (1730) + call LD1730(X,Y,Z,W,n_points_integration_angular) + case (2030) + call LD2030(X,Y,Z,W,n_points_integration_angular) + case (2354) + call LD2354(X,Y,Z,W,n_points_integration_angular) + case (2702) + call LD2702(X,Y,Z,W,n_points_integration_angular) + case (3074) + call LD3074(X,Y,Z,W,n_points_integration_angular) + case (3470) + call LD3470(X,Y,Z,W,n_points_integration_angular) + case (3890) + call LD3890(X,Y,Z,W,n_points_integration_angular) + case (4334) + call LD4334(X,Y,Z,W,n_points_integration_angular) + case (4802) + call LD4802(X,Y,Z,W,n_points_integration_angular) + case (5294) + call LD5294(X,Y,Z,W,n_points_integration_angular) + case (5810) + call LD5810(X,Y,Z,W,n_points_integration_angular) + case default + print *, irp_here//': wrong n_points_integration_angular. See in ${QP_ROOT}/src/becke_numerical_grid/list_angular_grid to see the possible angular grid points. Ex: ' + print *, '[ 50 | 74 | 170 | 194 | 266 | 302 | 590 | 1202 | 2030 | 5810 ]' + stop -1 + end select + + do i = 1, n_points_integration_angular + angular_quadrature_points(i,1) = x(i) + angular_quadrature_points(i,2) = y(i) + angular_quadrature_points(i,3) = z(i) + weights_angular_points(i) = w(i) * 4.d0 * pi + accu += w(i) + enddo + +END_PROVIDER