From 7f570517390420b577105eda58751629e17754f4 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Wed, 9 Oct 2024 01:07:57 +0200 Subject: [PATCH] Fixed cplx Boys fct for real arguments --- .../one_e_Coul_integrals_cosgtos.irp.f | 20 +- .../two_e_Coul_integrals_cosgtos.irp.f | 150 ++-- src/ao_two_e_ints/two_e_integrals.irp.f | 2 + src/hartree_fock/deb_ao_2e_int.irp.f | 188 +++++ src/hartree_fock/print_scf_int.irp.f | 114 --- src/utils/cgtos_utils.irp.f | 278 +------ src/utils/cpx_boys.irp.f | 686 ++++++++++++++++++ src/utils/cpx_erf.irp.f | 412 ++++++++++- 8 files changed, 1335 insertions(+), 515 deletions(-) create mode 100644 src/hartree_fock/deb_ao_2e_int.irp.f delete mode 100644 src/hartree_fock/print_scf_int.irp.f create mode 100644 src/utils/cpx_boys.irp.f 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 7f94f226..3f271447 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 @@ -112,7 +112,7 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a complex*16 :: d(0:n_pt_in) complex*16 :: V_n_e_cosgtos - complex*16 :: crint + complex*16 :: crint_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 @@ -158,7 +158,7 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) ) if(n_pt == 0) then - NAI_pol_mult_cosgtos = coeff * crint(0, const) + NAI_pol_mult_cosgtos = coeff * crint_2(0, const) return endif @@ -172,13 +172,13 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a accu = (0.d0, 0.d0) do i = 0, n_pt_out, 2 - accu += crint(shiftr(i, 1), const) * d(i) + accu += crint_2(shiftr(i, 1), const) * d(i) -! print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const)) +! print *, shiftr(i, 1), real(const), real(d(i)), real(crint_2(shiftr(i, 1), const)) enddo NAI_pol_mult_cosgtos = accu * coeff -end function NAI_pol_mult_cosgtos +end ! --- @@ -312,7 +312,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta & d(i) = d1(i) enddo -end subroutine give_cpolynomial_mult_center_one_e +end ! --- @@ -405,7 +405,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_ endif -end subroutine I_x1_pol_mult_one_e_cosgtos +end ! --- @@ -467,7 +467,7 @@ recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim) endif -end subroutine I_x2_pol_mult_one_e_cosgtos +end ! --- @@ -502,7 +502,7 @@ complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta) * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1) endif -end function V_n_e_cosgtos +end ! --- @@ -529,7 +529,7 @@ complex*16 function V_r_cosgtos(n, alpha) V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1) endif -end function V_r_cosgtos +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 ea9ff009..48aa0a40 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 @@ -72,72 +72,22 @@ double precision function ao_two_e_integral_cosgtos(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 ) + 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 ) + 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 ) + 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 ) + 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 - !integer :: ii - !do ii = 1, 3 - ! print *, 'fact_p1', fact_p1 - ! print *, 'fact_p2', fact_p2 - ! print *, 'fact_p3', fact_p3 - ! print *, 'fact_p4', fact_p4 - ! !print *, pp1, p1_inv - ! !print *, pp2, p2_inv - ! !print *, pp3, p3_inv - ! !print *, pp4, p4_inv - !enddo - ! if( abs(aimag(P1_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_1 is complex !!' - ! print *, P1_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - ! if( abs(aimag(P2_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_2 is complex !!' - ! print *, P2_center - ! print *, ' old expos:' - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! print *, ' new expo:' - ! print *, pp2, p2_inv - ! print *, ' factor:' - ! print *, fact_p2 - ! print *, ' old centers:' - ! print *, I_center, J_center - ! print *, ' powers:' - ! print *, I_power, J_power - ! stop - ! endif - ! if( abs(aimag(P3_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_3 is complex !!' - ! print *, P3_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - ! if( abs(aimag(P4_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_4 is complex !!' - ! print *, P4_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - !enddo - do r = 1, ao_prim_num(k) coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k) expo3 = ao_expo_ord_transp_cosgtos(r,k) @@ -146,71 +96,53 @@ 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) - 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 ) + 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 ) + 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 - !do ii = 1, 3 - ! !print *, qq1, q1_inv - ! !print *, qq2, q2_inv - ! print *, 'fact_q1', fact_q1 - ! print *, 'fact_q2', fact_q2 - !enddo - ! if( abs(aimag(Q1_center(ii))) .gt. 0.d0 ) then - ! print *, ' Q_1 is complex !!' - ! print *, Q1_center - ! print *, expo3, expo4 - ! print *, conjg(expo3), conjg(expo4) - ! stop - ! endif - ! if( abs(aimag(Q2_center(ii))) .gt. 0.d0 ) then - ! print *, ' Q_2 is complex !!' - ! print *, Q2_center - ! print *, expo3, expo4 - ! print *, conjg(expo3), conjg(expo4) - ! stop - ! endif - !enddo + 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) - 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 ) + 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) - 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 ) + 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) - 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 ) + 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) - 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 ) + 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) - 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 ) + 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) - 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, 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, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) 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 - enddo ! q - enddo ! p + enddo ! r + enddo ! q + enddo ! p else !print *, ' the same center' @@ -739,8 +671,8 @@ END_PROVIDER ! --- -complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fact_p, p, p_inv, iorder_p & - , Q_new, Q_center, fact_q, q, q_inv, iorder_q ) +complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fact_p, p, p_inv, iorder_p, & + Q_new, Q_center, fact_q, q, q_inv, iorder_q) BEGIN_DOC ! @@ -765,7 +697,7 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa complex*16 :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim) complex*16 :: d1(0:max_dim), d_poly(0:max_dim) - complex*16 :: crint_sum + complex*16 :: crint_sum_2 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol @@ -912,13 +844,13 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa !DIR$ FORCEINLINE call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) - accu = crint_sum(n_pt_out, const, d1) + 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 -end function general_primitive_integral_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 7b4a2e5a..ca989255 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -101,6 +101,8 @@ 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/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f new file mode 100644 index 00000000..469eb654 --- /dev/null +++ b/src/hartree_fock/deb_ao_2e_int.irp.f @@ -0,0 +1,188 @@ + +program deb_ao_2e_int + + !call check_ao_two_e_integral_cosgtos() + call check_crint1() + !call check_crint2() + +end + +! --- + +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, external :: ao_two_e_integral + double precision, external :: 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 + + 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 + print*, ' error on:', i, j, k, l + print*, tmp1, tmp2, dif + stop + endif + +! acc += dif +! nrm += dabs(tmp1) +! enddo +! enddo +! enddo +! enddo + + print *, ' acc (%) = ', dif * 100.d0 / nrm + +end + +! --- + +subroutine check_crint1() + + implicit none + integer :: i, n, i_rho + double precision :: dif_thr + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), & + (+1d-9, +1d-6), & + (-1d-6, -1d-5), & + (+1d-3, -1d-2), & + (-1d-1, +1d-1), & + (+1d-0, +1d-1), & + (-1d+1, +1d+1), & + (+1d+2, +1d+1), & + (-1d+3, +1d+2), & + (+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 + + n = 10 + dif_thr = 1d-7 + + do i_rho = 8, 10 + !do i_rho = 7, 7 + + !rho = (-10.d0, 0.1d0) + !rho = (+10.d0, 0.1d0) + rho = rho_test(i_rho) + print*, "rho = ", real(rho), aimag(rho) + + acc_re = 0.d0 + nrm_re = 0.d0 + 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) + + dif_re = dabs(real(int_an) - real(int_nm)) + dif_im = dabs(aimag(int_an) - aimag(int_nm)) + + if((dif_re .gt. dif_thr) .or. (dif_im .gt. dif_thr)) then + print*, ' error on i =', i + print*, real(int_an), real(int_nm), dif_re + print*, aimag(int_an), aimag(int_nm), dif_im + !print*, rint(i, real(rho)) + print*, crint_1(i, rho) + !print*, crint_2(i, rho) + stop + endif + acc_re += dif_re + nrm_re += dabs(real(int_nm)) + acc_im += dif_im + nrm_im += dabs(aimag(int_nm)) + 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) + enddo + +end + +! --- + +subroutine check_crint2() + + implicit none + + integer :: i, n, i_rho + double precision :: dif_thr + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), & + (+1d-9, +1d-6), & + (-1d-6, -1d-5), & + (+1d-3, -1d-2), & + (-1d-1, +1d-1), & + (+1d-0, +1d-1), & + (-1d+1, +1d+1), & + (+1d+2, +1d+1), & + (-1d+3, +1d+2), & + (+1d+4, +1d+4) /) + complex*16 :: rho + complex*16 :: int_an, int_nm + complex*16, external :: crint_1, crint_2 + + n = 30 + dif_thr = 1d-12 + + do i_rho = 1, 10 + rho = rho_test(i_rho) + print*, "rho = ", real(rho), aimag(rho) + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i = 0, n + int_an = crint_1(i, rho) + int_nm = crint_2(i, rho) + + dif_re = dabs(real(int_an) - real(int_nm)) + !if(dif_re .gt. dif_thr) then + ! print*, ' error in real part:', i + ! print*, real(int_an), real(int_nm), dif_re + ! stop + !endif + acc_re += dif_re + nrm_re += dabs(real(int_nm)) + + dif_im = dabs(aimag(int_an) - aimag(int_nm)) + !if(dif_im .gt. dif_thr) then + ! print*, ' error in imag part:', i + ! print*, aimag(int_an), aimag(int_nm), dif_im + ! stop + !endif + acc_im += dif_im + nrm_im += dabs(aimag(int_nm)) + 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) + enddo + +end + +! --- + + diff --git a/src/hartree_fock/print_scf_int.irp.f b/src/hartree_fock/print_scf_int.irp.f deleted file mode 100644 index ee7590f6..00000000 --- a/src/hartree_fock/print_scf_int.irp.f +++ /dev/null @@ -1,114 +0,0 @@ - -program print_scf_int - - call main() - -end - -subroutine main() - - implicit none - integer :: i, j, k, l - - print *, " Hcore:" - do j = 1, ao_num - do i = 1, ao_num - print *, i, j, ao_one_e_integrals(i,j) - enddo - enddo - - print *, " P:" - do j = 1, ao_num - do i = 1, ao_num - print *, i, j, SCF_density_matrix_ao_alpha(i,j) - enddo - enddo - - - double precision :: integ, density_a, density_b, density - double precision :: J_scf(ao_num, ao_num) - double precision :: K_scf(ao_num, ao_num) - - - double precision, external :: get_ao_two_e_integral - PROVIDE ao_integrals_map - - print *, " J:" - !do j = 1, ao_num - ! do l = 1, ao_num - ! do i = 1, ao_num - ! do k = 1, ao_num - ! ! < 1:k, 2:l | 1:i, 2:j > - ! print *, '< k l | i j >', k, l, i, j - ! print *, get_ao_two_e_integral(i, j, k, l, ao_integrals_map) - ! enddo - ! enddo - ! enddo - !enddo - - !do k = 1, ao_num - ! do i = 1, ao_num - ! do j = 1, ao_num - ! do l = 1, ao_num - ! ! ( 1:k, 1:i | 2:l, 2:j ) - ! print *, '(k i | l j)', k, i, l, j - ! print *, get_ao_two_e_integral(l, j, k, i, ao_integrals_map) - ! enddo - ! enddo - ! print *, '' - ! enddo - !enddo - - J_scf = 0.d0 - K_scf = 0.d0 - do i = 1, ao_num - do k = 1, ao_num - do j = 1, ao_num - do l = 1, ao_num - - density_a = SCF_density_matrix_ao_alpha(l,j) - density_b = SCF_density_matrix_ao_beta (l,j) - density = density_a + density_b - - integ = get_ao_two_e_integral(l, j, k, i, ao_integrals_map) - J_scf(k,i) += density * integ - integ = get_ao_two_e_integral(l, i, k, j, ao_integrals_map) - K_scf(k,i) -= density_a * integ - enddo - enddo - enddo - enddo - - print *, 'J x P' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, J_scf(k,i) - enddo - enddo - - print *, '' - print *, 'K x P' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, K_scf(k,i) - enddo - enddo - - print *, '' - print *, 'F in AO' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, Fock_matrix_ao(k,i) - enddo - enddo - - print *, '' - print *, 'F in MO' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, 2.d0 * Fock_matrix_mo_alpha(k,i) - enddo - enddo - -end - diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index a820d5f2..9c25edc2 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -56,7 +56,7 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new) iorder = a + b -end subroutine give_explicit_cpoly_and_cgaussian_x +end ! --- @@ -141,7 +141,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, !DIR$ FORCEINLINE call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new) -end subroutine give_explicit_cpoly_and_cgaussian +end ! --- @@ -249,7 +249,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv -end subroutine cgaussian_product +end ! --- @@ -290,7 +290,7 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) k = zexp(-k) xp = (a*xa + b*xb) * p_inv -end subroutine cgaussian_product_x +end ! --- @@ -338,7 +338,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) exit enddo -end subroutine multiply_cpoly +end ! --- @@ -373,7 +373,7 @@ subroutine add_cpoly(b, nb, c, nc, d, nd) if(nd < 0) exit enddo -end subroutine add_cpoly +end ! --- @@ -413,7 +413,7 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) endif -end subroutine add_cpoly_multiply +end ! --- @@ -475,7 +475,7 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b) P_B(i) = binom_func(b,b-i) * pows_b(b-i) enddo -end subroutine recentered_cpoly2 +end ! --- @@ -514,267 +514,7 @@ complex*16 function Fc_integral(n, inv_sq_p) Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1)) -end function Fc_integral - -! --- - -complex*16 function crint(n, rho) - - implicit none - include 'constants.include.F' - - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer :: i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - double precision :: n_tmp - complex*16 :: sq_rho, rho_inv, rho_exp - - complex*16 :: crint_smallz, cpx_erf - - rho_re = REAL (rho) - rho_im = AIMAG(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - - if(rho_mod < 10.d0) then - ! small z - - if(rho_mod .lt. 1.d-10) then - crint = 1.d0 / dble(n + n + 1) - else - crint = crint_smallz(n, rho) - endif - - else - ! large z - - if(rho_mod .gt. 40.d0) then - - n_tmp = dble(n) + 0.5d0 - crint = 0.5d0 * gamma(n_tmp) / (rho**n_tmp) - - else - - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - crint = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho - mmax = n - if(mmax .gt. 0) then - do i = 0, mmax-1 - crint = ((dble(i) + 0.5d0) * crint - rho_exp) * rho_inv - enddo - endif - - ! *** - - endif - - endif - -! print *, n, real(rho), real(crint) - -end function crint - -! --- - -complex*16 function crint_sum(n_pt_out, rho, d1) - - implicit none - include 'constants.include.F' - - integer, intent(in) :: n_pt_out - complex*16, intent(in) :: rho, d1(0:n_pt_out) - - integer :: n, i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - complex*16 :: sq_rho, F0 - complex*16 :: rho_tmp, rho_inv, rho_exp - complex*16, allocatable :: Fm(:) - - complex*16 :: crint_smallz, cpx_erf - - rho_re = REAL (rho) - rho_im = AIMAG(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - - if(rho_mod < 10.d0) then - ! small z - - if(rho_mod .lt. 1.d-10) then - -! print *, ' 111' -! print *, ' rho = ', rho - - crint_sum = d1(0) -! print *, 0, 1 - - do i = 2, n_pt_out, 2 - - n = shiftr(i, 1) - crint_sum = crint_sum + d1(i) / dble(n+n+1) - -! print *, n, 1.d0 / dble(n+n+1) - enddo - - ! *** - - else - -! print *, ' 222' -! print *, ' rho = ', real(rho) -! if(abs(aimag(rho)) .gt. 1d-15) then -! print *, ' complex rho', rho -! stop -! endif - - crint_sum = d1(0) * crint_smallz(0, rho) - -! print *, 0, real(d1(0)), real(crint_smallz(0, rho)) -! if(abs(aimag(d1(0))) .gt. 1d-15) then -! print *, ' complex d1(0)', d1(0) -! stop -! endif - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum = crint_sum + d1(i) * crint_smallz(n, rho) - -! print *, n, real(d1(i)), real(crint_smallz(n, rho)) -! if(abs(aimag(d1(i))) .gt. 1d-15) then -! print *, ' complex d1(i)', i, d1(i) -! stop -! endif - - enddo - -! print *, 'sum = ', real(crint_sum) -! if(abs(aimag(crint_sum)) .gt. 1d-15) then -! print *, ' complex crint_sum', crint_sum -! stop -! endif - - ! *** - - endif - - else - ! large z - - if(rho_mod .gt. 40.d0) then - -! print *, ' 333' -! print *, ' rho = ', rho - - rho_inv = (1.d0, 0.d0) / rho - rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv) - crint_sum = rho_tmp * d1(0) -! print *, 0, rho_tmp - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv - crint_sum = crint_sum + rho_tmp * d1(i) -! print *, n, rho_tmp - enddo - - ! *** - - else - -! print *, ' 444' -! print *, ' rho = ', rho - - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - !sq_rho = zsqrt(rho) - - - F0 = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho - crint_sum = F0 * d1(0) -! print *, 0, F0 - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - mmax = shiftr(n_pt_out, 1) - if(mmax .gt. 0) then - - allocate( Fm(mmax) ) - Fm(1:mmax) = (0.d0, 0.d0) - - do n = 0, mmax-1 - F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv - Fm(n+1) = F0 -! print *, n, F0 - enddo - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum = crint_sum + Fm(n) * d1(i) - enddo - deallocate(Fm) - endif - - ! *** - - endif - - endif - -end function crint_sum - -! --- - -complex*16 function crint_smallz(n, rho) - - BEGIN_DOC - ! Standard version of rint - END_DOC - - implicit none - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer, parameter :: kmax = 40 - double precision, parameter :: eps = 1.d-13 - - integer :: k - double precision :: delta_mod - 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) - - do k = 1, kmax - - rho_k = rho_k * rho - delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0) - crint_smallz = crint_smallz + delta_k - - delta_mod = dsqrt(REAL(delta_k)*REAL(delta_k) + AIMAG(delta_k)*AIMAG(delta_k)) - if(delta_mod .lt. eps) return - enddo - - if(delta_mod > eps) then - write(*,*) ' pb in crint_smallz !' - write(*,*) ' n, rho = ', n, rho - write(*,*) ' delta_mod = ', delta_mod - stop 1 - endif - -end function crint_smallz +end ! --- diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f new file mode 100644 index 00000000..5f9e8642 --- /dev/null +++ b/src/utils/cpx_boys.irp.f @@ -0,0 +1,686 @@ + +! --- + +complex*16 function crint_1(n, rho) + + implicit none + include 'constants.include.F' + + integer, intent(in) :: n + complex*16, intent(in) :: rho + + integer :: i, mmax + double precision :: rho_mod, rho_re, rho_im + double precision :: sq_rho_re, sq_rho_im + double precision :: n_tmp + complex*16 :: sq_rho, rho_inv, rho_exp + + complex*16 :: crint_smallz, cpx_erf_1 + complex*16 :: cpx_erf_2 + + rho_re = real (rho) + rho_im = aimag(rho) + rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) + + if(rho_mod < 10.d0) then + ! small z + if(rho_mod .lt. 1.d-15) then + crint_1 = 1.d0 / dble(n + n + 1) + else + crint_1 = crint_smallz(n, rho) + endif + + else + ! large z + + if(rho_mod .gt. 40.d0) then + + n_tmp = dble(n) + 0.5d0 + crint_1 = 0.5d0 * gamma(n_tmp) / (rho**n_tmp) + + else + + ! get \sqrt(rho) + sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) + sq_rho_im = 0.5d0 * rho_im / sq_rho_re + sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im + + rho_exp = 0.5d0 * zexp(-rho) + rho_inv = (1.d0, 0.d0) / rho + + !print*, sq_rho_re, sq_rho_im + !print*, cpx_erf_1(sq_rho_re, sq_rho_im) + !print*, cpx_erf_2(sq_rho_re, sq_rho_im) + + crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + mmax = n + if(mmax .gt. 0) then + do i = 0, mmax-1 + crint_1 = ((dble(i) + 0.5d0) * crint_1 - rho_exp) * rho_inv + enddo + endif + + endif + + endif + +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 + include 'constants.include.F' + + integer, intent(in) :: n_pt_out + complex*16, intent(in) :: rho, d1(0:n_pt_out) + + integer :: n, i, mmax + double precision :: rho_mod, rho_re, rho_im + double precision :: sq_rho_re, sq_rho_im + complex*16 :: sq_rho, F0 + complex*16 :: rho_tmp, rho_inv, rho_exp + complex*16, allocatable :: Fm(:) + + complex*16 :: crint_smallz, cpx_erf_1 + + + rho_re = real (rho) + rho_im = aimag(rho) + rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) + +! ! debug +! double precision :: d1_real(0:n_pt_out) +! double precision :: rint_sum +! do i = 0, n_pt_out +! d1_real(i) = real(d1(i)) +! enddo +! crint_sum_1 = rint_sum(n_pt_out, rho_re, d1_real) +! return + + if(rho_mod < 10.d0) then + ! small z + + if(rho_mod .lt. 1.d-15) then + + crint_sum_1 = d1(0) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_1 = crint_sum_1 + d1(i) / dble(n+n+1) + enddo + + else + + crint_sum_1 = d1(0) * crint_smallz(0, rho) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_1 = crint_sum_1 + d1(i) * crint_smallz(n, rho) + enddo + + endif + + else + ! large z + + if(rho_mod .gt. 40.d0) then + + rho_inv = (1.d0, 0.d0) / rho + rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv) + + crint_sum_1 = rho_tmp * d1(0) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv + crint_sum_1 = crint_sum_1 + rho_tmp * d1(i) + enddo + + else + + ! get \sqrt(rho) + sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) + sq_rho_im = 0.5d0 * rho_im / sq_rho_re + sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im + + F0 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + crint_sum_1 = F0 * d1(0) + + rho_exp = 0.5d0 * zexp(-rho) + rho_inv = (1.d0, 0.d0) / rho + + mmax = shiftr(n_pt_out, 1) + if(mmax .gt. 0) then + + allocate(Fm(mmax)) + Fm(1:mmax) = (0.d0, 0.d0) + + do n = 0, mmax-1 + F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv + Fm(n+1) = F0 + enddo + + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_1 = crint_sum_1 + Fm(n) * d1(i) + enddo + + deallocate(Fm) + endif + + endif ! rho_mod + endif ! rho_mod + +end + +! --- + +complex*16 function crint_smallz(n, rho) + + BEGIN_DOC + ! Standard version of rint + END_DOC + + implicit none + integer, intent(in) :: n + complex*16, intent(in) :: rho + + integer, parameter :: kmax = 40 + double precision, parameter :: eps = 1.d-13 + + integer :: k + double precision :: delta_mod + 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) + + do k = 1, kmax + + rho_k = rho_k * rho + delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0) + crint_smallz = crint_smallz + delta_k + + delta_mod = dsqrt(real(delta_k)*real(delta_k) + aimag(delta_k)*aimag(delta_k)) + if(delta_mod .lt. eps) return + enddo + + if(delta_mod > eps) then + write(*,*) ' pb in crint_smallz !' + write(*,*) ' n, rho = ', n, rho + write(*,*) ' delta_mod = ', delta_mod + stop 1 + endif + +end + +! --- + +complex*16 function crint_2(n, rho) + + implicit none + include 'constants.include.F' + + 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 + + if(abs(rho) < 10.d0) then + + if(abs(rho) .lt. 1d-6) then + tmp = 2.d0 * dble(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) + else + crint_2 = crint_smallz(n, rho) + endif + + else + + if(real(rho) .ge. 0.d0) then + call zboysfun(n, rho, vals) + crint_2 = vals(n) + else + call zboysfunnrp(n, rho, vals) + crint_2 = vals(n) * zexp(-rho) + endif + + endif + + return +end + +! --- + +subroutine zboysfun(n_max, x, vals) + + BEGIN_DOC + ! + ! Computes values of the Boys function for n = 0, 1, ..., n_max + ! for a complex valued argument + ! + ! Input: x --- argument, complex*16, Re(x) >= 0 + ! Output: vals --- values of the Boys function, n = 0, 1, ..., n_max + ! + 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 + 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 + + y = zexp(-x) +! x0_nmax = x0(n_max) + +! 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 + +end + +! --- + +subroutine zboysfunnrp(n_max, x, vals) + + BEGIN_DOC + ! + ! Computes values of e^z F(n,z) for n = 0, 1, ..., n_max + ! (where F(n,z) are the Boys functions) + ! for a complex valued argument WITH NEGATIVE REAL PART + ! + ! Input: x --- argument, complex *16 Re(x)<=0 + ! Output: vals --- values of e^z F(n,z), n = 0, 1, ..., n_max + ! + 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 + complex*16, intent(out) :: vals(0:n_max) + + integer :: n, k + double precision :: x0_nmax + complex*16 :: y, rtmp + complex*16 :: p, q, tmp + + y = exp(x) +! x0_nmax = x0(n_max) + +! 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 + +end + +! --- + +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 :: n, i + + 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 + +end + +! --- + diff --git a/src/utils/cpx_erf.irp.f b/src/utils/cpx_erf.irp.f index 61f81055..1c5fa61d 100644 --- a/src/utils/cpx_erf.irp.f +++ b/src/utils/cpx_erf.irp.f @@ -1,7 +1,7 @@ ! --- -complex*16 function cpx_erf(x, y) +complex*16 function cpx_erf_1(x, y) BEGIN_DOC ! @@ -25,25 +25,25 @@ complex*16 function cpx_erf(x, y) if(yabs .lt. 1.d-15) then - cpx_erf = (1.d0, 0.d0) * derf(x) + cpx_erf_1 = (1.d0, 0.d0) * derf(x) return else erf_tmp1 = (1.d0, 0.d0) * derf(x) erf_tmp2 = erf_E(x, yabs) + erf_F(x, yabs) - erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * ( erf_G(x, yabs) + erf_H(x, yabs) ) + erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * (erf_G(x, yabs) + erf_H(x, yabs)) erf_tot = erf_tmp1 + erf_tmp2 - erf_tmp3 endif if(y .gt. 0.d0) then - cpx_erf = erf_tot + cpx_erf_1 = erf_tot else - cpx_erf = CONJG(erf_tot) + cpx_erf_1 = conjg(erf_tot) endif -end function cpx_erf +end ! --- @@ -54,7 +54,7 @@ complex*16 function erf_E(x, yabs) double precision, intent(in) :: x, yabs - if( (dabs(x).gt.6.d0) .or. (x==0.d0) ) then + if((dabs(x).gt.6.d0) .or. (x==0.d0)) then erf_E = (0.d0, 0.d0) return endif @@ -70,7 +70,7 @@ complex*16 function erf_E(x, yabs) endif -end function erf_E +end ! --- @@ -109,7 +109,7 @@ double precision function erf_F(x, yabs) endif -end function erf_F +end ! --- @@ -149,7 +149,7 @@ complex*16 function erf_G(x, yabs) enddo -end function erf_G +end ! --- @@ -172,7 +172,7 @@ complex*16 function erf_H(x, yabs) endif - if( (dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0) ) then + if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then x2 = x * x ct = 0.5d0 * inv_pi @@ -186,7 +186,7 @@ complex*16 function erf_H(x, yabs) tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1 erf_H = erf_H + tmp2 - tmp_mod = dsqrt(REAL(tmp2)*REAL(tmp2) + AIMAG(tmp2)*AIMAG(tmp2)) + tmp_mod = dsqrt(real(tmp2)*real(tmp2) + aimag(tmp2)*aimag(tmp2)) if(tmp_mod .lt. 1d-15) exit enddo erf_H = ct * erf_H @@ -197,8 +197,394 @@ complex*16 function erf_H(x, yabs) endif -end function erf_H +end ! --- +complex*16 function cpx_erf_2(x, y) + + BEGIN_DOC + ! + ! compute erf(z) for z = x + i y + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, intent(in) :: x, y + + double precision :: yabs + complex*16 :: z + + yabs = dabs(y) + + if(yabs .lt. 1.d-15) then + + cpx_erf_2 = (1.d0, 0.d0) * derf(x) + return + + else + + z = x + (0.d0, 1.d0) * y + + if(x .ge. 0.d0) then + call zboysfun00(z, cpx_erf_2) + else + call zboysfun00nrp(z, cpx_erf_2) + cpx_erf_2 = cpx_erf_2 * zexp(-z) + endif + + endif + + return +end + +! --- + +subroutine zboysfun00(z, val) + + BEGIN_DOC + ! + ! Computes values of the Boys function for n=0 + ! for a complex valued argument + ! + ! Input: z --- argument, complex*16, Real(z) >= 0 + ! Output: val --- value of the Boys function n=0 + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & + 0.249999999999993161d0, & + -0.374999999999766599d0, & + 0.937499999992027020d0, & + -3.28124999972738868d0, & + 14.7656249906697030d0, & + -81.2109371803307752d0 /) + + double precision, parameter :: taylcoef(0:10) = (/ 1.0d0, & + -0.333333333333333333d0, & + 0.1d0, & + -0.238095238095238095d-01, & + 0.462962962962962963d-02, & + -0.757575757575757576d-03, & + 0.106837606837606838d-03, & + -0.132275132275132275d-04, & + 1.458916900093370682d-06, & + -1.450385222315046877d-07, & + 1.3122532963802805073d-08 /) + + double precision, parameter :: sqpio2 = 0.886226925452758014d0 + + double precision, parameter :: pp(1:22) = (/ 0.001477878263796956477d0, & + 0.013317276413725817441d0, & + 0.037063591452052541530d0, & + 0.072752512422882761543d0, & + 0.120236941228785688896d0, & + 0.179574293958937717967d0, & + 0.253534046984087292596d0, & + 0.350388652780721927513d0, & + 0.482109575931276669313d0, & + 0.663028993158374107103d0, & + 0.911814736856590885929d0, & + 1.2539502287919293d0, & + 1.7244634233573395d0, & + 2.3715248262781863d0, & + 3.2613796996078355d0, & + 4.485130169059591d0, & + 6.168062135122484d0, & + 8.48247187231787d0, & + 11.665305486296793d0, & + 16.042417132288328d0, & + 22.06192951814709d0, & + 30.340112094708307d0 /) + + double precision, parameter :: ff(1:22) = (/ 0.0866431027201416556d0, & + 0.0857720608434394764d0, & + 0.0839350436829178814d0, & + 0.0809661970413229146d0, & + 0.0769089548492978618d0, & + 0.0731552078711821626d0, & + 0.0726950035163157228d0, & + 0.0752842556089304050d0, & + 0.0770943953645196145d0, & + 0.0754250625677530441d0, & + 0.0689686192650315305d0, & + 0.05744480422143023d0, & + 0.04208199434694545d0, & + 0.025838539448223282d0, & + 0.012445024157255563d0, & + 0.004292541592599837d0, & + 0.0009354342987735969d0, & + 0.10840885466502504d-03, & + 5.271867966761674d-06, & + 7.765974039750418d-08, & + 2.2138172422680093d-10, & + 6.594161760037707d-14 /) + + complex*16, intent(in) :: z + complex*16, intent(out) :: val + + integer :: k + complex*16 :: z1, zz, y + + zz = zexp(-z) + + if(abs(z) .ge. 100.0d0) then + + ! large |z| + z1 = 1.0d0 / zsqrt(z) + y = 1.0d0 / z + val = asymcoef(7) + do k = 6, 1, -1 + val = val * y + asymcoef(k) + enddo + val = zz * val * y + z1 * sqpio2 + + else if(abs(z) .le. 0.35d0) then + + ! small |z| + val = taylcoef(10) * (1.d0, 0.d0) + do k = 9, 0, -1 + val = val * z + taylcoef(k) + enddo + + else + + ! intermediate |z| + val = sqpio2 / zsqrt(z) - 0.5d0 * zz * sum(ff(1:22)/(z+pp(1:22))) + !val = (0.d0, 0.d0) + !do k = 1, 22 + ! val += ff(k) / (z + pp(k)) + !enddo + !val = sqpio2 / zsqrt(z) - 0.5d0 * zz * val + + endif + + return +end + +! --- + +subroutine zboysfun00nrp(z, val) + + BEGIN_DOC + ! + ! Computes values of the exp(z) F(0,z) + ! (where F(0,z) is the Boys function) + ! for a complex valued argument with Real(z)<=0 + ! + ! Input: z --- argument, complex*16, !!! Real(z)<=0 !!! + ! Output: val --- value of the function !!! exp(z) F(0,z) !!!, where F(0,z) is the Boys function + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & + 0.249999999999993161d0, & + -0.374999999999766599d0, & + 0.937499999992027020d0, & + -3.28124999972738868d0, & + 14.7656249906697030d0, & + -81.2109371803307752d0 /) + + double precision, parameter :: taylcoef(0:10) = (/ 1.0d0, & + -0.333333333333333333d0, & + 0.1d0, & + -0.238095238095238095d-01, & + 0.462962962962962963d-02, & + -0.757575757575757576d-03, & + 0.106837606837606838d-03, & + -0.132275132275132275d-04, & + 1.458916900093370682d-06, & + -1.450385222315046877d-07, & + 1.3122532963802805073d-08 /) + + 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) + + double precision, parameter :: pp(1:16) = (/ 0.005299532504175031d0, & + 0.0277124884633837d0, & + 0.06718439880608407d0, & + 0.12229779582249845d0, & + 0.19106187779867811d0, & + 0.27099161117138637d0, & + 0.35919822461037054d0, & + 0.45249374508118123d0, & + 0.5475062549188188d0, & + 0.6408017753896295d0, & + 0.7290083888286136d0, & + 0.8089381222013219d0, & + 0.8777022041775016d0, & + 0.9328156011939159d0, & + 0.9722875115366163d0, & + 0.994700467495825d0 /) + + double precision, parameter :: ww(1:16) = (/ 0.013576229705876844d0, & + 0.03112676196932382d0, & + 0.04757925584124612d0, & + 0.062314485627766904d0, & + 0.07479799440828848d0, & + 0.08457825969750153d0, & + 0.09130170752246194d0, & + 0.0947253052275344d0, & + 0.0947253052275344d0, & + 0.09130170752246194d0, & + 0.08457825969750153d0, & + 0.07479799440828848d0, & + 0.062314485627766904d0, & + 0.04757925584124612d0, & + 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 :: 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, & + 0.8901093330366746d0, & + 0.6799475005849274d0, & + 0.3900551639790145d0, & + 0.15047608763371934d0, & + 0.0358806749968974d0, & + 0.005089440900100864d0, & + 0.00043900830706867264d0, & + 0.000025161192619824898d0, & + 1.1153308427285078d-6, & + 4.68317018372038d-8, & + 2.3522908467181876d-9, & + 1.7941242138648815d-10, & + 2.5798173021885247d-11, & + 8.27559122014575d-12 /) + + + double precision, parameter :: uu1(1:16) = (/ 0.999247875092057d0, & + 0.979635711599488d0, & + 0.8861006617341018d0, & + 0.6698533710831932d0, & + 0.3760730980014839d0, & + 0.13982165701683388d0, & + 0.031537442321301304d0, & + 0.004147133581658446d0, & + 0.0003253024081883165d0, & + 0.000016687766678889653d0, & + 6.555359391864376d-7, & + 2.4341421258295026d-8, & + 1.0887481200652014d-9, & + 7.51542178140961d-11, & + 1.002378402152542d-11, & + 3.0767730761654096d-12 /) + + complex*16, intent(in) :: z + complex*16, intent(out) :: val + + integer :: k + complex*16 :: z1, zz, y, zsum, tmp, zt, q, p + + zz = zexp(z) + + if(abs(z) .ge. 100.0d0) then + ! large |z| + z1 = 1.0d0 / zsqrt(z) + y = 1.0d0 / z + val = asymcoef(7) + do k = 6, 1, -1 + val = val * y + asymcoef(k) + enddo + val = val * y + z1 * sqpio2 * zz + return + endif + + if(abs(z) .le. 0.35d0) then + ! small |z| + val = taylcoef(10) * (1.d0, 0.d0) + do k = 9, 0, -1 + val = val * z + taylcoef(k) + enddo + val = val * zz + return + endif + + if(abs(etmax+z) .ge. 0.5d0) then + ! intermediate |z| + zsum = (0.d0, 0.d0) + do k = 1, 16 + if(abs(z + qq(k)) .ge. tol) then + 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 + endif + enddo + zt = ima * sqrt(z / etmax) + tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) + val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z) + else + zsum = (0.d0, 0.d0) + do k = 1, 16 + if(abs(z + qq1(k)) .ge. tol) then + 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 + 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) + endif + + return +end + +! ---