LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 03:53:29 +01:00
Code Releases Activity

Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable
Some checks failed
continuous-integration/drone/push Build is failing
Details

Browse Source
This commit is contained in:
Anthony Scemama 2024-10-11 15:43:11 +02:00
commit b7e0ae621a
8 changed files with 1334 additions and 668 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
scripts/import_champ_jastrow.py Executable file → Normal file
View File

@ -45,7 +45,7 @@ if __name__ == '__main__':
jastrow_file = sys.argv[2] jastrow_file = sys.argv[2]
jastrow = import_jastrow(jastrow_file) jastrow = import_jastrow(jastrow_file)
print (jastrow) print (jastrow)
ezfio.set_jastrow_jast_type("Qmckl") ezfio.set_jastrow_j2e_type("Qmckl")
ezfio.set_jastrow_jast_qmckl_type_nucl_num(jastrow['type_num']) ezfio.set_jastrow_jast_qmckl_type_nucl_num(jastrow['type_num'])
charges = ezfio.get_nuclei_nucl_charge() charges = ezfio.get_nuclei_nucl_charge()
types = {} types = {}

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

@ -111,8 +111,9 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a
complex*16 :: accu, P_center(3) complex*16 :: accu, P_center(3)
complex*16 :: d(0:n_pt_in) complex*16 :: d(0:n_pt_in)
complex*16 :: V_n_e_cosgtos complex*16, external :: V_n_e_cosgtos
complex*16 :: crint complex*16, external :: crint_2
complex*16, external :: crint_sum_2
if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. & 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 (A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then
@ -158,27 +159,27 @@ 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)) ) 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 if(n_pt == 0) then
NAI_pol_mult_cosgtos = coeff * crint(0, const) NAI_pol_mult_cosgtos = coeff * crint_2(0, const)
return return
endif endif
call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta & call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, &
, power_A, power_B, C_center, n_pt_in, d, n_pt_out) power_A, power_B, C_center, n_pt_in, d, n_pt_out)
if(n_pt_out < 0) then if(n_pt_out < 0) then
NAI_pol_mult_cosgtos = (0.d0, 0.d0) NAI_pol_mult_cosgtos = (0.d0, 0.d0)
return return
endif endif
accu = (0.d0, 0.d0) !accu = (0.d0, 0.d0)
do i = 0, n_pt_out, 2 !do i = 0, n_pt_out, 2
accu += crint(shiftr(i, 1), const) * d(i) ! accu += crint_2(shiftr(i, 1), const) * d(i)
!enddo
! print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const)) accu = crint_sum_2(n_pt_out, const, d)
enddo
NAI_pol_mult_cosgtos = accu * coeff NAI_pol_mult_cosgtos = accu * coeff
end function NAI_pol_mult_cosgtos return
end
! --- ! ---
@ -312,7 +313,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
d(i) = d1(i) d(i) = d1(i)
enddo enddo
end subroutine give_cpolynomial_mult_center_one_e end
! --- ! ---
@ -405,7 +406,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
endif endif
end subroutine I_x1_pol_mult_one_e_cosgtos end
! --- ! ---
@ -467,7 +468,7 @@ recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
endif endif
end subroutine I_x2_pol_mult_one_e_cosgtos end
! --- ! ---
@ -502,7 +503,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) * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
endif endif
end function V_n_e_cosgtos end
! --- ! ---
@ -529,7 +530,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) V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
endif endif
end function V_r_cosgtos end
! --- ! ---

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

@ -35,11 +35,9 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then
!print *, ' with shwartz acc '
ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l) ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l)
else else
!print *, ' without shwartz acc '
dim1 = n_pt_max_integrals dim1 = n_pt_max_integrals
@ -51,7 +49,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
ao_two_e_integral_cosgtos = 0.d0 ao_two_e_integral_cosgtos = 0.d0
if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then
!print *, ' not the same center'
do p = 1, 3 do p = 1, 3
I_power(p) = ao_power(i,p) I_power(p) = ao_power(i,p)
@ -72,72 +69,22 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j)
expo2 = ao_expo_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 & 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 ) expo1, expo2, I_power, J_power, I_center, J_center, dim1)
p1_inv = (1.d0,0.d0) / pp1 p1_inv = (1.d0,0.d0) / pp1
call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & 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 ) conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1)
p2_inv = (1.d0,0.d0) / pp2 p2_inv = (1.d0,0.d0) / pp2
call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & 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 ) expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1)
p3_inv = (1.d0,0.d0) / pp3 p3_inv = (1.d0,0.d0) / pp3
call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & 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 ) conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1)
p4_inv = (1.d0,0.d0) / pp4 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) do r = 1, ao_prim_num(k)
coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k) coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k)
expo3 = ao_expo_ord_transp_cosgtos(r,k) expo3 = ao_expo_ord_transp_cosgtos(r,k)
@ -146,66 +93,40 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l)
expo4 = ao_expo_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 & 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 ) expo3, expo4, K_power, L_power, K_center, L_center, dim1)
q1_inv = (1.d0,0.d0) / qq1 q1_inv = (1.d0,0.d0) / qq1
call give_explicit_cpoly_and_cgaussian( Q2_new, Q2_center, qq2, fact_q2, iorder_q2 & 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 ) conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1)
q2_inv = (1.d0,0.d0) / qq2 q2_inv = (1.d0,0.d0) / qq2
!do ii = 1, 3 integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
! !print *, qq1, q1_inv Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
! !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
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 & 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 ) 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 & 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 ) 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 & 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 ) 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 & 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 ) 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 & 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 ) 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 & 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 ) 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 )
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
!integral_tot = integral1
!print*, integral_tot
ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot) ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot)
enddo ! s enddo ! s
enddo ! r enddo ! r
@ -213,7 +134,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
enddo ! p enddo ! p
else else
!print *, ' the same center'
do p = 1, 3 do p = 1, 3
I_power(p) = ao_power(i,p) I_power(p) = ao_power(i,p)
@ -290,7 +210,7 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
endif endif
endif endif
end function ao_two_e_integral_cosgtos end
! --- ! ---
@ -326,8 +246,8 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
double precision :: thr double precision :: thr
double precision :: schwartz_ij double precision :: schwartz_ij
complex*16 :: ERI_cosgtos complex*16, external :: ERI_cosgtos
complex*16 :: general_primitive_integral_cosgtos complex*16, external :: general_primitive_integral_cosgtos
ao_2e_cosgtos_schwartz_accel = 0.d0 ao_2e_cosgtos_schwartz_accel = 0.d0
@ -366,45 +286,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l)
expo2 = ao_expo_ord_transp_cosgtos(s,l) expo2 = ao_expo_ord_transp_cosgtos(s,l)
call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 & call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, &
, expo1, expo2, K_power, L_power, K_center, L_center, dim1 ) expo1, expo2, K_power, L_power, K_center, L_center, dim1)
p1_inv = (1.d0,0.d0) / pp1 p1_inv = (1.d0,0.d0) / pp1
call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, &
, conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1 ) conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1)
p2_inv = (1.d0,0.d0) / pp2 p2_inv = (1.d0,0.d0) / pp2
call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, &
, expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1 ) expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1)
p3_inv = (1.d0,0.d0) / pp3 p3_inv = (1.d0,0.d0) / pp3
call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, &
, conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1 ) conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1)
p4_inv = (1.d0,0.d0) / pp4 p4_inv = (1.d0,0.d0) / pp4
integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -544,41 +464,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l)
expo2 = ao_expo_ord_transp_cosgtos(s,l) expo2 = ao_expo_ord_transp_cosgtos(s,l)
integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, &
, K_power(1), L_power(1), K_power(1), L_power(1) & K_power(1), L_power(1), K_power(1), L_power(1), &
, K_power(2), L_power(2), K_power(2), L_power(2) & K_power(2), L_power(2), K_power(2), L_power(2), &
, K_power(3), L_power(3), K_power(3), L_power(3) ) K_power(3), L_power(3), K_power(3), L_power(3))
integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 &
, K_power(1), L_power(1), K_power(1), L_power(1) &
, K_power(2), L_power(2), K_power(2), L_power(2) &
, K_power(3), L_power(3), K_power(3), L_power(3) )
integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, &
, K_power(1), L_power(1), K_power(1), L_power(1) & K_power(1), L_power(1), K_power(1), L_power(1), &
, K_power(2), L_power(2), K_power(2), L_power(2) & K_power(2), L_power(2), K_power(2), L_power(2), &
, K_power(3), L_power(3), K_power(3), L_power(3) ) K_power(3), L_power(3), K_power(3), L_power(3))
integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 &
, K_power(1), L_power(1), K_power(1), L_power(1) &
, K_power(2), L_power(2), K_power(2), L_power(2) &
, K_power(3), L_power(3), K_power(3), L_power(3) )
integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, &
, K_power(1), L_power(1), K_power(1), L_power(1) & K_power(1), L_power(1), K_power(1), L_power(1), &
, K_power(2), L_power(2), K_power(2), L_power(2) & K_power(2), L_power(2), K_power(2), L_power(2), &
, K_power(3), L_power(3), K_power(3), L_power(3) ) K_power(3), L_power(3), K_power(3), L_power(3))
integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 &
, K_power(1), L_power(1), K_power(1), L_power(1) &
, K_power(2), L_power(2), K_power(2), L_power(2) &
, K_power(3), L_power(3), K_power(3), L_power(3) )
integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, &
, K_power(1), L_power(1), K_power(1), L_power(1) & K_power(1), L_power(1), K_power(1), L_power(1), &
, K_power(2), L_power(2), K_power(2), L_power(2) & K_power(2), L_power(2), K_power(2), L_power(2), &
, K_power(3), L_power(3), K_power(3), L_power(3) ) K_power(3), L_power(3), K_power(3), L_power(3))
integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 &
, K_power(1), L_power(1), K_power(1), L_power(1) & integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, &
, K_power(2), L_power(2), K_power(2), L_power(2) & K_power(1), L_power(1), K_power(1), L_power(1), &
, K_power(3), L_power(3), K_power(3), L_power(3) ) K_power(2), L_power(2), K_power(2), L_power(2), &
K_power(3), L_power(3), K_power(3), L_power(3))
integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, &
K_power(1), L_power(1), K_power(1), L_power(1), &
K_power(2), L_power(2), K_power(2), L_power(2), &
K_power(3), L_power(3), K_power(3), L_power(3))
integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, &
K_power(1), L_power(1), K_power(1), L_power(1), &
K_power(2), L_power(2), K_power(2), L_power(2), &
K_power(3), L_power(3), K_power(3), L_power(3))
integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, &
K_power(1), L_power(1), K_power(1), L_power(1), &
K_power(2), L_power(2), K_power(2), L_power(2), &
K_power(3), L_power(3), K_power(3), L_power(3))
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -598,45 +522,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j)
expo2 = ao_expo_ord_transp_cosgtos(q,j) expo2 = ao_expo_ord_transp_cosgtos(q,j)
integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, &
, I_power(1), J_power(1), I_power(1), J_power(1) & I_power(1), J_power(1), I_power(1), J_power(1), &
, I_power(2), J_power(2), I_power(2), J_power(2) & I_power(2), J_power(2), I_power(2), J_power(2), &
, I_power(3), J_power(3), I_power(3), J_power(3) ) I_power(3), J_power(3), I_power(3), J_power(3))
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -655,45 +579,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l)
expo4 = ao_expo_ord_transp_cosgtos(s,l) expo4 = ao_expo_ord_transp_cosgtos(s,l)
integral1 = ERI_cosgtos( expo1, expo2, expo3, expo4 & integral1 = ERI_cosgtos(expo1, expo2, expo3, expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral2 = ERI_cosgtos( expo1, expo2, conjg(expo3), expo4 & integral2 = ERI_cosgtos(expo1, expo2, conjg(expo3), expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral3 = ERI_cosgtos( conjg(expo1), expo2, expo3, expo4 & integral3 = ERI_cosgtos(conjg(expo1), expo2, expo3, expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo3), expo4 & integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo3), expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral5 = ERI_cosgtos( expo1, conjg(expo2), expo3, expo4 & integral5 = ERI_cosgtos(expo1, conjg(expo2), expo3, expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo3), expo4 & integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo3), expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo3, expo4 & integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo3, expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo3), expo4 & integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, &
, I_power(1), J_power(1), K_power(1), L_power(1) & I_power(1), J_power(1), K_power(1), L_power(1), &
, I_power(2), J_power(2), K_power(2), L_power(2) & I_power(2), J_power(2), K_power(2), L_power(2), &
, I_power(3), J_power(3), K_power(3), L_power(3) ) I_power(3), J_power(3), K_power(3), L_power(3))
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -707,7 +631,7 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
deallocate(schwartz_kl) deallocate(schwartz_kl)
end function ao_2e_cosgtos_schwartz_accel end
! --- ! ---
@ -739,8 +663,8 @@ END_PROVIDER
! --- ! ---
complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fact_p, p, p_inv, iorder_p & 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 ) Q_new, Q_center, fact_q, q, q_inv, iorder_q)
BEGIN_DOC BEGIN_DOC
! !
@ -765,7 +689,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 :: 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 :: 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 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol
@ -912,13 +836,11 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) 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 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
! --- ! ---
@ -994,7 +916,7 @@ complex*16 function ERI_cosgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_
ERI_cosgtos = I_f * coeff ERI_cosgtos = I_f * coeff
end function ERI_cosgtos end
! --- ! ---
@ -1076,7 +998,7 @@ subroutine integrale_new_cosgtos(I_f, a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_
I_f += gauleg_w(i, j) * t1(i) I_f += gauleg_w(i, j) * t1(i)
enddo enddo
end subroutine integrale_new_cosgtos end
! --- ! ---
@ -1123,7 +1045,7 @@ recursive subroutine I_x1_new_cosgtos(a, c, B_10, B_01, B_00, res, n_pt)
endif endif
end subroutine I_x1_new_cosgtos end
! --- ! ---
@ -1163,7 +1085,7 @@ recursive subroutine I_x2_new_cosgtos(c, B_10, B_01, B_00, res, n_pt)
endif endif
end subroutine I_x2_new_cosgtos end
! --- ! ---
@ -1234,7 +1156,7 @@ subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt
return return
endif endif
end subroutine give_cpolynom_mult_center_x end
! --- ! ---
@ -1277,7 +1199,7 @@ subroutine I_x1_pol_mult_cosgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt
endif endif
end subroutine I_x1_pol_mult_cosgtos end
! --- ! ---
@ -1370,7 +1292,7 @@ recursive subroutine I_x1_pol_mult_recurs_cosgtos(a, c, B_10, B_01, B_00, C_00,
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call multiply_cpoly(Y, ny, C_00, 2, d, nd) call multiply_cpoly(Y, ny, C_00, 2, d, nd)
end subroutine I_x1_pol_mult_recurs_cosgtos end
! --- ! ---
@ -1426,7 +1348,7 @@ recursive subroutine I_x1_pol_mult_a1_cosgtos(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call multiply_cpoly(Y, ny, C_00, 2, d, nd) call multiply_cpoly(Y, ny, C_00, 2, d, nd)
end subroutine I_x1_pol_mult_a1_cosgtos end
! --- ! ---
@ -1490,7 +1412,7 @@ recursive subroutine I_x1_pol_mult_a2_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call multiply_cpoly(Y, ny, C_00, 2, d, nd) call multiply_cpoly(Y, ny, C_00, 2, d, nd)
end subroutine I_x1_pol_mult_a2_cosgtos end
! --- ! ---
@ -1575,7 +1497,7 @@ recursive subroutine I_x2_pol_mult_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d, n
end select end select
end subroutine I_x2_pol_mult_cosgtos end
! --- ! ---

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) call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new)
iorder = a + b 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 !DIR$ FORCEINLINE
call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new) 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(2) = ( a * xa(2) + b * xb(2) ) * p_inv
xp(3) = ( a * xa(3) + b * xb(3) ) * 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) k = zexp(-k)
xp = (a*xa + b*xb) * p_inv 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 exit
enddo enddo
end subroutine multiply_cpoly end
! --- ! ---
@ -373,7 +373,7 @@ subroutine add_cpoly(b, nb, c, nc, d, nd)
if(nd < 0) exit if(nd < 0) exit
enddo enddo
end subroutine add_cpoly end
! --- ! ---
@ -413,7 +413,7 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd)
endif 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) P_B(i) = binom_func(b,b-i) * pows_b(b-i)
enddo 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)) Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1))
end function Fc_integral end
! ---
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
! --- ! ---

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

@ -0,0 +1,543 @@
! ---
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)
crint_smallz = ct / gamma(dble(n) + 1.5d0)
rho_k = (1.d0, 0.d0)
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
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
integer, intent(in) :: n_max
complex*16, intent(in) :: x
complex*16, intent(out) :: vals(0:n_max)
integer :: n
complex*16 :: yy, x_inv
call zboysfun00(x, vals(0))
yy = 0.5d0 * zexp(-x)
x_inv = (1.d0, 0.d0) / x
do n = 1, n_max
vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - yy) * x_inv
enddo
return
end
! ---
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
integer, intent(in) :: n_max
complex*16, intent(in) :: x
complex*16, intent(out) :: vals(0:n_max)
integer :: n
complex*16 :: x_inv
call zboysfun00nrp(x, vals(0))
x_inv = (1.d0, 0.d0) / x
do n = 1, n_max
vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - 0.5d0) * x_inv
enddo
return
end
! ---
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
integer :: n_max
complex*16, allocatable :: vals(:)
!complex*16, external :: crint_2
!crint_sum_2 = (0.d0, 0.d0)
!do i = 0, n_pt_out, 2
! n = shiftr(i, 1)
! crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho)
!enddo
n_max = shiftr(n_pt_out, 1)
allocate(vals(0:n_max))
call crint_2_vec(n_max, rho, vals)
crint_sum_2 = d1(0) * vals(0)
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum_2 += d1(i) * vals(n)
enddo
deallocate(vals)
return
end
! ---
subroutine crint_2_vec(n_max, rho, vals)
implicit none
integer, intent(in) :: n_max
complex*16, intent(in) :: rho
complex*16, intent(out) :: vals(0:n_max)
integer :: n
double precision :: tmp, abs_rho
complex*16 :: rho2, rho3, erho
abs_rho = abs(rho)
if(abs_rho < 10.d0) then
if(abs_rho .lt. 1d-6) then
! use finite expansion for very small rho
! rho^2 / 2
rho2 = 0.5d0 * rho * rho
! rho^3 / 6
rho3 = 0.3333333333333333d0 * rho * rho2
vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3
do n = 1, n_max
tmp = 2.d0 * dble(n)
vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) &
+ rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0)
enddo
else
call crint_smallz_vec(n_max, rho, vals)
endif
else
if(real(rho) .ge. 0.d0) then
call zboysfun(n_max, rho, vals)
else
call zboysfunnrp(n_max, rho, vals)
erho = zexp(-rho)
do n = 0, n_max
vals(n) = vals(n) * erho
enddo
endif
endif
return
end
! ---
subroutine crint_smallz_vec(n_max, rho, vals)
BEGIN_DOC
! Standard version of rint
END_DOC
implicit none
integer, intent(in) :: n_max
complex*16, intent(in) :: rho
complex*16, intent(out) :: vals(0:n_max)
integer, parameter :: kmax = 40
double precision, parameter :: eps = 1.d-13
integer :: k, n
complex*16 :: ct, delta_k
complex*16 :: rhoe
complex*16, allocatable :: rho_k(:)
allocate(rho_k(0:kmax))
rho_k(0) = (1.d0, 0.d0)
do k = 1, kmax
rho_k(k) = rho_k(k-1) * rho
enddo
rhoe = 0.5d0 * zexp(-rho)
do n = 0, n_max
ct = rhoe * gamma(dble(n) + 0.5d0)
vals(n) = ct / gamma(dble(n) + 1.5d0)
do k = 1, kmax
delta_k = ct * rho_k(k) / gamma(dble(n+k) + 1.5d0)
vals(n) += delta_k
if(abs(delta_k) .lt. eps) then
exit
endif
enddo
!if(abs(delta_k) > eps) then
! write(*,*) ' pb in crint_smallz_vec !'
! write(*,*) ' n, rho = ', n, rho
! write(*,*) ' |delta_k| = ', abs(delta_k)
!endif
enddo
deallocate(rho_k)
return
end
! ---

406
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 BEGIN_DOC
! !
@ -25,7 +25,7 @@ complex*16 function cpx_erf(x, y)
if(yabs .lt. 1.d-15) then 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 return
else else
@ -38,12 +38,12 @@ complex*16 function cpx_erf(x, y)
endif endif
if(y .gt. 0.d0) then if(y .gt. 0.d0) then
cpx_erf = erf_tot cpx_erf_1 = erf_tot
else else
cpx_erf = CONJG(erf_tot) cpx_erf_1 = conjg(erf_tot)
endif endif
end function cpx_erf end
! --- ! ---
@ -70,7 +70,7 @@ complex*16 function erf_E(x, yabs)
endif endif
end function erf_E end
! --- ! ---
@ -109,7 +109,7 @@ double precision function erf_F(x, yabs)
endif endif
end function erf_F end
! --- ! ---
@ -149,7 +149,7 @@ complex*16 function erf_G(x, yabs)
enddo enddo
end function erf_G end
! --- ! ---
@ -186,7 +186,7 @@ complex*16 function erf_H(x, yabs)
tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1 tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1
erf_H = erf_H + tmp2 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 if(tmp_mod .lt. 1d-15) exit
enddo enddo
erf_H = ct * erf_H erf_H = ct * erf_H
@ -197,8 +197,394 @@ complex*16 function erf_H(x, yabs)
endif 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
! ---