PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 12:23:43 +01:00
Code Releases Activity

Fixed cplx Boys fct for real arguments

Browse Source
This commit is contained in:
Abdallah Ammar 2024-10-09 01:07:57 +02:00
parent 1f1321c962
commit 7f57051739
8 changed files with 1335 additions and 515 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f
View File

@ -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
! ---

144
src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f
View File

@ -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,63 +96,45 @@ 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
@ -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
! ---

2
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -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

188
src/hartree_fock/deb_ao_2e_int.irp.f Normal file
View File

@ -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
! ---

114
src/hartree_fock/print_scf_int.irp.f
View File

@ -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

278
src/utils/cgtos_utils.irp.f
View File

@ -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
! ---

686
src/utils/cpx_boys.irp.f Normal file
View File

@ -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
! ---

412
src/utils/cpx_erf.irp.f
View File

@ -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
! ---