diff --git a/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f b/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f index 3f271447..88fe4c85 100644 --- a/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f +++ b/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f @@ -111,8 +111,9 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a complex*16 :: accu, P_center(3) complex*16 :: d(0:n_pt_in) - complex*16 :: V_n_e_cosgtos - complex*16 :: crint_2 + complex*16, external :: V_n_e_cosgtos + complex*16, external :: crint_2 + complex*16, external :: crint_sum_2 if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. & (A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then @@ -162,22 +163,22 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a return endif - call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta & - , power_A, power_B, C_center, n_pt_in, d, n_pt_out) + call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & + power_A, power_B, C_center, n_pt_in, d, n_pt_out) if(n_pt_out < 0) then NAI_pol_mult_cosgtos = (0.d0, 0.d0) return endif - accu = (0.d0, 0.d0) - do i = 0, n_pt_out, 2 - accu += crint_2(shiftr(i, 1), const) * d(i) - -! print *, shiftr(i, 1), real(const), real(d(i)), real(crint_2(shiftr(i, 1), const)) - enddo + !accu = (0.d0, 0.d0) + !do i = 0, n_pt_out, 2 + ! accu += crint_2(shiftr(i, 1), const) * d(i) + !enddo + accu = crint_sum_2(n_pt_out, const, d) NAI_pol_mult_cosgtos = accu * coeff + return end ! --- 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 48aa0a40..df402ff1 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 @@ -35,11 +35,9 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then - !print *, ' with shwartz acc ' ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l) else - !print *, ' without shwartz acc ' dim1 = n_pt_max_integrals @@ -51,7 +49,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) ao_two_e_integral_cosgtos = 0.d0 if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then - !print *, ' not the same center' do p = 1, 3 I_power(p) = ao_power(i,p) @@ -130,14 +127,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 - !write(*,"(8(F15.7,2X))") real(integral1), real(integral2), real(integral3), real(integral4), & - ! real(integral5), real(integral6), real(integral7), real(integral8) - !write(33,"(5(F22.15,2X))") real(expo1), real(expo2), real(expo3), real(expo4), coef4*16.d0 - !write(43,"(1(F22.15,2X))") coef4 * 2.d0 * real(integral_tot) - - !integral_tot = integral1 - !print*, integral_tot - ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot) enddo ! s enddo ! r @@ -145,7 +134,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) enddo ! p else - !print *, ' the same center' do p = 1, 3 I_power(p) = ao_power(i,p) @@ -222,7 +210,7 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) endif endif -end function ao_two_e_integral_cosgtos +end ! --- @@ -258,8 +246,8 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) double precision :: thr double precision :: schwartz_ij - complex*16 :: ERI_cosgtos - complex*16 :: general_primitive_integral_cosgtos + complex*16, external :: ERI_cosgtos + complex*16, external :: general_primitive_integral_cosgtos ao_2e_cosgtos_schwartz_accel = 0.d0 @@ -273,7 +261,7 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) thr = ao_integrals_threshold*ao_integrals_threshold - allocate( schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)) ) + allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then @@ -298,45 +286,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l) expo2 = ao_expo_ord_transp_cosgtos(s,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 ) + 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 - 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 ) + 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 - 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 ) + 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 - 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 ) + 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 - 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 ) + 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) - 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 ) + 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) - 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 ) + 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) - 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 ) + 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) - 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 ) + 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) - 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 ) + 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 ) + 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, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -476,41 +464,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l) expo2 = ao_expo_ord_transp_cosgtos(s,l) - integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -530,45 +522,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) - integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -587,45 +579,45 @@ 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) - 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 @@ -639,11 +631,11 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) deallocate(schwartz_kl) -end function ao_2e_cosgtos_schwartz_accel +end ! --- -BEGIN_PROVIDER [ double precision, ao_2e_cosgtos_schwartz, (ao_num,ao_num)] +BEGIN_PROVIDER [double precision, ao_2e_cosgtos_schwartz, (ao_num, ao_num)] BEGIN_DOC ! Needed to compute Schwartz inequalities @@ -845,8 +837,6 @@ complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fac call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) accu = crint_sum_2(n_pt_out, const, d1) -! print *, n_pt_out, real(d1(0:n_pt_out)) -! print *, real(accu) general_primitive_integral_cosgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq @@ -926,7 +916,7 @@ complex*16 function ERI_cosgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_ ERI_cosgtos = I_f * coeff -end function ERI_cosgtos +end ! --- @@ -1008,7 +998,7 @@ subroutine integrale_new_cosgtos(I_f, a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_ I_f += gauleg_w(i, j) * t1(i) enddo -end subroutine integrale_new_cosgtos +end ! --- @@ -1055,7 +1045,7 @@ recursive subroutine I_x1_new_cosgtos(a, c, B_10, B_01, B_00, res, n_pt) endif -end subroutine I_x1_new_cosgtos +end ! --- @@ -1095,7 +1085,7 @@ recursive subroutine I_x2_new_cosgtos(c, B_10, B_01, B_00, res, n_pt) endif -end subroutine I_x2_new_cosgtos +end ! --- @@ -1166,7 +1156,7 @@ subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt return endif -end subroutine give_cpolynom_mult_center_x +end ! --- @@ -1209,7 +1199,7 @@ subroutine I_x1_pol_mult_cosgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt endif -end subroutine I_x1_pol_mult_cosgtos +end ! --- @@ -1302,7 +1292,7 @@ recursive subroutine I_x1_pol_mult_recurs_cosgtos(a, c, B_10, B_01, B_00, C_00, !DIR$ FORCEINLINE call multiply_cpoly(Y, ny, C_00, 2, d, nd) -end subroutine I_x1_pol_mult_recurs_cosgtos +end ! --- @@ -1358,7 +1348,7 @@ recursive subroutine I_x1_pol_mult_a1_cosgtos(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_ !DIR$ FORCEINLINE call multiply_cpoly(Y, ny, C_00, 2, d, nd) -end subroutine I_x1_pol_mult_a1_cosgtos +end ! --- @@ -1422,7 +1412,7 @@ recursive subroutine I_x1_pol_mult_a2_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d !DIR$ FORCEINLINE call multiply_cpoly(Y, ny, C_00, 2, d, nd) -end subroutine I_x1_pol_mult_a2_cosgtos +end ! --- @@ -1507,7 +1497,7 @@ recursive subroutine I_x2_pol_mult_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d, n end select -end subroutine I_x2_pol_mult_cosgtos +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 ca989255..7b4a2e5a 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -101,8 +101,6 @@ double precision function ao_two_e_integral(i, j, k, l) integral = general_primitive_integral(dim1, & P_new,P_center,fact_p,pp,p_inv,iorder_p, & Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - !write(32,"(5(F22.15,2X))") ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), coef4 - !write(42,"(1(F22.15,2X))") coef4*integral ao_two_e_integral = ao_two_e_integral + coef4 * integral enddo ! s enddo ! r diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f index 5f9e8642..9ffcc817 100644 --- a/src/utils/cpx_boys.irp.f +++ b/src/utils/cpx_boys.irp.f @@ -252,9 +252,9 @@ complex*16 function crint_smallz(n, rho) complex*16 :: rho_k, ct, delta_k ct = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0) - rho_k = (1.d0, 0.d0) - crint_smallz = ct * rho_k / gamma(dble(n) + 1.5d0) + crint_smallz = ct / gamma(dble(n) + 1.5d0) + rho_k = (1.d0, 0.d0) do k = 1, kmax rho_k = rho_k * rho @@ -269,7 +269,7 @@ complex*16 function crint_smallz(n, rho) write(*,*) ' pb in crint_smallz !' write(*,*) ' n, rho = ', n, rho write(*,*) ' delta_mod = ', delta_mod - stop 1 + !stop 1 endif end @@ -279,7 +279,6 @@ end complex*16 function crint_2(n, rho) implicit none - include 'constants.include.F' integer, intent(in) :: n complex*16, intent(in) :: rho @@ -332,166 +331,23 @@ subroutine zboysfun(n_max, x, vals) END_DOC implicit none - double precision, parameter :: tol = 1.0d-03 - double precision, parameter :: x0(1:50) = (/ 0.5d0, & - 0.8660254037844386d0, & - 1.2331060371652351d0, & - 1.6005429364718398d0, & - 1.968141713517676d0, & - 2.3358274254733784d0, & - 2.703565149324218d0, & - 3.0713364398393708d0, & - 3.4391306442915526d0, & - 3.8069411832467615d0, & - 4.174763774674144d0, & - 4.5425955121971775d0, & - 4.9104343539723185d0, & - 5.278278823606476d0, & - 5.646127827158555d0, & - 6.0139805368686865d0, & - 6.381836314862427d0, & - 6.749694661668329d0, & - 7.117555180620059d0, & - 7.48541755270593d0, & - 7.853281518456119d0, & - 8.221146864672829d0, & - 8.589013414557305d0, & - 8.956881020260978d0, & - 9.324749557193652d0, & - 9.692618919623587d0, & - 10.060489017239897d0, & - 10.428359772440366d0, & - 10.796231118172152d0, & - 11.164102996198254d0, & - 11.531975355694907d0, & - 11.899848152108484d0, & - 12.26772134621756d0, & - 12.63559490335844d0, & - 13.003468792781872d0, & - 13.371342987115684d0, & - 13.739217461913594d0, & - 14.107092195274445d0, & - 14.474967167519445d0, & - 14.842842360917235d0, & - 15.210717759448817d0, & - 15.578593348605757d0, & - 15.94646911521622d0, & - 16.31434504729452d0, & - 16.68222113391055d0, & - 17.05009736507609d0, & - 17.4179737316455d0, & - 17.78585022522868d0, & - 18.153726838114668d0, & - 18.521603563204277d0 /) - - double precision, parameter :: t(0:11) = (/ 0.20000000000000000d+01, & - 0.66666666666666663d+00, & - 0.40000000000000002d+00, & - 0.28571428571428570d+00, & - 0.22222222222222221d+00, & - 0.18181818181818182d+00, & - 0.15384615384615385d+00, & - 0.13333333333333333d+00, & - 0.11764705882352941d+00, & - 0.10526315789473684d+00, & - 0.95238095238095233d-01, & - 0.86956521739130432d-01 /) - - complex*16, parameter :: zz(1:10) = (/ ( 0.64304020652330500d+01, 0.18243694739308491d+02), & - ( 0.64304020652330500d+01,-0.18243694739308491d+02), & - (-0.12572081889410178d+01, 0.14121366415342502d+02), & - (-0.12572081889410178d+01,-0.14121366415342502d+02), & - (-0.54103079551670268d+01, 0.10457909575828442d+02), & - (-0.54103079551670268d+01,-0.10457909575828442d+02), & - (-0.78720025594983341d+01, 0.69309284623985663d+01), & - (-0.78720025594983341d+01,-0.69309284623985663d+01), & - (-0.92069621609035313d+01, 0.34559308619699376d+01), & - (-0.92069621609035313d+01,-0.34559308619699376d+01) /) - - complex*16, parameter :: fact(1:10) = (/ ( 0.13249210991966042d-02, 0.91787356295447745d-03), & - ( 0.13249210991966042d-02,-0.91787356295447745d-03), & - ( 0.55545905103006735d-01,-0.35151540664451613d+01), & - ( 0.55545905103006735d-01, 0.35151540664451613d+01), & - (-0.11456407675096416d+03, 0.19213789620924834d+03), & - (-0.11456407675096416d+03,-0.19213789620924834d+03), & - ( 0.20915556220686653d+04,-0.15825742912360638d+04), & - ( 0.20915556220686653d+04, 0.15825742912360638d+04), & - (-0.94779394228935325d+04, 0.30814443710192086d+04), & - (-0.94779394228935325d+04,-0.30814443710192086d+04) /) - - complex*16, parameter :: ww(1:10) = (/ (-0.83418049867878959d-08,-0.70958810331788253d-08), & - (-0.83418050437598581d-08, 0.70958810084577824d-08), & - ( 0.82436739552884774d-07,-0.27704117936134414d-06), & - ( 0.82436739547688584d-07, 0.27704117938414886d-06), & - ( 0.19838416382728666d-05, 0.78321058613942770d-06), & - ( 0.19838416382681279d-05,-0.78321058613180811d-06), & - (-0.47372729839268780d-05, 0.58076919074212929d-05), & - (-0.47372729839287016d-05,-0.58076919074154416d-05), & - (-0.68186014282131608d-05,-0.13515261354290787d-04), & - (-0.68186014282138385d-05, 0.13515261354295612d-04) /) - - double precision, parameter :: rzz(1:1) = (/ -0.96321934290343840d+01 /) - double precision, parameter :: rfact(1:1) = (/ 0.15247844519077540d+05 /) - double precision, parameter :: rww(1:1) = (/ 0.18995875677635889d-04 /) - - integer, intent(in) :: n_max - complex*16, intent(in) :: x + integer, intent(in) :: n_max + complex*16, intent(in) :: x complex*16, intent(out) :: vals(0:n_max) - integer :: n, k - double precision :: x0_nmax - complex*16 :: y, yy, rtmp - complex*16 :: p, q, tmp + integer :: n + complex*16 :: yy, x_inv - y = zexp(-x) -! x0_nmax = x0(n_max) + call zboysfun00(x, vals(0)) -! if(abs(x) .ge. x0_nmax) then - -!print*,'check' - call zboysfun00(x, vals(0)) -!print*, vals(0) - yy = 0.5d0 * y - do n = 1, n_max - vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - yy) / x -!print*, n, x -!print*, vals(n) - enddo - -! else -! -! rtmp = (0.d0, 0.d0) -! do k = 1, 10 -! if(abs(x + zz(k)) .ge. tol) then -! rtmp = rtmp + ww(k) * (1.0d0 - fact(k)*y) / (x + zz(k)) -! else -! q = x + zz(k) -! p = 1.0d0 - 0.5d0 * q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 -! rtmp = rtmp + ww(k) * p -! endif -! enddo -! -! tmp = (0.d0, 0.d0) -! do k = 1, 1 -! if(abs(x + rzz(k)) .ge. tol) then -! tmp = tmp + rww(k) * (1.0d0 - rfact(k)*y) / (x + rzz(k)) -! else -! q = x + rzz(k) -! p = 1.0d0 - 0.5d0 * q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 -! tmp = tmp + rww(k) * p -! endif -! enddo -! -! vals(n_max) = rtmp + tmp -! print*, vals(n_max) -! yy = 0.5d0 * y -! do n = n_max-1, 0, -1 -! vals(n) = (x * vals(n+1) + yy) * t(n) -! enddo -! -! endif + yy = 0.5d0 * zexp(-x) + x_inv = (1.d0, 0.d0) / x + do n = 1, n_max + vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - yy) * x_inv + enddo + return end ! --- @@ -510,155 +366,22 @@ subroutine zboysfunnrp(n_max, x, vals) END_DOC implicit none - double precision, parameter :: tol = 1.0d-03 - double precision, parameter :: x0(1:50) = (/ 0.5d0, & - 0.8660254037844386d0, & - 1.2331060371652351d0, & - 1.6005429364718398d0, & - 1.968141713517676d0, & - 2.3358274254733784d0, & - 2.703565149324218d0, & - 3.0713364398393708d0, & - 3.4391306442915526d0, & - 3.8069411832467615d0, & - 4.174763774674144d0, & - 4.5425955121971775d0, & - 4.9104343539723185d0, & - 5.278278823606476d0, & - 5.646127827158555d0, & - 6.0139805368686865d0, & - 6.381836314862427d0, & - 6.749694661668329d0, & - 7.117555180620059d0, & - 7.48541755270593d0, & - 7.853281518456119d0, & - 8.221146864672829d0, & - 8.589013414557305d0, & - 8.956881020260978d0, & - 9.324749557193652d0, & - 9.692618919623587d0, & - 10.060489017239897d0, & - 10.428359772440366d0, & - 10.796231118172152d0, & - 11.164102996198254d0, & - 11.531975355694907d0, & - 11.899848152108484d0, & - 12.26772134621756d0, & - 12.63559490335844d0, & - 13.003468792781872d0, & - 13.371342987115684d0, & - 13.739217461913594d0, & - 14.107092195274445d0, & - 14.474967167519445d0, & - 14.842842360917235d0, & - 15.210717759448817d0, & - 15.578593348605757d0, & - 15.94646911521622d0, & - 16.31434504729452d0, & - 16.68222113391055d0, & - 17.05009736507609d0, & - 17.4179737316455d0, & - 17.78585022522868d0, & - 18.153726838114668d0, & - 18.521603563204277d0 /) - - double precision, parameter :: t(0:11) = (/ 0.20000000000000000D+01, & - 0.66666666666666663D+00, & - 0.40000000000000002D+00, & - 0.28571428571428570D+00, & - 0.22222222222222221D+00, & - 0.18181818181818182D+00, & - 0.15384615384615385D+00, & - 0.13333333333333333D+00, & - 0.11764705882352941D+00, & - 0.10526315789473684D+00, & - 0.95238095238095233D-01, & - 0.86956521739130432D-01 /) - - complex*16, parameter :: zz(1:10) = (/ ( 0.64304020652330500D+01, 0.18243694739308491D+02), & - ( 0.64304020652330500D+01,-0.18243694739308491D+02), & - (-0.12572081889410178D+01, 0.14121366415342502D+02), & - (-0.12572081889410178D+01,-0.14121366415342502D+02), & - (-0.54103079551670268D+01, 0.10457909575828442D+02), & - (-0.54103079551670268D+01,-0.10457909575828442D+02), & - (-0.78720025594983341D+01, 0.69309284623985663D+01), & - (-0.78720025594983341D+01,-0.69309284623985663D+01), & - (-0.92069621609035313D+01, 0.34559308619699376D+01), & - (-0.92069621609035313D+01,-0.34559308619699376D+01) /) - - complex*16, parameter :: fact(1:10) = (/ ( 0.13249210991966042D-02, 0.91787356295447745D-03), & - ( 0.13249210991966042D-02,-0.91787356295447745D-03), & - ( 0.55545905103006735D-01,-0.35151540664451613D+01), & - ( 0.55545905103006735D-01, 0.35151540664451613D+01), & - (-0.11456407675096416D+03, 0.19213789620924834D+03), & - (-0.11456407675096416D+03,-0.19213789620924834D+03), & - ( 0.20915556220686653D+04,-0.15825742912360638D+04), & - ( 0.20915556220686653D+04, 0.15825742912360638D+04), & - (-0.94779394228935325D+04, 0.30814443710192086D+04), & - (-0.94779394228935325D+04,-0.30814443710192086D+04) /) - - complex*16, parameter :: ww(1:10) = (/ (-0.83418049867878959D-08,-0.70958810331788253D-08), & - (-0.83418050437598581D-08, 0.70958810084577824D-08), & - ( 0.82436739552884774D-07,-0.27704117936134414D-06), & - ( 0.82436739547688584D-07, 0.27704117938414886D-06), & - ( 0.19838416382728666D-05, 0.78321058613942770D-06), & - ( 0.19838416382681279D-05,-0.78321058613180811D-06), & - (-0.47372729839268780D-05, 0.58076919074212929D-05), & - (-0.47372729839287016D-05,-0.58076919074154416D-05), & - (-0.68186014282131608D-05,-0.13515261354290787D-04), & - (-0.68186014282138385D-05, 0.13515261354295612D-04) /) - - double precision, parameter :: rzz(1:1) = (/ -0.96321934290343840D+01 /) - double precision, parameter :: rfact(1:1) = (/ 0.15247844519077540D+05 /) - double precision, parameter :: rww(1:1) = (/ 0.18995875677635889D-04 /) - - integer, intent(in) :: n_max - complex*16, intent(in) :: x + integer, intent(in) :: n_max + complex*16, intent(in) :: x complex*16, intent(out) :: vals(0:n_max) - integer :: n, k - double precision :: x0_nmax - complex*16 :: y, rtmp - complex*16 :: p, q, tmp + integer :: n + complex*16 :: x_inv - y = exp(x) -! x0_nmax = x0(n_max) + call zboysfun00nrp(x, vals(0)) -! if(abs(x) .ge. x0_nmax) then - call zboysfun00nrp(x, vals(0)) - do n = 1, n_max - vals(n) = ((n - 0.5d0) * vals(n-1) - 0.5d0) / x - enddo - return -! endif -! -! rtmp = (0.d0, 0.d0) -! do k = 1, 10 -! if(abs(x + zz(k)) .ge. tol) then -! rtmp = rtmp + ww(k) * (y - fact(k)) / (x + zz(k)) -! else -! q = x+zz(k) -! p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 -! rtmp = rtmp + ww(k)*p -! endif -! enddo -! tmp = (0.d0, 0.d0) -! do k = 1, 1 -! if (abs(x + rzz(k)) .ge. tol) then -! tmp = tmp + rww(k)*(y-rfact(k))/(x + rzz(k)) -! else -! q = x+rzz(k) -! p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 -! tmp = tmp + rww(k) * p -! endif -! enddo -! vals(n_max) = rtmp+tmp -! do n = n_max-1, 0, -1 -! vals(n) = (x * vals(n+1) + 0.5d0) * t(n) -! enddo -! return + x_inv = (1.d0, 0.d0) / x + do n = 1, n_max + vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - 0.5d0) * x_inv + enddo + return end ! --- @@ -667,19 +390,153 @@ complex*16 function crint_sum_2(n_pt_out, rho, d1) implicit none - integer, intent(in) :: n_pt_out - complex*16, intent(in) :: rho, d1(0:n_pt_out) + integer, intent(in) :: n_pt_out + complex*16, intent(in) :: rho, d1(0:n_pt_out) + + integer :: n, i + integer :: n_max - integer :: n, i + complex*16, allocatable :: vals(:) - complex*16, external :: crint_2 + !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 - crint_sum_2 = (0.d0, 0.d0) - do i = 0, n_pt_out, 2 + n_max = shiftr(n_pt_out, 1) + + allocate(vals(0:n_max)) + call crint_2_vec(n_max, rho, vals) + + crint_sum_2 = d1(0) * vals(0) + do i = 2, n_pt_out, 2 n = shiftr(i, 1) - crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho) + crint_sum_2 += d1(i) * vals(n) enddo + deallocate(vals) + + return +end + +! --- + +subroutine crint_2_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 + + + abs_rho = abs(rho) + + if(abs_rho < 10.d0) then + + if(abs_rho .lt. 1d-6) then + + ! use finite expansion for very small rho + + ! 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 + + else + + call crint_smallz_vec(n_max, rho, vals) + + endif + + else + + if(real(rho) .ge. 0.d0) then + + call zboysfun(n_max, rho, vals) + + else + + call zboysfunnrp(n_max, rho, vals) + erho = zexp(-rho) + do n = 0, n_max + vals(n) = vals(n) * erho + enddo + + endif + + endif + + return +end + +! --- + +subroutine crint_smallz_vec(n_max, rho, vals) + + BEGIN_DOC + ! Standard version of rint + END_DOC + + implicit none + integer, intent(in) :: n_max + complex*16, intent(in) :: rho + complex*16, intent(out) :: vals(0:n_max) + + integer, parameter :: kmax = 40 + double precision, parameter :: eps = 1.d-13 + + integer :: k, n + complex*16 :: ct, delta_k + complex*16 :: rhoe + complex*16, allocatable :: rho_k(:) + + + allocate(rho_k(0:kmax)) + + rho_k(0) = (1.d0, 0.d0) + do k = 1, kmax + rho_k(k) = rho_k(k-1) * rho + enddo + + rhoe = 0.5d0 * zexp(-rho) + + do n = 0, n_max + + ct = rhoe * gamma(dble(n) + 0.5d0) + vals(n) = ct / gamma(dble(n) + 1.5d0) + + do k = 1, kmax + delta_k = ct * rho_k(k) / gamma(dble(n+k) + 1.5d0) + vals(n) += delta_k + if(abs(delta_k) .lt. eps) then + exit + endif + enddo + + !if(abs(delta_k) > eps) then + ! write(*,*) ' pb in crint_smallz_vec !' + ! write(*,*) ' n, rho = ', n, rho + ! write(*,*) ' |delta_k| = ', abs(delta_k) + !endif + enddo + + deallocate(rho_k) + + return end ! ---