From 23a242052d53a5a3226236f17b037b41361daa73 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Sat, 12 Oct 2024 16:44:41 +0200 Subject: [PATCH] few opt in cGTOs --- src/ao_two_e_ints/screening.irp.f | 2 +- .../two_e_Coul_integrals_cosgtos.irp.f | 367 ++++++++--------- src/ao_two_e_ints/two_e_integrals.irp.f | 15 +- src/hartree_fock/deb_ao_2e_int.irp.f | 194 +++++++-- src/utils/cgtos_utils.irp.f | 121 ++---- src/utils/cpx_boys.irp.f | 381 ++++++++++++++---- src/utils/cpx_erf.irp.f | 90 +++-- src/utils/integration.irp.f | 112 +++-- 8 files changed, 792 insertions(+), 490 deletions(-) diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f index d3230370..0efaaef6 100644 --- a/src/ao_two_e_ints/screening.irp.f +++ b/src/ao_two_e_ints/screening.irp.f @@ -3,7 +3,7 @@ logical function ao_two_e_integral_zero(i,j,k,l) integer, intent(in) :: i,j,k,l ao_two_e_integral_zero = .False. - if (.not.(read_ao_two_e_integrals.or.is_periodic)) then + if (.not.(read_ao_two_e_integrals.or.is_periodic.or.use_cosgtos)) then if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then ao_two_e_integral_zero = .True. return diff --git a/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f b/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f index df402ff1..5aa9fea3 100644 --- a/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f +++ b/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f @@ -11,27 +11,25 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) implicit none include 'utils/constants.include.F' - integer, intent(in) :: i, j, k, l + integer, intent(in) :: i, j, k, l + + integer :: p, q, r, s + integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3) + integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3) + double precision :: coef1, coef2, coef3, coef4 + complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3) + complex*16 :: expo1, expo2, expo3, expo4 + complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv + complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv + complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv + complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv + complex*16 :: integral1, integral2, integral3, integral4 + complex*16 :: integral5, integral6, integral7, integral8 + complex*16 :: integral_tot - integer :: p, q, r, s - integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3) - integer :: iorder_p1(3), iorder_p2(3), iorder_p3(3), iorder_p4(3), iorder_q1(3), iorder_q2(3) - double precision :: coef1, coef2, coef3, coef4 - complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3) - complex*16 :: expo1, expo2, expo3, expo4 - complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv - complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv - complex*16 :: P3_new(0:max_dim,3), P3_center(3), fact_p3, pp3, p3_inv - complex*16 :: P4_new(0:max_dim,3), P4_center(3), fact_p4, pp4, p4_inv - complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv - complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv - complex*16 :: integral1, integral2, integral3, integral4 - complex*16 :: integral5, integral6, integral7, integral8 - complex*16 :: integral_tot - - double precision :: ao_2e_cosgtos_schwartz_accel - complex*16 :: ERI_cosgtos - complex*16 :: general_primitive_integral_cosgtos + double precision, external :: ao_2e_cosgtos_schwartz_accel + complex*16, external :: ERI_cosgtos + complex*16, external :: general_primitive_integral_cosgtos if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then @@ -71,19 +69,11 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & expo1, expo2, I_power, J_power, I_center, J_center, dim1) - p1_inv = (1.d0,0.d0) / pp1 + p1_inv = (1.d0, 0.d0) / pp1 call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1) - p2_inv = (1.d0,0.d0) / pp2 - - call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, & - expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1) - p3_inv = (1.d0,0.d0) / pp3 - - call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, & - conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1) - p4_inv = (1.d0,0.d0) / pp4 + p2_inv = (1.d0, 0.d0) / pp2 do r = 1, ao_prim_num(k) coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k) @@ -95,35 +85,43 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & expo3, expo4, K_power, L_power, K_center, L_center, dim1) - q1_inv = (1.d0,0.d0) / qq1 + q1_inv = (1.d0, 0.d0) / qq1 call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1) - q2_inv = (1.d0,0.d0) / qq2 + q2_inv = (1.d0, 0.d0) / qq2 - integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + integral1 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + integral2 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + integral3 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + integral4 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + integral5 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + integral6 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + integral7 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + integral8 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -158,57 +156,57 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) expo4 = ao_expo_ord_transp_cosgtos(s,l) - integral1 = ERI_cosgtos( expo1, expo2, expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot) enddo ! s - enddo ! r - enddo ! q - enddo ! p + enddo ! r + enddo ! q + enddo ! p - endif - endif + endif ! same centers + endif ! do schwartz end @@ -228,14 +226,12 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) integer :: p, q, r, s integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3) - integer :: iorder_p1(3), iorder_p2(3), iorder_p3(3), iorder_p4(3), iorder_q1(3), iorder_q2(3) + integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3) double precision :: coef1, coef2, coef3, coef4 complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3) complex*16 :: expo1, expo2, expo3, expo4 complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv - complex*16 :: P3_new(0:max_dim,3), P3_center(3), fact_p3, pp3, p3_inv - complex*16 :: P4_new(0:max_dim,3), P4_center(3), fact_p4, pp4, p4_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv complex*16 :: integral1, integral2, integral3, integral4 @@ -288,47 +284,46 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & expo1, expo2, K_power, L_power, K_center, L_center, dim1) - p1_inv = (1.d0,0.d0) / pp1 + p1_inv = (1.d0, 0.d0) / pp1 call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1) - p2_inv = (1.d0,0.d0) / pp2 + p2_inv = (1.d0, 0.d0) / pp2 - call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, & - expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1) - p3_inv = (1.d0,0.d0) / pp3 + integral1 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, & - conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1) - p4_inv = (1.d0,0.d0) / pp4 + integral2 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) + integral3 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) + integral4 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) + integral5 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) + integral6 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) + integral7 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - - integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - - integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) + integral8 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 - schwartz_kl(s,r) = coef2 * 2.d0 * real(integral_tot) schwartz_kl(0,r) = max(schwartz_kl(0,r), schwartz_kl(s,r)) @@ -346,45 +341,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j) expo2 = ao_expo_ord_transp_cosgtos(q,j) - call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 & - , expo1, expo2, I_power, J_power, I_center, J_center, dim1 ) - p1_inv = (1.d0,0.d0) / pp1 + call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & + expo1, expo2, I_power, J_power, I_center, J_center, dim1) + p1_inv = (1.d0, 0.d0) / pp1 - call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & - , conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1 ) - p2_inv = (1.d0,0.d0) / pp2 + call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & + conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1) + p2_inv = (1.d0, 0.d0) / pp2 - call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & - , expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1 ) - p3_inv = (1.d0,0.d0) / pp3 + integral1 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & - , conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1 ) - p4_inv = (1.d0,0.d0) / pp4 + integral2 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral3 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral4 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral5 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral6 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral7 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) - - integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) - - integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral8 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -404,46 +399,53 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) expo4 = ao_expo_ord_transp_cosgtos(s,l) - call give_explicit_cpoly_and_cgaussian( Q1_new, Q1_center, qq1, fact_q1, iorder_q1 & - , expo3, expo4, K_power, L_power, K_center, L_center, dim1 ) - q1_inv = (1.d0,0.d0) / qq1 + call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & + expo3, expo4, K_power, L_power, K_center, L_center, dim1) + q1_inv = (1.d0, 0.d0) / qq1 - call give_explicit_cpoly_and_cgaussian( Q2_new, Q2_center, qq2, fact_q2, iorder_q2 & - , conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1 ) - q2_inv = (1.d0,0.d0) / qq2 + call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & + conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1) + q2_inv = (1.d0, 0.d0) / qq2 - integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral1 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral2 = general_primitive_integral_cosgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral3 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral4 = general_primitive_integral_cosgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + integral5 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral6 = general_primitive_integral_cosgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral7 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) - - integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral8 = general_primitive_integral_cosgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 ao_2e_cosgtos_schwartz_accel = ao_2e_cosgtos_schwartz_accel + coef4 * 2.d0 * real(integral_tot) enddo ! s - enddo ! r - enddo ! q - enddo ! p + enddo ! r + enddo ! q + enddo ! p else @@ -705,14 +707,6 @@ complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fac p10_2 = pq_inv_2 * p10_1 * q ! 0.5d0*q/(pq + p*p) p01_2 = pq_inv_2 * p01_1 * p ! 0.5d0*p/(q*q + pq) - ! get \sqrt(p + q) - !ppq = p + q - !ppq_re = REAL (ppq) - !ppq_im = AIMAG(ppq) - !ppq_mod = dsqrt(ppq_re*ppq_re + ppq_im*ppq_im) - !sq_ppq_re = sq_op5 * dsqrt(ppq_re + ppq_mod) - !sq_ppq_im = 0.5d0 * ppq_im / sq_ppq_re - !sq_ppq = sq_ppq_re + (0.d0, 1.d0) * sq_ppq_im sq_ppq = zsqrt(p + q) ! --- @@ -727,13 +721,13 @@ complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fac do i = 0, iorder_p(1) tmp_p = P_new(i,1) - tmp_mod = dsqrt(REAL(tmp_p)*REAL(tmp_p) + AIMAG(tmp_p)*AIMAG(tmp_p)) + tmp_mod = dsqrt(real(tmp_p)*real(tmp_p) + aimag(tmp_p)*aimag(tmp_p)) if(tmp_mod < thresh) cycle do j = 0, iorder_q(1) tmp_q = tmp_p * Q_new(j,1) - tmp_mod = dsqrt(REAL(tmp_q)*REAL(tmp_q) + AIMAG(tmp_q)*AIMAG(tmp_q)) + tmp_mod = dsqrt(real(tmp_q)*real(tmp_q) + aimag(tmp_q)*aimag(tmp_q)) if(tmp_mod < thresh) cycle !DIR$ FORCEINLINE @@ -1089,8 +1083,8 @@ end ! --- -subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt_in & - , pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, d, n_pt_out) +subroutine give_cpolynom_mult_center_x(P_center, Q_center, a_x, d_x, p, q, n_pt_in, & + pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, d, n_pt_out) BEGIN_DOC ! subroutine that returns the explicit polynom in term of the "t" @@ -1142,12 +1136,6 @@ subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt call I_x1_pol_mult_cosgtos(a_x, d_x, B10, B01, B00, C00, D00, d, n_pt1, n_pt_in) n_pt_out = n_pt1 -! print *, ' ' -! print *, a_x, d_x -! print *, real(B10), real(B01), real(B00), real(C00), real(D00) -! print *, n_pt1, real(d(0:n_pt1)) -! print *, ' ' - if(n_pt1 < 0) then n_pt_out = -1 do i = 0, n_pt_in @@ -1501,4 +1489,3 @@ end ! --- - diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 7b4a2e5a..ab0c7edc 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -44,7 +44,6 @@ double precision function ao_two_e_integral(i, j, k, l) logical, external :: do_schwartz_accel if(use_cosgtos) then - !print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l) @@ -54,17 +53,17 @@ double precision function ao_two_e_integral(i, j, k, l) else if (do_schwartz_accel(i,j,k,l)) then - ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l) + ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l) else - dim1 = n_pt_max_integrals + dim1 = n_pt_max_integrals - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral = 0.d0 + num_i = ao_nucl(i) + num_j = ao_nucl(j) + num_k = ao_nucl(k) + num_l = ao_nucl(l) + ao_two_e_integral = 0.d0 if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then do p = 1, 3 diff --git a/src/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f index 469eb654..816e34bd 100644 --- a/src/hartree_fock/deb_ao_2e_int.irp.f +++ b/src/hartree_fock/deb_ao_2e_int.irp.f @@ -1,9 +1,12 @@ program deb_ao_2e_int - !call check_ao_two_e_integral_cosgtos() - call check_crint1() + implicit none + + call check_ao_two_e_integral_cosgtos() + !call check_crint1() !call check_crint2() + !call check_crint3() end @@ -14,8 +17,8 @@ subroutine check_ao_two_e_integral_cosgtos() implicit none integer :: i, j, k, l - double precision :: tmp1, tmp2 double precision :: acc, nrm, dif + double precision :: tmp1, tmp2 double precision, external :: ao_two_e_integral double precision, external :: ao_two_e_integral_cosgtos @@ -23,31 +26,31 @@ subroutine check_ao_two_e_integral_cosgtos() acc = 0.d0 nrm = 0.d0 - i = 1 - j = 6 - k = 1 - l = 16 -! do i = 1, ao_num -! do k = 1, ao_num -! do j = 1, ao_num -! do l = 1, ao_num + i = 11 + j = 100 + k = 74 + l = 104 + ! do i = 1, ao_num + ! do k = 1, ao_num + ! do j = 1, ao_num + ! do l = 1, ao_num tmp1 = ao_two_e_integral (i, j, k, l) tmp2 = ao_two_e_integral_cosgtos(i, j, k, l) - dif = dabs(tmp1 - tmp2) - if(dif .gt. 1d-12) then + dif = abs(tmp1 - tmp2) + !if(dif .gt. 1d-10) then print*, ' error on:', i, j, k, l print*, tmp1, tmp2, dif - stop - endif + !stop + !endif -! acc += dif -! nrm += dabs(tmp1) -! enddo -! enddo -! enddo -! enddo + acc += dif + nrm += abs(tmp1) + ! enddo + ! enddo + ! enddo + ! enddo print *, ' acc (%) = ', dif * 100.d0 / nrm @@ -73,8 +76,9 @@ subroutine check_crint1() (+1d+4, +1d+4) /) complex*16 :: rho complex*16 :: int_an, int_nm + double precision, external :: rint - complex*16, external :: crint_1, crint_2, crint_quad + complex*16, external :: crint_1, crint_2 n = 10 dif_thr = 1d-7 @@ -92,9 +96,9 @@ subroutine check_crint1() acc_im = 0.d0 nrm_im = 0.d0 do i = 0, n - !int_an = crint_1 (i, rho) - int_an = crint_2 (i, rho) - int_nm = crint_quad(i, rho) + !int_an = crint_1(i, rho) + int_an = crint_2(i, rho) + call crint_quad_1(i, rho, 100000000, int_nm) dif_re = dabs(real(int_an) - real(int_nm)) dif_im = dabs(aimag(int_an) - aimag(int_nm)) @@ -185,4 +189,144 @@ end ! --- +subroutine check_crint3() + + implicit none + + integer :: i_test, n_test + integer :: nx, ny, n, n_quad + integer :: i, seed_size, clock_time + double precision :: xr(1:4), x + double precision :: yr(1:4), y + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + double precision :: delta_ref + double precision :: t1, t2, t_int1, t_int2 + complex*16 :: rho + complex*16 :: int1_old, int1_ref, int2_old, int2_ref + integer, allocatable :: seed(:) + + complex*16, external :: crint_2 + + call random_seed(size=seed_size) + allocate(seed(seed_size)) + call system_clock(count=clock_time) + seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /) + !seed = 123456789 + call random_seed(put=seed) + + + t_int1 = 0.d0 + t_int2 = 0.d0 + + n_test = 5 + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i_test = 1, n_test + + ! Re(rho) + call random_number(xr) + x = xr(1) + if(xr(2) .gt. 0.5d0) x = -x + nx = int(15.d0 * xr(3)) + if(xr(4) .gt. 0.5d0) nx = -nx + x = x * 10.d0**nx + + ! Im(rho) + call random_number(yr) + y = yr(1) + if(yr(2) .gt. 0.5d0) y = -y + ny = int(5.d0 * yr(3)) + if(yr(4) .gt. 0.5d0) ny = -ny + y = y * 10.d0**ny + + rho = x + (0.d0, 1.d0) * y + + call random_number(x) + x = 31.d0 * x + n = int(x) + !if(n.eq.0) cycle + + print*, " n = ", n + print*, " rho = ", real(rho), aimag(rho) + + call wall_time(t1) + int1_old = crint_2(n, rho) + !n_quad = 10000000 + !call crint_quad_1(n, rho, n_quad, int1_old) + !!delta_ref = 1.d0 + !!do while(delta_ref .gt. 1d-12) + !! n_quad = n_quad * 10 + !! !print*, " delta = ", delta_ref + !! !print*, " increasing n_quad to:", n_quad + !! call crint_quad_1(n, rho, n_quad, int1_ref) + !! delta_ref = abs(int1_ref - int1_old) + !! int1_old = int1_ref + !! if(n_quad .ge. 1000000000) then + !! print*, ' convergence was not reached for crint_quad_1' + !! print*, " delta = ", delta_ref + !! exit + !! endif + !!enddo + call wall_time(t2) + t_int1 = t_int1 + t2 - t1 + !print*, " n_quad for crint_quad_1:", n_quad + + call wall_time(t1) + n_quad = 10000000 + call crint_quad_12(n, rho, n_quad, int2_old) + !delta_ref = 1.d0 + !do while(delta_ref .gt. 1d-12) + ! n_quad = n_quad * 10 + ! !print*, " delta = ", delta_ref + ! !print*, " increasing n_quad to:", n_quad + ! call crint_quad_12(n, rho, n_quad, int2_ref) + ! delta_ref = abs(int2_ref - int2_old) + ! int2_old = int2_ref + ! if(n_quad .ge. 1000000000) then + ! print*, ' convergence was not reached for crint_quad_2' + ! print*, " delta = ", delta_ref + ! exit + ! endif + !enddo + call wall_time(t2) + t_int2 = t_int2 + t2 - t1 + !print*, " n_quad for crint_quad_2:", n_quad + + dif_re = dabs(real(int1_old) - real(int2_old)) + dif_im = dabs(aimag(int1_old) - aimag(int2_old)) + if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then + print*, ' important error found: ' + print*, " n = ", n + print*, " rho = ", real(rho), aimag(rho) + print*, real(int1_old), real(int2_old), dif_re + print*, aimag(int1_old), aimag(int2_old), dif_im + !stop + endif + + if((real(int1_old) /= real(int1_old)) .or. (aimag(int1_old) /= aimag(int1_old)) .or. & + (real(int2_old) /= real(int2_old)) .or. (aimag(int2_old) /= aimag(int2_old)) ) then + cycle + else + acc_re += dif_re + acc_im += dif_im + nrm_re += dabs(real(int1_old)) + nrm_im += dabs(aimag(int1_old)) + endif + enddo + + print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re + 1d-15) + print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im + 1d-15) + + print*, "crint_quad_1 wall time (sec) = ", t_int1 + print*, "crint_quad_2 wall time (sec) = ", t_int2 + + + deallocate(seed) + +end + +! --- diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index 9c25edc2..adeda59e 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -145,68 +145,6 @@ end ! --- -!subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim) -! BEGIN_DOC -! ! Transforms the product of -! ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) -! ! exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) exp(-(r-Nucl_center)^2 gama -! ! -! ! into -! ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) -! ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) -! ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) -! END_DOC -! implicit none -! include 'constants.include.F' -! integer, intent(in) :: dim -! integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) -! double precision, intent(in) :: alpha, beta, gama ! exponents -! double precision, intent(in) :: A_center(3) ! A center -! double precision, intent(in) :: B_center (3) ! B center -! double precision, intent(in) :: Nucl_center(3) ! B center -! double precision, intent(out) :: P_center(3) ! new center -! double precision, intent(out) :: p ! new exponent -! double precision, intent(out) :: fact_k ! constant factor -! double precision, intent(out) :: P_new(0:max_dim,3)! polynomial -! integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials -! -! double precision :: P_center_tmp(3) ! new center -! double precision :: p_tmp ! new exponent -! double precision :: fact_k_tmp,fact_k_bis ! constant factor -! double precision :: P_new_tmp(0:max_dim,3)! polynomial -! integer :: i,j -! double precision :: binom_func -! -! ! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp -! call give_explicit_cpoly_and_cgaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim) -! ! Then you create the new gaussian from the product of the new one per the Nuclei one -! call cgaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center) -! fact_k = fact_k_bis * fact_k_tmp -! -! ! Then you build the coefficient of the new polynom -! do i = 0, iorder(1) -! P_new(i,1) = 0.d0 -! do j = i,iorder(1) -! P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i) -! enddo -! enddo -! do i = 0, iorder(2) -! P_new(i,2) = 0.d0 -! do j = i,iorder(2) -! P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i) -! enddo -! enddo -! do i = 0, iorder(3) -! P_new(i,3) = 0.d0 -! do j = i,iorder(3) -! P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i) -! enddo -! enddo -! -!end - -! --- - subroutine cgaussian_product(a, xa, b, xb, k, p, xp) BEGIN_DOC @@ -237,7 +175,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) ab = a * b * p_inv k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3)) - tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k)) + tmp_mod = dsqrt(real(k)*real(k) + aimag(k)*aimag(k)) if(tmp_mod .gt. 40.d0) then k = (0.d0, 0.d0) xp(1:3) = (0.d0, 0.d0) @@ -245,9 +183,9 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) endif k = zexp(-k) - xp(1) = ( a * xa(1) + b * xb(1) ) * p_inv - xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv - xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv + xp(1) = (a * xa(1) + b * xb(1)) * p_inv + xp(2) = (a * xa(2) + b * xb(2)) * p_inv + xp(3) = (a * xa(3) + b * xb(3)) * p_inv end @@ -309,8 +247,6 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) integer, intent(out) :: nd integer :: ndtmp, ib, ic - double precision :: tmp_mod - complex*16 :: tmp if(ior(nc, nb) >= 0) then ! True if nc>=0 and nb>=0 continue @@ -332,9 +268,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) enddo do nd = ndtmp, 0, -1 - tmp = d(nd) - tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) - if(tmp_mod .lt. 1.d-15) cycle + if(abs(d(nd)) .lt. 1.d-15) cycle exit enddo @@ -432,47 +366,42 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b) complex*16, intent(in) :: x_A, x_P, x_B, x_Q complex*16, intent(out) :: P_A(0:a), P_B(0:b) - integer :: i, minab, maxab - complex*16 :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) + integer :: i + integer :: maxbinom + complex*16 :: pows_a(0:a+b+2), pows_b(0:a+b+2) double precision :: binom_func - if((a<0) .or. (b<0)) return + if((a < 0) .or. (b < 0)) return - maxab = max(a, b) - minab = max(min(a, b), 0) + maxbinom = size(binom_transp, 1) pows_a(0) = (1.d0, 0.d0) pows_a(1) = x_P - x_A + do i = 2, a + pows_a(i) = pows_a(i-1) * pows_a(1) + enddo pows_b(0) = (1.d0, 0.d0) pows_b(1) = x_Q - x_B - - do i = 2, maxab - pows_a(i) = pows_a(i-1) * pows_a(1) + do i = 2, b pows_b(i) = pows_b(i-1) * pows_b(1) enddo P_A(0) = pows_a(a) + do i = 1, min(a, maxbinom) + P_A(i) = binom_transp(i,a) * pows_a(a-i) + enddo + do i = maxbinom+1, a + P_A(i) = binom_func(a, i) * pows_a(a-i) + enddo + P_B(0) = pows_b(b) - - do i = 1, min(minab, 20) - P_A(i) = binom_transp(a-i,a) * pows_a(a-i) - P_B(i) = binom_transp(b-i,b) * pows_b(b-i) + do i = 1, min(b, maxbinom) + P_B(i) = binom_transp(i,b) * pows_b(b-i) enddo - - do i = minab+1, min(a, 20) - P_A(i) = binom_transp(a-i,a) * pows_a(a-i) - enddo - do i = minab+1, min(b, 20) - P_B(i) = binom_transp(b-i,b) * pows_b(b-i) - enddo - - do i = 101, a - P_A(i) = binom_func(a,a-i) * pows_a(a-i) - enddo - do i = 101, b - P_B(i) = binom_func(b,b-i) * pows_b(b-i) + do i = maxbinom+1, b + P_B(i) = binom_func(b, i) * pows_b(b-i) enddo end diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f index 9ffcc817..aa1be811 100644 --- a/src/utils/cpx_boys.irp.f +++ b/src/utils/cpx_boys.irp.f @@ -68,66 +68,6 @@ end ! --- -complex*16 function crint_quad(n, rho) - - implicit none - - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer :: i_quad, n_quad - double precision :: tmp_inv, tmp - - n_quad = 1000000000 - tmp_inv = 1.d0 / dble(n_quad) - - !crint_quad = 0.5d0 * zexp(-rho) - !do i_quad = 1, n_quad - 1 - ! tmp = tmp_inv * dble(i_quad) - ! tmp = tmp * tmp - ! crint_quad += zexp(-rho*tmp) * tmp**n - !enddo - !crint_quad = crint_quad * tmp_inv - - !crint_quad = 0.5d0 * zexp(-rho) - !do i_quad = 1, n_quad - 1 - ! tmp = tmp_inv * dble(i_quad) - ! crint_quad += zexp(-rho*tmp) * tmp**n / dsqrt(tmp) - !enddo - !crint_quad = crint_quad * 0.5d0 * tmp_inv - - ! Composite Boole's Rule - crint_quad = 7.d0 * zexp(-rho) - do i_quad = 1, n_quad - 1 - tmp = tmp_inv * dble(i_quad) - tmp = tmp * tmp - if(modulo(i_quad, 4) .eq. 0) then - crint_quad += 14.d0 * zexp(-rho*tmp) * tmp**n - else if(modulo(i_quad, 2) .eq. 0) then - crint_quad += 12.d0 * zexp(-rho*tmp) * tmp**n - else - crint_quad += 32.d0 * zexp(-rho*tmp) * tmp**n - endif - enddo - crint_quad = crint_quad * 2.d0 * tmp_inv / 45.d0 - - ! Composite Simpson's 3/8 rule - !crint_quad = zexp(-rho) - !do i_quad = 1, n_quad - 1 - ! tmp = tmp_inv * dble(i_quad) - ! tmp = tmp * tmp - ! if(modulo(i_quad, 3) .eq. 0) then - ! crint_quad += 2.d0 * zexp(-rho*tmp) * tmp**n - ! else - ! crint_quad += 3.d0 * zexp(-rho*tmp) * tmp**n - ! endif - !enddo - !crint_quad = crint_quad * 3.d0 * tmp_inv / 8.d0 - -end - -! --- - complex*16 function crint_sum_1(n_pt_out, rho, d1) implicit none @@ -283,20 +223,21 @@ complex*16 function crint_2(n, rho) integer, intent(in) :: n complex*16, intent(in) :: rho - double precision :: tmp - complex*16 :: rho2 - complex*16 :: vals(0:n) - complex*16, external :: crint_smallz + double precision :: tmp + complex*16 :: rho2 + complex*16 :: vals(0:n) + + complex*16, external :: crint_smallz if(abs(rho) < 10.d0) then if(abs(rho) .lt. 1d-6) then - tmp = 2.d0 * dble(n) + tmp = dble(n + n) rho2 = rho * rho - crint_2 = 1.d0 / (tmp + 1.d0) & - - rho / (tmp + 3.d0) & - + 0.5d0 * rho2 / (tmp + 5.d0) & - - 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0) + crint_2 = - 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0) & + + 0.5d0 * rho2 / (tmp + 5.d0) & + - rho / (tmp + 3.d0) & + + 1.d0 / (tmp + 1.d0) else crint_2 = crint_smallz(n, rho) endif @@ -328,6 +269,9 @@ subroutine zboysfun(n_max, x, vals) ! Input: x --- argument, complex*16, Re(x) >= 0 ! Output: vals --- values of the Boys function, n = 0, 1, ..., n_max ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! END_DOC implicit none @@ -363,6 +307,9 @@ subroutine zboysfunnrp(n_max, x, vals) ! Input: x --- argument, complex *16 Re(x)<=0 ! Output: vals --- values of e^z F(n,z), n = 0, 1, ..., n_max ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! END_DOC implicit none @@ -398,17 +345,12 @@ complex*16 function crint_sum_2(n_pt_out, rho, d1) complex*16, allocatable :: vals(:) - !complex*16, external :: crint_2 - !crint_sum_2 = (0.d0, 0.d0) - !do i = 0, n_pt_out, 2 - ! n = shiftr(i, 1) - ! crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho) - !enddo - n_max = shiftr(n_pt_out, 1) - allocate(vals(0:n_max)) + call crint_2_vec(n_max, rho, vals) + ! FOR DEBUG + !call crint_quad_12_vec(n_max, rho, vals) crint_sum_2 = d1(0) * vals(0) do i = 2, n_pt_out, 2 @@ -444,16 +386,22 @@ subroutine crint_2_vec(n_max, rho, vals) ! use finite expansion for very small rho - ! rho^2 / 2 - rho2 = 0.5d0 * rho * rho - ! rho^3 / 6 - rho3 = 0.3333333333333333d0 * rho * rho2 + !! rho^2 / 2 + !rho2 = 0.5d0 * rho * rho + !! rho^3 / 6 + !rho3 = 0.3333333333333333d0 * rho * rho2 + !vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3 + !do n = 1, n_max + ! tmp = 2.d0 * dble(n) + ! vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) & + ! + rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0) + !enddo - vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3 + ! TODO (last term) + vals(0) = 1.d0 + rho * (-0.3333333333333333d0 + rho*(0.1d0 - 0.047619047619047616d0*rho)) do n = 1, n_max tmp = 2.d0 * dble(n) - vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) & - + rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0) + vals(n) = 1.d0 / (tmp + 1.d0) + rho * (-1.d0/(tmp+3.d0) + rho*(0.5d0/(tmp+5.d0) - 0.3333333333333333d0*rho/(tmp+7.d0))) enddo else @@ -541,3 +489,268 @@ end ! --- +subroutine crint_quad_1(n, rho, n_quad, crint_quad) + + implicit none + + integer, intent(in) :: n, n_quad + complex*16, intent(in) :: rho + complex*16, intent(out) :: crint_quad + + integer :: i_quad + double precision :: tmp_inv, tmp0, tmp1, tmp2 + double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /) + + tmp_inv = 1.d0 / dble(n_quad) + + crint_quad = 7.d0 * zexp(-rho) + + tmp0 = 0.d0 + select case (n) + + case (0) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + + case (1) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + + case (2) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 * tmp1 + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + + case (3) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 * tmp1 * tmp1 + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + + case (4) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + tmp2 = tmp1 * tmp1 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp2 * tmp2 + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + + case default + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1**n + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + end select + +end + +! --- + +subroutine crint_quad_2(n, rho, n_quad, crint_quad) + + implicit none + + integer, intent(in) :: n, n_quad + complex*16, intent(in) :: rho + complex*16, intent(out) :: crint_quad + + integer :: i_quad + double precision :: tmp_inv, tmp0, tmp1, tmp2 + double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /) + complex*16 :: rhoc, rhoe + + tmp_inv = 1.d0 / dble(n_quad) + + crint_quad = 7.d0 * zexp(-rho) + + tmp0 = 0.d0 + rhoc = zexp(-rho*tmp_inv) + rhoe = (1.d0, 0.d0) + select case (n) + + case (0) + !do i_quad = 1, n_quad - 1 + ! tmp0 = tmp0 + tmp_inv + ! rhoe = rhoe * rhoc + ! tmp1 = (rhoe - 1.d0) / dsqrt(tmp0) + ! crint_quad = crint_quad + coef(iand(i_quad, 3)) * tmp1 + !enddo + !crint_quad = 1.d0 + crint_quad * 0.022222222222222223d0 * tmp_inv + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe / dsqrt(tmp0) + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (1) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (2) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (3) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 * tmp0 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (4) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + tmp2 = tmp1 * tmp1 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp2 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case default + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0**n / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + end select + +end + +! --- + +subroutine crint_quad_12(n, rho, n_quad, crint_quad) + + implicit none + + integer, intent(in) :: n, n_quad + complex*16, intent(in) :: rho + complex*16, intent(out) :: crint_quad + + integer :: i_quad + double precision :: tmp_inv, tmp0, tmp1, tmp2 + double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /) + complex*16 :: rhoc, rhoe + + tmp_inv = 1.d0 / dble(n_quad) + + crint_quad = 7.d0 * zexp(-rho) + + tmp0 = 0.d0 + rhoc = zexp(-rho*tmp_inv) + rhoe = (1.d0, 0.d0) + select case (n) + + case (0) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) + enddo + crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv + + case (1) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (2) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (3) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 * tmp0 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case (4) + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0 * tmp0 + tmp2 = tmp1 * tmp1 / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp2 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + case default + do i_quad = 1, n_quad - 1 + tmp0 = tmp0 + tmp_inv + tmp1 = tmp0**n / dsqrt(tmp0) + rhoe = rhoe * rhoc + crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1 + enddo + crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv + + end select + +end + +! --- + +subroutine crint_quad_12_vec(n_max, rho, vals) + + implicit none + + integer, intent(in) :: n_max + complex*16, intent(in) :: rho + complex*16, intent(out) :: vals(0:n_max) + + integer :: n + double precision :: tmp, abs_rho + complex*16 :: rho2, rho3, erho + + do n = 1, n_max + call crint_quad_12(n, rho, 10000000, vals(n)) + enddo + + return +end + +! --- + diff --git a/src/utils/cpx_erf.irp.f b/src/utils/cpx_erf.irp.f index 1c5fa61d..e1e8bfa6 100644 --- a/src/utils/cpx_erf.irp.f +++ b/src/utils/cpx_erf.irp.f @@ -143,7 +143,7 @@ complex*16 function erf_G(x, yabs) idble = dble(i) tmp0 = 0.5d0 * idble tmp1 = tmp0 * tmp0 + x2 - tmp2 = dexp( idble * yabs - tmp1 - dlog(tmp1) - dlog_2pi) * (x - (0.d0, 1.d0)*tmp0) + tmp2 = dexp(idble * yabs - tmp1 - dlog(tmp1) - dlog_2pi) * (x - (0.d0, 1.d0)*tmp0) erf_G = erf_G + tmp2 @@ -340,8 +340,8 @@ subroutine zboysfun00(z, val) if(abs(z) .ge. 100.0d0) then ! large |z| - z1 = 1.0d0 / zsqrt(z) - y = 1.0d0 / z + z1 = (1.0d0, 0.d0) / zsqrt(z) + y = (1.0d0, 0.d0) / z val = asymcoef(7) do k = 6, 1, -1 val = val * y + asymcoef(k) @@ -389,6 +389,8 @@ subroutine zboysfun00nrp(z, val) ! END_DOC + include 'constants.include.F' + implicit none double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & @@ -413,7 +415,6 @@ subroutine zboysfun00nrp(z, val) double precision, parameter :: tol = 1.0d-03 double precision, parameter :: sqpio2 = 0.886226925452758014d0 ! sqrt(pi)/2 - double precision, parameter :: pi = 3.14159265358979324d0 double precision, parameter :: etmax = 25.7903399171930621d0 double precision, parameter :: etmax1 = 26.7903399171930621d0 complex*16, parameter :: ima = (0.d0, 1.d0) @@ -452,40 +453,40 @@ subroutine zboysfun00nrp(z, val) 0.03112676196932382d0, & 0.013576229705876844d0 /) - double precision, parameter :: qq (1:16) = (/ 0.0007243228510223928d0, & - 0.01980651726441906d0, & - 0.11641097769229371d0, & - 0.38573968881461146d0, & - 0.9414671037609641d0, & - 1.8939510935716377d0, & - 3.3275564293459383d0, & - 5.280587297262129d0, & - 7.730992222360452d0, & - 10.590207725831563d0, & - 13.706359477128965d0, & - 16.876705473663804d0, & - 19.867876155236257d0, & - 22.441333930203022d0, & - 24.380717439613566d0, & - 25.51771075067431d0 /) + double precision, parameter :: qq(1:16) = (/ 0.0007243228510223928d0, & + 0.01980651726441906d0, & + 0.11641097769229371d0, & + 0.38573968881461146d0, & + 0.9414671037609641d0, & + 1.8939510935716377d0, & + 3.3275564293459383d0, & + 5.280587297262129d0, & + 7.730992222360452d0, & + 10.590207725831563d0, & + 13.706359477128965d0, & + 16.876705473663804d0, & + 19.867876155236257d0, & + 22.441333930203022d0, & + 24.380717439613566d0, & + 25.51771075067431d0 /) - double precision, parameter :: qq1 (1:16) = (/ 0.0007524078957852004d0,& - 0.020574499281252233d0, & - 0.12092472113522865d0, & - 0.40069643967765295d0, & - 0.9779717449089211d0, & - 1.9673875468969015d0, & - 3.4565797939091802d0, & - 5.485337886599723d0, & - 8.030755321535683d0, & - 11.000834641174064d0, & - 14.237812708111456d0, & - 17.531086359214406d0, & - 20.6382373144543d0, & - 23.31147887603379d0, & - 25.326060444703632d0, & - 26.507139770710722d0 /) + double precision, parameter :: qq1(1:16) = (/ 0.0007524078957852004d0,& + 0.020574499281252233d0, & + 0.12092472113522865d0, & + 0.40069643967765295d0, & + 0.9779717449089211d0, & + 1.9673875468969015d0, & + 3.4565797939091802d0, & + 5.485337886599723d0, & + 8.030755321535683d0, & + 11.000834641174064d0, & + 14.237812708111456d0, & + 17.531086359214406d0, & + 20.6382373144543d0, & + 23.31147887603379d0, & + 25.326060444703632d0, & + 26.507139770710722d0 /) double precision, parameter :: uu(1:16) = (/ 0.9992759394074501d0, & 0.9803883431758104d0, & @@ -532,8 +533,8 @@ subroutine zboysfun00nrp(z, val) if(abs(z) .ge. 100.0d0) then ! large |z| - z1 = 1.0d0 / zsqrt(z) - y = 1.0d0 / z + z1 = (1.0d0, 0.d0) / zsqrt(z) + y = (1.0d0, 0.d0) / z val = asymcoef(7) do k = 6, 1, -1 val = val * y + asymcoef(k) @@ -560,13 +561,13 @@ subroutine zboysfun00nrp(z, val) zsum = zsum + ww(k) * (zz - uu(k)) / (qq(k) + z) else q = z + qq(k) - p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 - zsum = zsum + ww(k) * p *zz + p = q * (0.041666666666666664d0*q * (q * (0.2d0*q - 1.d0) + 4.d0) - 0.5d0) + 1.d0 + zsum = zsum + ww(k) * p * zz endif enddo - zt = ima * sqrt(z / etmax) + zt = ima * zsqrt(z / etmax) tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) - val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z) + val = dsqrt(etmax) * zsum * inv_sq_pi + zz * tmp / zsqrt(pi*z) else zsum = (0.d0, 0.d0) do k = 1, 16 @@ -574,13 +575,14 @@ subroutine zboysfun00nrp(z, val) zsum = zsum + ww(k) * (zz - uu1(k)) / (qq1(k) + z) else q = z + qq1(k) - p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 + !p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 + p = q * (0.041666666666666664d0*q * (q * (0.2d0*q - 1.d0) + 4.d0) - 0.5d0) + 1.d0 zsum = zsum + ww(k) * p * zz endif enddo zt = ima * zsqrt(z / etmax1) tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) - val = dsqrt(etmax1) * zsum / dsqrt(pi) + zz * tmp / zsqrt(pi*z) + val = dsqrt(etmax1) * zsum * inv_sq_pi + zz * tmp / zsqrt(pi*z) endif return diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index 72029c73..4bbc89cb 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -1008,50 +1008,78 @@ subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, n_points) deallocate(pows_a) !deallocate(pows_b) -end subroutine recentered_poly2_v0 - -subroutine recentered_poly2(P_new,x_A,x_P,a,P_new2,x_B,x_Q,b) - implicit none - BEGIN_DOC - ! Recenter two polynomials - END_DOC - integer, intent(in) :: a,b - double precision, intent(in) :: x_A,x_P,x_B,x_Q - double precision, intent(out) :: P_new(0:a),P_new2(0:b) - double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) - double precision :: binom_func - integer :: i,j,k,l, minab, maxab - if ((a<0).or.(b<0) ) return - maxab = max(a,b) - minab = max(min(a,b),0) - pows_a(0) = 1.d0 - pows_a(1) = (x_P - x_A) - pows_b(0) = 1.d0 - pows_b(1) = (x_Q - x_B) - do i = 2,maxab - pows_a(i) = pows_a(i-1)*pows_a(1) - pows_b(i) = pows_b(i-1)*pows_b(1) - enddo - P_new (0) = pows_a(a) - P_new2(0) = pows_b(b) - do i = 1,min(minab,20) - P_new (i) = binom_transp(a-i,a) * pows_a(a-i) - P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) - enddo - do i = minab+1,min(a,20) - P_new (i) = binom_transp(a-i,a) * pows_a(a-i) - enddo - do i = minab+1,min(b,20) - P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) - enddo - do i = 101,a - P_new(i) = binom_func(a,a-i) * pows_a(a-i) - enddo - do i = 101,b - P_new2(i) = binom_func(b,b-i) * pows_b(b-i) - enddo end +! --- + +subroutine recentered_poly2(P_new, x_A, x_P, a, Q_new, x_B, x_Q, b) + + BEGIN_DOC + ! + ! Recenter two polynomials: + ! + ! (x - x_A)^a -> \sum_{i=0}^{a} \binom{a}{i} (x_A - x_P)^{a-i} (x - x_P)^i + ! (x - x_B)^b -> \sum_{i=0}^{b} \binom{b}{i} (x_B - x_Q)^{b-i} (x - x_Q)^i + ! + ! where: + ! \binom{a}{i} = \binom{a}{a-i} = a! / [i! (a-i)!] + ! + END_DOC + + implicit none + + integer, intent(in) :: a, b + double precision, intent(in) :: x_A, x_P, x_B, x_Q + double precision, intent(out) :: P_new(0:a), Q_new(0:b) + + double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) + integer :: i, minab, maxab + integer :: maxbinom + + double precision, external :: binom_func + + if((a < 0) .or. (b < 0)) return + + maxbinom = size(binom_transp, 1) + + maxab = max(a, b) + minab = min(min(a, b), maxbinom) + + pows_a(0) = 1.d0 + pows_a(1) = x_P - x_A + pows_b(0) = 1.d0 + pows_b(1) = x_Q - x_B + + do i = 2, maxab + pows_a(i) = pows_a(i-1) * pows_a(1) + pows_b(i) = pows_b(i-1) * pows_b(1) + enddo + + P_new(0) = pows_a(a) + Q_new(0) = pows_b(b) + do i = 1, minab + P_new(i) = binom_transp(i,a) * pows_a(a-i) + Q_new(i) = binom_transp(i,b) * pows_b(b-i) + enddo + + do i = minab+1, min(a, maxbinom) + P_new(i) = binom_transp(i,a) * pows_a(a-i) + enddo + do i = minab+1, min(b, maxbinom) + Q_new(i) = binom_transp(i,b) * pows_b(b-i) + enddo + + do i = maxbinom+1, a + P_new(i) = binom_func(a, i) * pows_a(a-i) + enddo + do i = maxbinom+1, b + Q_new(i) = binom_func(b, i) * pows_b(b-i) + enddo + +end + +! --- + subroutine pol_modif_center(A_center, B_center, iorder, A_pol, B_pol) BEGIN_DOC