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

few opt in cGTOs

This commit is contained in:
AbdAmmar 2024-10-12 16:44:41 +02:00
parent ac171d49f6
commit 23a242052d
8 changed files with 792 additions and 490 deletions

View File

@ -3,7 +3,7 @@ logical function ao_two_e_integral_zero(i,j,k,l)
integer, intent(in) :: i,j,k,l
ao_two_e_integral_zero = .False.
if (.not.(read_ao_two_e_integrals.or.is_periodic)) then
if (.not.(read_ao_two_e_integrals.or.is_periodic.or.use_cosgtos)) then
if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then
ao_two_e_integral_zero = .True.
return

View File

@ -15,23 +15,21 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
integer :: p, q, r, s
integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3)
integer :: iorder_p1(3), iorder_p2(3), iorder_p3(3), iorder_p4(3), iorder_q1(3), iorder_q2(3)
integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3)
double precision :: coef1, coef2, coef3, coef4
complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3)
complex*16 :: expo1, expo2, expo3, expo4
complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv
complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv
complex*16 :: P3_new(0:max_dim,3), P3_center(3), fact_p3, pp3, p3_inv
complex*16 :: P4_new(0:max_dim,3), P4_center(3), fact_p4, pp4, p4_inv
complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv
complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv
complex*16 :: integral1, integral2, integral3, integral4
complex*16 :: integral5, integral6, integral7, integral8
complex*16 :: integral_tot
double precision :: ao_2e_cosgtos_schwartz_accel
complex*16 :: ERI_cosgtos
complex*16 :: general_primitive_integral_cosgtos
double precision, external :: ao_2e_cosgtos_schwartz_accel
complex*16, external :: ERI_cosgtos
complex*16, external :: general_primitive_integral_cosgtos
if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then
@ -77,14 +75,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1)
p2_inv = (1.d0, 0.d0) / pp2
call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, &
expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1)
p3_inv = (1.d0,0.d0) / pp3
call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, &
conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1)
p4_inv = (1.d0,0.d0) / pp4
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)
@ -101,28 +91,36 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1)
q2_inv = (1.d0, 0.d0) / qq2
integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
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, &
integral2 = general_primitive_integral_cosgtos(dim1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
integral3 = general_primitive_integral_cosgtos(dim1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
integral4 = general_primitive_integral_cosgtos(dim1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
integral5 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
integral6 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
integral7 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
integral8 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -158,45 +156,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)
integral1 = ERI_cosgtos( expo1, expo2, expo3, expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral1 = ERI_cosgtos(expo1, expo2, expo3, expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral2 = ERI_cosgtos( expo1, expo2, conjg(expo3), expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral2 = ERI_cosgtos(expo1, expo2, conjg(expo3), expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral3 = ERI_cosgtos( conjg(expo1), expo2, expo3, expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral3 = ERI_cosgtos(conjg(expo1), expo2, expo3, expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo3), expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo3), expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral5 = ERI_cosgtos( expo1, conjg(expo2), expo3, expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral5 = ERI_cosgtos(expo1, conjg(expo2), expo3, expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo3), expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo3), expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo3, expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo3, expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo3), expo4 &
, I_power(1), J_power(1), K_power(1), L_power(1) &
, I_power(2), J_power(2), K_power(2), L_power(2) &
, I_power(3), J_power(3), K_power(3), L_power(3) )
integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), &
I_power(2), J_power(2), K_power(2), L_power(2), &
I_power(3), J_power(3), K_power(3), L_power(3))
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -207,8 +205,8 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
enddo ! q
enddo ! p
endif
endif
endif ! same centers
endif ! do schwartz
end
@ -228,14 +226,12 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
integer :: p, q, r, s
integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3)
integer :: iorder_p1(3), iorder_p2(3), iorder_p3(3), iorder_p4(3), iorder_q1(3), iorder_q2(3)
integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3)
double precision :: coef1, coef2, coef3, coef4
complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3)
complex*16 :: expo1, expo2, expo3, expo4
complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv
complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv
complex*16 :: P3_new(0:max_dim,3), P3_center(3), fact_p3, pp3, p3_inv
complex*16 :: P4_new(0:max_dim,3), P4_center(3), fact_p4, pp4, p4_inv
complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv
complex*16 :: Q2_new(0:max_dim,3), Q2_center(3), fact_q2, qq2, q2_inv
complex*16 :: integral1, integral2, integral3, integral4
@ -294,41 +290,40 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1)
p2_inv = (1.d0, 0.d0) / pp2
call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, &
expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1)
p3_inv = (1.d0,0.d0) / pp3
call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, &
conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1)
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)
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)
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)
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)
integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
integral5 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
integral6 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
integral7 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
integral8 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
schwartz_kl(s,r) = coef2 * 2.d0 * real(integral_tot)
schwartz_kl(0,r) = max(schwartz_kl(0,r), schwartz_kl(s,r))
@ -346,45 +341,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j)
expo2 = ao_expo_ord_transp_cosgtos(q,j)
call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 &
, expo1, expo2, I_power, J_power, I_center, J_center, dim1 )
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 )
p3_inv = (1.d0,0.d0) / pp3
integral1 = general_primitive_integral_cosgtos(dim1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 &
, conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1 )
p4_inv = (1.d0,0.d0) / pp4
integral2 = general_primitive_integral_cosgtos(dim1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
integral3 = general_primitive_integral_cosgtos(dim1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
integral4 = general_primitive_integral_cosgtos(dim1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
integral5 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
integral6 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
integral7 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
integral8 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -404,38 +399,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l)
expo4 = ao_expo_ord_transp_cosgtos(s,l)
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
integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
, Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
integral1 = general_primitive_integral_cosgtos(dim1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
, Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
integral2 = general_primitive_integral_cosgtos(dim1, &
P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
, Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
integral3 = general_primitive_integral_cosgtos(dim1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
, Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
integral4 = general_primitive_integral_cosgtos(dim1, &
P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral5 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
, Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
integral6 = general_primitive_integral_cosgtos(dim1, &
conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
, Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
integral7 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
, Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
, Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
integral8 = general_primitive_integral_cosgtos(dim1, &
conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, &
Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
@ -705,14 +707,6 @@ complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fac
p10_2 = pq_inv_2 * p10_1 * q ! 0.5d0*q/(pq + p*p)
p01_2 = pq_inv_2 * p01_1 * p ! 0.5d0*p/(q*q + pq)
! get \sqrt(p + q)
!ppq = p + q
!ppq_re = REAL (ppq)
!ppq_im = AIMAG(ppq)
!ppq_mod = dsqrt(ppq_re*ppq_re + ppq_im*ppq_im)
!sq_ppq_re = sq_op5 * dsqrt(ppq_re + ppq_mod)
!sq_ppq_im = 0.5d0 * ppq_im / sq_ppq_re
!sq_ppq = sq_ppq_re + (0.d0, 1.d0) * sq_ppq_im
sq_ppq = zsqrt(p + q)
! ---
@ -727,13 +721,13 @@ complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fac
do i = 0, iorder_p(1)
tmp_p = P_new(i,1)
tmp_mod = dsqrt(REAL(tmp_p)*REAL(tmp_p) + AIMAG(tmp_p)*AIMAG(tmp_p))
tmp_mod = dsqrt(real(tmp_p)*real(tmp_p) + aimag(tmp_p)*aimag(tmp_p))
if(tmp_mod < thresh) cycle
do j = 0, iorder_q(1)
tmp_q = tmp_p * Q_new(j,1)
tmp_mod = dsqrt(REAL(tmp_q)*REAL(tmp_q) + AIMAG(tmp_q)*AIMAG(tmp_q))
tmp_mod = dsqrt(real(tmp_q)*real(tmp_q) + aimag(tmp_q)*aimag(tmp_q))
if(tmp_mod < thresh) cycle
!DIR$ FORCEINLINE
@ -1089,8 +1083,8 @@ end
! ---
subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt_in &
, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, d, n_pt_out)
subroutine give_cpolynom_mult_center_x(P_center, Q_center, a_x, d_x, p, q, n_pt_in, &
pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, d, n_pt_out)
BEGIN_DOC
! subroutine that returns the explicit polynom in term of the "t"
@ -1142,12 +1136,6 @@ subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt
call I_x1_pol_mult_cosgtos(a_x, d_x, B10, B01, B00, C00, D00, d, n_pt1, n_pt_in)
n_pt_out = n_pt1
! print *, ' '
! print *, a_x, d_x
! print *, real(B10), real(B01), real(B00), real(C00), real(D00)
! print *, n_pt1, real(d(0:n_pt1))
! print *, ' '
if(n_pt1 < 0) then
n_pt_out = -1
do i = 0, n_pt_in
@ -1501,4 +1489,3 @@ end
! ---

View File

@ -44,7 +44,6 @@ double precision function ao_two_e_integral(i, j, k, l)
logical, external :: do_schwartz_accel
if(use_cosgtos) then
!print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos
ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l)

View File

@ -1,9 +1,12 @@
program deb_ao_2e_int
!call check_ao_two_e_integral_cosgtos()
call check_crint1()
implicit none
call check_ao_two_e_integral_cosgtos()
!call check_crint1()
!call check_crint2()
!call check_crint3()
end
@ -14,8 +17,8 @@ subroutine check_ao_two_e_integral_cosgtos()
implicit none
integer :: i, j, k, l
double precision :: tmp1, tmp2
double precision :: acc, nrm, dif
double precision :: tmp1, tmp2
double precision, external :: ao_two_e_integral
double precision, external :: ao_two_e_integral_cosgtos
@ -23,10 +26,10 @@ subroutine check_ao_two_e_integral_cosgtos()
acc = 0.d0
nrm = 0.d0
i = 1
j = 6
k = 1
l = 16
i = 11
j = 100
k = 74
l = 104
! do i = 1, ao_num
! do k = 1, ao_num
! do j = 1, ao_num
@ -35,15 +38,15 @@ subroutine check_ao_two_e_integral_cosgtos()
tmp1 = ao_two_e_integral (i, j, k, l)
tmp2 = ao_two_e_integral_cosgtos(i, j, k, l)
dif = dabs(tmp1 - tmp2)
if(dif .gt. 1d-12) then
dif = abs(tmp1 - tmp2)
!if(dif .gt. 1d-10) then
print*, ' error on:', i, j, k, l
print*, tmp1, tmp2, dif
stop
endif
!stop
!endif
! acc += dif
! nrm += dabs(tmp1)
acc += dif
nrm += abs(tmp1)
! enddo
! enddo
! enddo
@ -73,8 +76,9 @@ subroutine check_crint1()
(+1d+4, +1d+4) /)
complex*16 :: rho
complex*16 :: int_an, int_nm
double precision, external :: rint
complex*16, external :: crint_1, crint_2, crint_quad
complex*16, external :: crint_1, crint_2
n = 10
dif_thr = 1d-7
@ -94,7 +98,7 @@ subroutine check_crint1()
do i = 0, n
!int_an = crint_1(i, rho)
int_an = crint_2(i, rho)
int_nm = crint_quad(i, rho)
call crint_quad_1(i, rho, 100000000, int_nm)
dif_re = dabs(real(int_an) - real(int_nm))
dif_im = dabs(aimag(int_an) - aimag(int_nm))
@ -185,4 +189,144 @@ end
! ---
subroutine check_crint3()
implicit none
integer :: i_test, n_test
integer :: nx, ny, n, n_quad
integer :: i, seed_size, clock_time
double precision :: xr(1:4), x
double precision :: yr(1:4), y
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
double precision :: delta_ref
double precision :: t1, t2, t_int1, t_int2
complex*16 :: rho
complex*16 :: int1_old, int1_ref, int2_old, int2_ref
integer, allocatable :: seed(:)
complex*16, external :: crint_2
call random_seed(size=seed_size)
allocate(seed(seed_size))
call system_clock(count=clock_time)
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
!seed = 123456789
call random_seed(put=seed)
t_int1 = 0.d0
t_int2 = 0.d0
n_test = 5
acc_re = 0.d0
nrm_re = 0.d0
acc_im = 0.d0
nrm_im = 0.d0
do i_test = 1, n_test
! Re(rho)
call random_number(xr)
x = xr(1)
if(xr(2) .gt. 0.5d0) x = -x
nx = int(15.d0 * xr(3))
if(xr(4) .gt. 0.5d0) nx = -nx
x = x * 10.d0**nx
! Im(rho)
call random_number(yr)
y = yr(1)
if(yr(2) .gt. 0.5d0) y = -y
ny = int(5.d0 * yr(3))
if(yr(4) .gt. 0.5d0) ny = -ny
y = y * 10.d0**ny
rho = x + (0.d0, 1.d0) * y
call random_number(x)
x = 31.d0 * x
n = int(x)
!if(n.eq.0) cycle
print*, " n = ", n
print*, " rho = ", real(rho), aimag(rho)
call wall_time(t1)
int1_old = crint_2(n, rho)
!n_quad = 10000000
!call crint_quad_1(n, rho, n_quad, int1_old)
!!delta_ref = 1.d0
!!do while(delta_ref .gt. 1d-12)
!! n_quad = n_quad * 10
!! !print*, " delta = ", delta_ref
!! !print*, " increasing n_quad to:", n_quad
!! call crint_quad_1(n, rho, n_quad, int1_ref)
!! delta_ref = abs(int1_ref - int1_old)
!! int1_old = int1_ref
!! if(n_quad .ge. 1000000000) then
!! print*, ' convergence was not reached for crint_quad_1'
!! print*, " delta = ", delta_ref
!! exit
!! endif
!!enddo
call wall_time(t2)
t_int1 = t_int1 + t2 - t1
!print*, " n_quad for crint_quad_1:", n_quad
call wall_time(t1)
n_quad = 10000000
call crint_quad_12(n, rho, n_quad, int2_old)
!delta_ref = 1.d0
!do while(delta_ref .gt. 1d-12)
! n_quad = n_quad * 10
! !print*, " delta = ", delta_ref
! !print*, " increasing n_quad to:", n_quad
! call crint_quad_12(n, rho, n_quad, int2_ref)
! delta_ref = abs(int2_ref - int2_old)
! int2_old = int2_ref
! if(n_quad .ge. 1000000000) then
! print*, ' convergence was not reached for crint_quad_2'
! print*, " delta = ", delta_ref
! exit
! endif
!enddo
call wall_time(t2)
t_int2 = t_int2 + t2 - t1
!print*, " n_quad for crint_quad_2:", n_quad
dif_re = dabs(real(int1_old) - real(int2_old))
dif_im = dabs(aimag(int1_old) - aimag(int2_old))
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
print*, ' important error found: '
print*, " n = ", n
print*, " rho = ", real(rho), aimag(rho)
print*, real(int1_old), real(int2_old), dif_re
print*, aimag(int1_old), aimag(int2_old), dif_im
!stop
endif
if((real(int1_old) /= real(int1_old)) .or. (aimag(int1_old) /= aimag(int1_old)) .or. &
(real(int2_old) /= real(int2_old)) .or. (aimag(int2_old) /= aimag(int2_old)) ) then
cycle
else
acc_re += dif_re
acc_im += dif_im
nrm_re += dabs(real(int1_old))
nrm_im += dabs(aimag(int1_old))
endif
enddo
print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re + 1d-15)
print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im + 1d-15)
print*, "crint_quad_1 wall time (sec) = ", t_int1
print*, "crint_quad_2 wall time (sec) = ", t_int2
deallocate(seed)
end
! ---

View File

@ -145,68 +145,6 @@ end
! ---
!subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
! BEGIN_DOC
! ! Transforms the product of
! ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3)
! ! exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) exp(-(r-Nucl_center)^2 gama
! !
! ! into
! ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
! ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
! ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
! END_DOC
! implicit none
! include 'constants.include.F'
! integer, intent(in) :: dim
! integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
! double precision, intent(in) :: alpha, beta, gama ! exponents
! double precision, intent(in) :: A_center(3) ! A center
! double precision, intent(in) :: B_center (3) ! B center
! double precision, intent(in) :: Nucl_center(3) ! B center
! double precision, intent(out) :: P_center(3) ! new center
! double precision, intent(out) :: p ! new exponent
! double precision, intent(out) :: fact_k ! constant factor
! double precision, intent(out) :: P_new(0:max_dim,3)! polynomial
! integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials
!
! double precision :: P_center_tmp(3) ! new center
! double precision :: p_tmp ! new exponent
! double precision :: fact_k_tmp,fact_k_bis ! constant factor
! double precision :: P_new_tmp(0:max_dim,3)! polynomial
! integer :: i,j
! double precision :: binom_func
!
! ! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp
! call give_explicit_cpoly_and_cgaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim)
! ! Then you create the new gaussian from the product of the new one per the Nuclei one
! call cgaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center)
! fact_k = fact_k_bis * fact_k_tmp
!
! ! Then you build the coefficient of the new polynom
! do i = 0, iorder(1)
! P_new(i,1) = 0.d0
! do j = i,iorder(1)
! P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i)
! enddo
! enddo
! do i = 0, iorder(2)
! P_new(i,2) = 0.d0
! do j = i,iorder(2)
! P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i)
! enddo
! enddo
! do i = 0, iorder(3)
! P_new(i,3) = 0.d0
! do j = i,iorder(3)
! P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i)
! enddo
! enddo
!
!end
! ---
subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
BEGIN_DOC
@ -237,7 +175,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
ab = a * b * p_inv
k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3))
tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k))
tmp_mod = dsqrt(real(k)*real(k) + aimag(k)*aimag(k))
if(tmp_mod .gt. 40.d0) then
k = (0.d0, 0.d0)
xp(1:3) = (0.d0, 0.d0)
@ -309,8 +247,6 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd)
integer, intent(out) :: nd
integer :: ndtmp, ib, ic
double precision :: tmp_mod
complex*16 :: tmp
if(ior(nc, nb) >= 0) then ! True if nc>=0 and nb>=0
continue
@ -332,9 +268,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd)
enddo
do nd = ndtmp, 0, -1
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(tmp_mod .lt. 1.d-15) cycle
if(abs(d(nd)) .lt. 1.d-15) cycle
exit
enddo
@ -432,47 +366,42 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b)
complex*16, intent(in) :: x_A, x_P, x_B, x_Q
complex*16, intent(out) :: P_A(0:a), P_B(0:b)
integer :: i, minab, maxab
complex*16 :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
integer :: i
integer :: maxbinom
complex*16 :: pows_a(0:a+b+2), pows_b(0:a+b+2)
double precision :: binom_func
if((a < 0) .or. (b < 0)) return
maxab = max(a, b)
minab = max(min(a, b), 0)
maxbinom = size(binom_transp, 1)
pows_a(0) = (1.d0, 0.d0)
pows_a(1) = x_P - x_A
do i = 2, a
pows_a(i) = pows_a(i-1) * pows_a(1)
enddo
pows_b(0) = (1.d0, 0.d0)
pows_b(1) = x_Q - x_B
do i = 2, maxab
pows_a(i) = pows_a(i-1) * pows_a(1)
do i = 2, b
pows_b(i) = pows_b(i-1) * pows_b(1)
enddo
P_A(0) = pows_a(a)
do i = 1, min(a, maxbinom)
P_A(i) = binom_transp(i,a) * pows_a(a-i)
enddo
do i = maxbinom+1, a
P_A(i) = binom_func(a, i) * pows_a(a-i)
enddo
P_B(0) = pows_b(b)
do i = 1, min(minab, 20)
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
do i = 1, min(b, maxbinom)
P_B(i) = binom_transp(i,b) * pows_b(b-i)
enddo
do i = minab+1, min(a, 20)
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
enddo
do i = minab+1, min(b, 20)
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = 101, a
P_A(i) = binom_func(a,a-i) * pows_a(a-i)
enddo
do i = 101, b
P_B(i) = binom_func(b,b-i) * pows_b(b-i)
do i = maxbinom+1, b
P_B(i) = binom_func(b, i) * pows_b(b-i)
enddo
end

View File

@ -68,66 +68,6 @@ end
! ---
complex*16 function crint_quad(n, rho)
implicit none
integer, intent(in) :: n
complex*16, intent(in) :: rho
integer :: i_quad, n_quad
double precision :: tmp_inv, tmp
n_quad = 1000000000
tmp_inv = 1.d0 / dble(n_quad)
!crint_quad = 0.5d0 * zexp(-rho)
!do i_quad = 1, n_quad - 1
! tmp = tmp_inv * dble(i_quad)
! tmp = tmp * tmp
! crint_quad += zexp(-rho*tmp) * tmp**n
!enddo
!crint_quad = crint_quad * tmp_inv
!crint_quad = 0.5d0 * zexp(-rho)
!do i_quad = 1, n_quad - 1
! tmp = tmp_inv * dble(i_quad)
! crint_quad += zexp(-rho*tmp) * tmp**n / dsqrt(tmp)
!enddo
!crint_quad = crint_quad * 0.5d0 * tmp_inv
! Composite Boole's Rule
crint_quad = 7.d0 * zexp(-rho)
do i_quad = 1, n_quad - 1
tmp = tmp_inv * dble(i_quad)
tmp = tmp * tmp
if(modulo(i_quad, 4) .eq. 0) then
crint_quad += 14.d0 * zexp(-rho*tmp) * tmp**n
else if(modulo(i_quad, 2) .eq. 0) then
crint_quad += 12.d0 * zexp(-rho*tmp) * tmp**n
else
crint_quad += 32.d0 * zexp(-rho*tmp) * tmp**n
endif
enddo
crint_quad = crint_quad * 2.d0 * tmp_inv / 45.d0
! Composite Simpson's 3/8 rule
!crint_quad = zexp(-rho)
!do i_quad = 1, n_quad - 1
! tmp = tmp_inv * dble(i_quad)
! tmp = tmp * tmp
! if(modulo(i_quad, 3) .eq. 0) then
! crint_quad += 2.d0 * zexp(-rho*tmp) * tmp**n
! else
! crint_quad += 3.d0 * zexp(-rho*tmp) * tmp**n
! endif
!enddo
!crint_quad = crint_quad * 3.d0 * tmp_inv / 8.d0
end
! ---
complex*16 function crint_sum_1(n_pt_out, rho, d1)
implicit none
@ -286,17 +226,18 @@ complex*16 function crint_2(n, 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)
tmp = dble(n + n)
rho2 = rho * rho
crint_2 = 1.d0 / (tmp + 1.d0) &
- rho / (tmp + 3.d0) &
crint_2 = - 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0) &
+ 0.5d0 * rho2 / (tmp + 5.d0) &
- 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0)
- rho / (tmp + 3.d0) &
+ 1.d0 / (tmp + 1.d0)
else
crint_2 = crint_smallz(n, rho)
endif
@ -328,6 +269,9 @@ subroutine zboysfun(n_max, x, vals)
! Input: x --- argument, complex*16, Re(x) >= 0
! Output: vals --- values of the Boys function, n = 0, 1, ..., n_max
!
! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
! https://doi.org/10.1063/5.0062444
!
END_DOC
implicit none
@ -363,6 +307,9 @@ subroutine zboysfunnrp(n_max, x, vals)
! Input: x --- argument, complex *16 Re(x)<=0
! Output: vals --- values of e^z F(n,z), n = 0, 1, ..., n_max
!
! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
! https://doi.org/10.1063/5.0062444
!
END_DOC
implicit none
@ -398,17 +345,12 @@ complex*16 function crint_sum_2(n_pt_out, rho, d1)
complex*16, allocatable :: vals(:)
!complex*16, external :: crint_2
!crint_sum_2 = (0.d0, 0.d0)
!do i = 0, n_pt_out, 2
! n = shiftr(i, 1)
! crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho)
!enddo
n_max = shiftr(n_pt_out, 1)
allocate(vals(0:n_max))
call crint_2_vec(n_max, rho, vals)
! FOR DEBUG
!call crint_quad_12_vec(n_max, rho, vals)
crint_sum_2 = d1(0) * vals(0)
do i = 2, n_pt_out, 2
@ -444,16 +386,22 @@ subroutine crint_2_vec(n_max, rho, vals)
! use finite expansion for very small rho
! rho^2 / 2
rho2 = 0.5d0 * rho * rho
! rho^3 / 6
rho3 = 0.3333333333333333d0 * rho * rho2
!! rho^2 / 2
!rho2 = 0.5d0 * rho * rho
!! rho^3 / 6
!rho3 = 0.3333333333333333d0 * rho * rho2
!vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3
!do n = 1, n_max
! tmp = 2.d0 * dble(n)
! vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) &
! + rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0)
!enddo
vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3
! TODO (last term)
vals(0) = 1.d0 + rho * (-0.3333333333333333d0 + rho*(0.1d0 - 0.047619047619047616d0*rho))
do n = 1, n_max
tmp = 2.d0 * dble(n)
vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) &
+ rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0)
vals(n) = 1.d0 / (tmp + 1.d0) + rho * (-1.d0/(tmp+3.d0) + rho*(0.5d0/(tmp+5.d0) - 0.3333333333333333d0*rho/(tmp+7.d0)))
enddo
else
@ -541,3 +489,268 @@ end
! ---
subroutine crint_quad_1(n, rho, n_quad, crint_quad)
implicit none
integer, intent(in) :: n, n_quad
complex*16, intent(in) :: rho
complex*16, intent(out) :: crint_quad
integer :: i_quad
double precision :: tmp_inv, tmp0, tmp1, tmp2
double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /)
tmp_inv = 1.d0 / dble(n_quad)
crint_quad = 7.d0 * zexp(-rho)
tmp0 = 0.d0
select case (n)
case (0)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1)
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (1)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (2)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 * tmp1
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (3)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 * tmp1 * tmp1
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (4)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
tmp2 = tmp1 * tmp1
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp2 * tmp2
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case default
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1**n
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
end select
end
! ---
subroutine crint_quad_2(n, rho, n_quad, crint_quad)
implicit none
integer, intent(in) :: n, n_quad
complex*16, intent(in) :: rho
complex*16, intent(out) :: crint_quad
integer :: i_quad
double precision :: tmp_inv, tmp0, tmp1, tmp2
double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /)
complex*16 :: rhoc, rhoe
tmp_inv = 1.d0 / dble(n_quad)
crint_quad = 7.d0 * zexp(-rho)
tmp0 = 0.d0
rhoc = zexp(-rho*tmp_inv)
rhoe = (1.d0, 0.d0)
select case (n)
case (0)
!do i_quad = 1, n_quad - 1
! tmp0 = tmp0 + tmp_inv
! rhoe = rhoe * rhoc
! tmp1 = (rhoe - 1.d0) / dsqrt(tmp0)
! crint_quad = crint_quad + coef(iand(i_quad, 3)) * tmp1
!enddo
!crint_quad = 1.d0 + crint_quad * 0.022222222222222223d0 * tmp_inv
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe / dsqrt(tmp0)
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (1)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (2)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (3)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (4)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
tmp2 = tmp1 * tmp1 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp2
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case default
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0**n / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
end select
end
! ---
subroutine crint_quad_12(n, rho, n_quad, crint_quad)
implicit none
integer, intent(in) :: n, n_quad
complex*16, intent(in) :: rho
complex*16, intent(out) :: crint_quad
integer :: i_quad
double precision :: tmp_inv, tmp0, tmp1, tmp2
double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /)
complex*16 :: rhoc, rhoe
tmp_inv = 1.d0 / dble(n_quad)
crint_quad = 7.d0 * zexp(-rho)
tmp0 = 0.d0
rhoc = zexp(-rho*tmp_inv)
rhoe = (1.d0, 0.d0)
select case (n)
case (0)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1)
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (1)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (2)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (3)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (4)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
tmp2 = tmp1 * tmp1 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp2
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case default
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0**n / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
end select
end
! ---
subroutine crint_quad_12_vec(n_max, rho, vals)
implicit none
integer, intent(in) :: n_max
complex*16, intent(in) :: rho
complex*16, intent(out) :: vals(0:n_max)
integer :: n
double precision :: tmp, abs_rho
complex*16 :: rho2, rho3, erho
do n = 1, n_max
call crint_quad_12(n, rho, 10000000, vals(n))
enddo
return
end
! ---

View File

@ -340,8 +340,8 @@ subroutine zboysfun00(z, val)
if(abs(z) .ge. 100.0d0) then
! large |z|
z1 = 1.0d0 / zsqrt(z)
y = 1.0d0 / z
z1 = (1.0d0, 0.d0) / zsqrt(z)
y = (1.0d0, 0.d0) / z
val = asymcoef(7)
do k = 6, 1, -1
val = val * y + asymcoef(k)
@ -389,6 +389,8 @@ subroutine zboysfun00nrp(z, val)
!
END_DOC
include 'constants.include.F'
implicit none
double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, &
@ -413,7 +415,6 @@ subroutine zboysfun00nrp(z, val)
double precision, parameter :: tol = 1.0d-03
double precision, parameter :: sqpio2 = 0.886226925452758014d0 ! sqrt(pi)/2
double precision, parameter :: pi = 3.14159265358979324d0
double precision, parameter :: etmax = 25.7903399171930621d0
double precision, parameter :: etmax1 = 26.7903399171930621d0
complex*16, parameter :: ima = (0.d0, 1.d0)
@ -532,8 +533,8 @@ subroutine zboysfun00nrp(z, val)
if(abs(z) .ge. 100.0d0) then
! large |z|
z1 = 1.0d0 / zsqrt(z)
y = 1.0d0 / z
z1 = (1.0d0, 0.d0) / zsqrt(z)
y = (1.0d0, 0.d0) / z
val = asymcoef(7)
do k = 6, 1, -1
val = val * y + asymcoef(k)
@ -560,13 +561,13 @@ subroutine zboysfun00nrp(z, val)
zsum = zsum + ww(k) * (zz - uu(k)) / (qq(k) + z)
else
q = z + qq(k)
p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
p = q * (0.041666666666666664d0*q * (q * (0.2d0*q - 1.d0) + 4.d0) - 0.5d0) + 1.d0
zsum = zsum + ww(k) * p * zz
endif
enddo
zt = ima * sqrt(z / etmax)
zt = ima * zsqrt(z / etmax)
tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt))
val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z)
val = dsqrt(etmax) * zsum * inv_sq_pi + zz * tmp / zsqrt(pi*z)
else
zsum = (0.d0, 0.d0)
do k = 1, 16
@ -574,13 +575,14 @@ subroutine zboysfun00nrp(z, val)
zsum = zsum + ww(k) * (zz - uu1(k)) / (qq1(k) + z)
else
q = z + qq1(k)
p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
!p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
p = q * (0.041666666666666664d0*q * (q * (0.2d0*q - 1.d0) + 4.d0) - 0.5d0) + 1.d0
zsum = zsum + ww(k) * p * zz
endif
enddo
zt = ima * zsqrt(z / etmax1)
tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt))
val = dsqrt(etmax1) * zsum / dsqrt(pi) + zz * tmp / zsqrt(pi*z)
val = dsqrt(etmax1) * zsum * inv_sq_pi + zz * tmp / zsqrt(pi*z)
endif
return

View File

@ -1008,50 +1008,78 @@ subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, n_points)
deallocate(pows_a)
!deallocate(pows_b)
end subroutine recentered_poly2_v0
end
! ---
subroutine recentered_poly2(P_new, x_A, x_P, a, Q_new, x_B, x_Q, b)
subroutine recentered_poly2(P_new,x_A,x_P,a,P_new2,x_B,x_Q,b)
implicit none
BEGIN_DOC
! Recenter two polynomials
!
! Recenter two polynomials:
!
! (x - x_A)^a -> \sum_{i=0}^{a} \binom{a}{i} (x_A - x_P)^{a-i} (x - x_P)^i
! (x - x_B)^b -> \sum_{i=0}^{b} \binom{b}{i} (x_B - x_Q)^{b-i} (x - x_Q)^i
!
! where:
! \binom{a}{i} = \binom{a}{a-i} = a! / [i! (a-i)!]
!
END_DOC
implicit none
integer, intent(in) :: a, b
double precision, intent(in) :: x_A, x_P, x_B, x_Q
double precision, intent(out) :: P_new(0:a),P_new2(0:b)
double precision, intent(out) :: P_new(0:a), Q_new(0:b)
double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
double precision :: binom_func
integer :: i,j,k,l, minab, maxab
integer :: i, minab, maxab
integer :: maxbinom
double precision, external :: binom_func
if((a < 0) .or. (b < 0)) return
maxbinom = size(binom_transp, 1)
maxab = max(a, b)
minab = max(min(a,b),0)
minab = min(min(a, b), maxbinom)
pows_a(0) = 1.d0
pows_a(1) = (x_P - x_A)
pows_a(1) = x_P - x_A
pows_b(0) = 1.d0
pows_b(1) = (x_Q - x_B)
pows_b(1) = x_Q - x_B
do i = 2, maxab
pows_a(i) = pows_a(i-1) * pows_a(1)
pows_b(i) = pows_b(i-1) * pows_b(1)
enddo
P_new(0) = pows_a(a)
P_new2(0) = pows_b(b)
do i = 1,min(minab,20)
P_new (i) = binom_transp(a-i,a) * pows_a(a-i)
P_new2(i) = binom_transp(b-i,b) * pows_b(b-i)
Q_new(0) = pows_b(b)
do i = 1, minab
P_new(i) = binom_transp(i,a) * pows_a(a-i)
Q_new(i) = binom_transp(i,b) * pows_b(b-i)
enddo
do i = minab+1,min(a,20)
P_new (i) = binom_transp(a-i,a) * pows_a(a-i)
do i = minab+1, min(a, maxbinom)
P_new(i) = binom_transp(i,a) * pows_a(a-i)
enddo
do i = minab+1,min(b,20)
P_new2(i) = binom_transp(b-i,b) * pows_b(b-i)
do i = minab+1, min(b, maxbinom)
Q_new(i) = binom_transp(i,b) * pows_b(b-i)
enddo
do i = 101,a
P_new(i) = binom_func(a,a-i) * pows_a(a-i)
do i = maxbinom+1, a
P_new(i) = binom_func(a, i) * pows_a(a-i)
enddo
do i = 101,b
P_new2(i) = binom_func(b,b-i) * pows_b(b-i)
do i = maxbinom+1, b
Q_new(i) = binom_func(b, i) * pows_b(b-i)
enddo
end
! ---
subroutine pol_modif_center(A_center, B_center, iorder, A_pol, B_pol)
BEGIN_DOC