diff --git a/scripts/import_champ_jastrow.py b/scripts/import_champ_jastrow.py old mode 100755 new mode 100644 index 489309b7..4bc58ff4 --- a/scripts/import_champ_jastrow.py +++ b/scripts/import_champ_jastrow.py @@ -45,7 +45,7 @@ if __name__ == '__main__': jastrow_file = sys.argv[2] jastrow = import_jastrow(jastrow_file) 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']) charges = ezfio.get_nuclei_nucl_charge() types = {} diff --git a/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f b/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f index 7f94f226..88fe4c85 100644 --- a/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f +++ b/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f @@ -111,8 +111,9 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a complex*16 :: accu, P_center(3) complex*16 :: d(0:n_pt_in) - complex*16 :: V_n_e_cosgtos - complex*16 :: crint + complex*16, external :: V_n_e_cosgtos + complex*16, external :: crint_2 + complex*16, external :: crint_sum_2 if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. & (A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then @@ -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)) ) if(n_pt == 0) then - NAI_pol_mult_cosgtos = coeff * crint(0, const) + NAI_pol_mult_cosgtos = coeff * crint_2(0, const) return endif - call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta & - , power_A, power_B, C_center, n_pt_in, d, n_pt_out) + call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & + power_A, power_B, C_center, n_pt_in, d, n_pt_out) if(n_pt_out < 0) then NAI_pol_mult_cosgtos = (0.d0, 0.d0) return endif - accu = (0.d0, 0.d0) - do i = 0, n_pt_out, 2 - accu += crint(shiftr(i, 1), const) * d(i) - -! print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const)) - enddo + !accu = (0.d0, 0.d0) + !do i = 0, n_pt_out, 2 + ! accu += crint_2(shiftr(i, 1), const) * d(i) + !enddo + accu = crint_sum_2(n_pt_out, const, d) NAI_pol_mult_cosgtos = accu * coeff -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) 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 -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 -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) 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) endif -end function V_r_cosgtos +end ! --- diff --git a/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f b/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f index ea9ff009..df402ff1 100644 --- a/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f +++ b/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f @@ -35,11 +35,9 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then - !print *, ' with shwartz acc ' ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l) else - !print *, ' without shwartz acc ' dim1 = n_pt_max_integrals @@ -51,7 +49,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) ao_two_e_integral_cosgtos = 0.d0 if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then - !print *, ' not the same center' do p = 1, 3 I_power(p) = ao_power(i,p) @@ -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) expo2 = ao_expo_ord_transp_cosgtos(q,j) - call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 & - , expo1, expo2, I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & + expo1, expo2, I_power, J_power, I_center, J_center, dim1) p1_inv = (1.d0,0.d0) / pp1 - call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & - , conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & + conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1) p2_inv = (1.d0,0.d0) / pp2 - call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & - , expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, & + expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1) p3_inv = (1.d0,0.d0) / pp3 - call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & - , conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, & + conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1) p4_inv = (1.d0,0.d0) / pp4 - !integer :: ii - !do ii = 1, 3 - ! print *, 'fact_p1', fact_p1 - ! print *, 'fact_p2', fact_p2 - ! print *, 'fact_p3', fact_p3 - ! print *, 'fact_p4', fact_p4 - ! !print *, pp1, p1_inv - ! !print *, pp2, p2_inv - ! !print *, pp3, p3_inv - ! !print *, pp4, p4_inv - !enddo - ! if( abs(aimag(P1_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_1 is complex !!' - ! print *, P1_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - ! if( abs(aimag(P2_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_2 is complex !!' - ! print *, P2_center - ! print *, ' old expos:' - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! print *, ' new expo:' - ! print *, pp2, p2_inv - ! print *, ' factor:' - ! print *, fact_p2 - ! print *, ' old centers:' - ! print *, I_center, J_center - ! print *, ' powers:' - ! print *, I_power, J_power - ! stop - ! endif - ! if( abs(aimag(P3_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_3 is complex !!' - ! print *, P3_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - ! if( abs(aimag(P4_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_4 is complex !!' - ! print *, P4_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - !enddo - do r = 1, ao_prim_num(k) coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k) expo3 = ao_expo_ord_transp_cosgtos(r,k) @@ -146,74 +93,47 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) expo4 = ao_expo_ord_transp_cosgtos(s,l) - call give_explicit_cpoly_and_cgaussian( Q1_new, Q1_center, qq1, fact_q1, iorder_q1 & - , expo3, expo4, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & + expo3, expo4, K_power, L_power, K_center, L_center, dim1) q1_inv = (1.d0,0.d0) / qq1 - call give_explicit_cpoly_and_cgaussian( Q2_new, Q2_center, qq2, fact_q2, iorder_q2 & - , conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & + conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1) q2_inv = (1.d0,0.d0) / qq2 - !do ii = 1, 3 - ! !print *, qq1, q1_inv - ! !print *, qq2, q2_inv - ! print *, 'fact_q1', fact_q1 - ! print *, 'fact_q2', fact_q2 - !enddo - ! if( abs(aimag(Q1_center(ii))) .gt. 0.d0 ) then - ! print *, ' Q_1 is complex !!' - ! print *, Q1_center - ! print *, expo3, expo4 - ! print *, conjg(expo3), conjg(expo4) - ! stop - ! endif - ! if( abs(aimag(Q2_center(ii))) .gt. 0.d0 ) then - ! print *, ' Q_2 is complex !!' - ! print *, Q2_center - ! print *, expo3, expo4 - ! print *, conjg(expo3), conjg(expo4) - ! stop - ! endif - !enddo + integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) - - integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) - - integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 - !integral_tot = integral1 - !print*, integral_tot - ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot) enddo ! s - enddo ! r - enddo ! q - enddo ! p + enddo ! r + enddo ! q + enddo ! p else - !print *, ' the same center' do p = 1, 3 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 -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 :: schwartz_ij - complex*16 :: ERI_cosgtos - complex*16 :: general_primitive_integral_cosgtos + complex*16, external :: ERI_cosgtos + complex*16, external :: general_primitive_integral_cosgtos ao_2e_cosgtos_schwartz_accel = 0.d0 @@ -341,7 +261,7 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) thr = ao_integrals_threshold*ao_integrals_threshold - allocate( schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)) ) + allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then @@ -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) expo2 = ao_expo_ord_transp_cosgtos(s,l) - call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 & - , expo1, expo2, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & + expo1, expo2, K_power, L_power, K_center, L_center, dim1) p1_inv = (1.d0,0.d0) / pp1 - call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & - , conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & + conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1) p2_inv = (1.d0,0.d0) / pp2 - call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & - , expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, & + expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1) p3_inv = (1.d0,0.d0) / pp3 - call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & - , conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, & + conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1) p4_inv = (1.d0,0.d0) / pp4 - integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -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) expo2 = ao_expo_ord_transp_cosgtos(s,l) - integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -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) expo2 = ao_expo_ord_transp_cosgtos(q,j) - integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -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) 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 @@ -707,11 +631,11 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) deallocate(schwartz_kl) -end function ao_2e_cosgtos_schwartz_accel +end ! --- -BEGIN_PROVIDER [ double precision, ao_2e_cosgtos_schwartz, (ao_num,ao_num)] +BEGIN_PROVIDER [double precision, ao_2e_cosgtos_schwartz, (ao_num, ao_num)] BEGIN_DOC ! Needed to compute Schwartz inequalities @@ -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 & - , Q_new, Q_center, fact_q, q, q_inv, iorder_q ) +complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fact_p, p, p_inv, iorder_p, & + Q_new, Q_center, fact_q, q, q_inv, iorder_q) BEGIN_DOC ! @@ -765,7 +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 :: d1(0:max_dim), d_poly(0:max_dim) - complex*16 :: crint_sum + complex*16 :: crint_sum_2 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol @@ -912,13 +836,11 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa !DIR$ FORCEINLINE call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) - accu = crint_sum(n_pt_out, const, d1) -! print *, n_pt_out, real(d1(0:n_pt_out)) -! print *, real(accu) + accu = crint_sum_2(n_pt_out, const, d1) 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 -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) 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 -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 -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 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 -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 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 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 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 subroutine I_x2_pol_mult_cosgtos +end ! --- diff --git a/src/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f new file mode 100644 index 00000000..469eb654 --- /dev/null +++ b/src/hartree_fock/deb_ao_2e_int.irp.f @@ -0,0 +1,188 @@ + +program deb_ao_2e_int + + !call check_ao_two_e_integral_cosgtos() + call check_crint1() + !call check_crint2() + +end + +! --- + +subroutine check_ao_two_e_integral_cosgtos() + + implicit none + + integer :: i, j, k, l + double precision :: tmp1, tmp2 + double precision :: acc, nrm, dif + + double precision, external :: ao_two_e_integral + double precision, external :: ao_two_e_integral_cosgtos + + acc = 0.d0 + nrm = 0.d0 + + i = 1 + j = 6 + k = 1 + l = 16 +! do i = 1, ao_num +! do k = 1, ao_num +! do j = 1, ao_num +! do l = 1, ao_num + + tmp1 = ao_two_e_integral (i, j, k, l) + tmp2 = ao_two_e_integral_cosgtos(i, j, k, l) + + dif = dabs(tmp1 - tmp2) + if(dif .gt. 1d-12) then + print*, ' error on:', i, j, k, l + print*, tmp1, tmp2, dif + stop + endif + +! acc += dif +! nrm += dabs(tmp1) +! enddo +! enddo +! enddo +! enddo + + print *, ' acc (%) = ', dif * 100.d0 / nrm + +end + +! --- + +subroutine check_crint1() + + implicit none + integer :: i, n, i_rho + double precision :: dif_thr + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), & + (+1d-9, +1d-6), & + (-1d-6, -1d-5), & + (+1d-3, -1d-2), & + (-1d-1, +1d-1), & + (+1d-0, +1d-1), & + (-1d+1, +1d+1), & + (+1d+2, +1d+1), & + (-1d+3, +1d+2), & + (+1d+4, +1d+4) /) + complex*16 :: rho + complex*16 :: int_an, int_nm + double precision, external :: rint + complex*16, external :: crint_1, crint_2, crint_quad + + n = 10 + dif_thr = 1d-7 + + do i_rho = 8, 10 + !do i_rho = 7, 7 + + !rho = (-10.d0, 0.1d0) + !rho = (+10.d0, 0.1d0) + rho = rho_test(i_rho) + print*, "rho = ", real(rho), aimag(rho) + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i = 0, n + !int_an = crint_1 (i, rho) + int_an = crint_2 (i, rho) + int_nm = crint_quad(i, rho) + + dif_re = dabs(real(int_an) - real(int_nm)) + dif_im = dabs(aimag(int_an) - aimag(int_nm)) + + if((dif_re .gt. dif_thr) .or. (dif_im .gt. dif_thr)) then + print*, ' error on i =', i + print*, real(int_an), real(int_nm), dif_re + print*, aimag(int_an), aimag(int_nm), dif_im + !print*, rint(i, real(rho)) + print*, crint_1(i, rho) + !print*, crint_2(i, rho) + stop + endif + acc_re += dif_re + nrm_re += dabs(real(int_nm)) + acc_im += dif_im + nrm_im += dabs(aimag(int_nm)) + enddo + + print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15) + print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15) + enddo + +end + +! --- + +subroutine check_crint2() + + implicit none + + integer :: i, n, i_rho + double precision :: dif_thr + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), & + (+1d-9, +1d-6), & + (-1d-6, -1d-5), & + (+1d-3, -1d-2), & + (-1d-1, +1d-1), & + (+1d-0, +1d-1), & + (-1d+1, +1d+1), & + (+1d+2, +1d+1), & + (-1d+3, +1d+2), & + (+1d+4, +1d+4) /) + complex*16 :: rho + complex*16 :: int_an, int_nm + complex*16, external :: crint_1, crint_2 + + n = 30 + dif_thr = 1d-12 + + do i_rho = 1, 10 + rho = rho_test(i_rho) + print*, "rho = ", real(rho), aimag(rho) + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i = 0, n + int_an = crint_1(i, rho) + int_nm = crint_2(i, rho) + + dif_re = dabs(real(int_an) - real(int_nm)) + !if(dif_re .gt. dif_thr) then + ! print*, ' error in real part:', i + ! print*, real(int_an), real(int_nm), dif_re + ! stop + !endif + acc_re += dif_re + nrm_re += dabs(real(int_nm)) + + dif_im = dabs(aimag(int_an) - aimag(int_nm)) + !if(dif_im .gt. dif_thr) then + ! print*, ' error in imag part:', i + ! print*, aimag(int_an), aimag(int_nm), dif_im + ! stop + !endif + acc_im += dif_im + nrm_im += dabs(aimag(int_nm)) + enddo + + print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15) + print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15) + enddo + +end + +! --- + + diff --git a/src/hartree_fock/print_scf_int.irp.f b/src/hartree_fock/print_scf_int.irp.f deleted file mode 100644 index ee7590f6..00000000 --- a/src/hartree_fock/print_scf_int.irp.f +++ /dev/null @@ -1,114 +0,0 @@ - -program print_scf_int - - call main() - -end - -subroutine main() - - implicit none - integer :: i, j, k, l - - print *, " Hcore:" - do j = 1, ao_num - do i = 1, ao_num - print *, i, j, ao_one_e_integrals(i,j) - enddo - enddo - - print *, " P:" - do j = 1, ao_num - do i = 1, ao_num - print *, i, j, SCF_density_matrix_ao_alpha(i,j) - enddo - enddo - - - double precision :: integ, density_a, density_b, density - double precision :: J_scf(ao_num, ao_num) - double precision :: K_scf(ao_num, ao_num) - - - double precision, external :: get_ao_two_e_integral - PROVIDE ao_integrals_map - - print *, " J:" - !do j = 1, ao_num - ! do l = 1, ao_num - ! do i = 1, ao_num - ! do k = 1, ao_num - ! ! < 1:k, 2:l | 1:i, 2:j > - ! print *, '< k l | i j >', k, l, i, j - ! print *, get_ao_two_e_integral(i, j, k, l, ao_integrals_map) - ! enddo - ! enddo - ! enddo - !enddo - - !do k = 1, ao_num - ! do i = 1, ao_num - ! do j = 1, ao_num - ! do l = 1, ao_num - ! ! ( 1:k, 1:i | 2:l, 2:j ) - ! print *, '(k i | l j)', k, i, l, j - ! print *, get_ao_two_e_integral(l, j, k, i, ao_integrals_map) - ! enddo - ! enddo - ! print *, '' - ! enddo - !enddo - - J_scf = 0.d0 - K_scf = 0.d0 - do i = 1, ao_num - do k = 1, ao_num - do j = 1, ao_num - do l = 1, ao_num - - density_a = SCF_density_matrix_ao_alpha(l,j) - density_b = SCF_density_matrix_ao_beta (l,j) - density = density_a + density_b - - integ = get_ao_two_e_integral(l, j, k, i, ao_integrals_map) - J_scf(k,i) += density * integ - integ = get_ao_two_e_integral(l, i, k, j, ao_integrals_map) - K_scf(k,i) -= density_a * integ - enddo - enddo - enddo - enddo - - print *, 'J x P' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, J_scf(k,i) - enddo - enddo - - print *, '' - print *, 'K x P' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, K_scf(k,i) - enddo - enddo - - print *, '' - print *, 'F in AO' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, Fock_matrix_ao(k,i) - enddo - enddo - - print *, '' - print *, 'F in MO' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, 2.d0 * Fock_matrix_mo_alpha(k,i) - enddo - enddo - -end - diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index a820d5f2..9c25edc2 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -56,7 +56,7 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new) iorder = a + b -end subroutine give_explicit_cpoly_and_cgaussian_x +end ! --- @@ -141,7 +141,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, !DIR$ FORCEINLINE call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new) -end subroutine give_explicit_cpoly_and_cgaussian +end ! --- @@ -249,7 +249,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv -end subroutine cgaussian_product +end ! --- @@ -290,7 +290,7 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) k = zexp(-k) xp = (a*xa + b*xb) * p_inv -end subroutine cgaussian_product_x +end ! --- @@ -338,7 +338,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) exit enddo -end subroutine multiply_cpoly +end ! --- @@ -373,7 +373,7 @@ subroutine add_cpoly(b, nb, c, nc, d, nd) if(nd < 0) exit enddo -end subroutine add_cpoly +end ! --- @@ -413,7 +413,7 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) endif -end subroutine add_cpoly_multiply +end ! --- @@ -475,7 +475,7 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b) P_B(i) = binom_func(b,b-i) * pows_b(b-i) enddo -end subroutine recentered_cpoly2 +end ! --- @@ -514,267 +514,7 @@ complex*16 function Fc_integral(n, inv_sq_p) Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1)) -end function Fc_integral - -! --- - -complex*16 function crint(n, rho) - - implicit none - include 'constants.include.F' - - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer :: i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - double precision :: n_tmp - complex*16 :: sq_rho, rho_inv, rho_exp - - complex*16 :: crint_smallz, cpx_erf - - rho_re = REAL (rho) - rho_im = AIMAG(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - - if(rho_mod < 10.d0) then - ! small z - - if(rho_mod .lt. 1.d-10) then - crint = 1.d0 / dble(n + n + 1) - else - crint = crint_smallz(n, rho) - endif - - else - ! large z - - if(rho_mod .gt. 40.d0) then - - n_tmp = dble(n) + 0.5d0 - crint = 0.5d0 * gamma(n_tmp) / (rho**n_tmp) - - else - - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - crint = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho - mmax = n - if(mmax .gt. 0) then - do i = 0, mmax-1 - crint = ((dble(i) + 0.5d0) * crint - rho_exp) * rho_inv - enddo - endif - - ! *** - - endif - - endif - -! print *, n, real(rho), real(crint) - -end function crint - -! --- - -complex*16 function crint_sum(n_pt_out, rho, d1) - - implicit none - include 'constants.include.F' - - integer, intent(in) :: n_pt_out - complex*16, intent(in) :: rho, d1(0:n_pt_out) - - integer :: n, i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - complex*16 :: sq_rho, F0 - complex*16 :: rho_tmp, rho_inv, rho_exp - complex*16, allocatable :: Fm(:) - - complex*16 :: crint_smallz, cpx_erf - - rho_re = REAL (rho) - rho_im = AIMAG(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - - if(rho_mod < 10.d0) then - ! small z - - if(rho_mod .lt. 1.d-10) then - -! print *, ' 111' -! print *, ' rho = ', rho - - crint_sum = d1(0) -! print *, 0, 1 - - do i = 2, n_pt_out, 2 - - n = shiftr(i, 1) - crint_sum = crint_sum + d1(i) / dble(n+n+1) - -! print *, n, 1.d0 / dble(n+n+1) - enddo - - ! *** - - else - -! print *, ' 222' -! print *, ' rho = ', real(rho) -! if(abs(aimag(rho)) .gt. 1d-15) then -! print *, ' complex rho', rho -! stop -! endif - - crint_sum = d1(0) * crint_smallz(0, rho) - -! print *, 0, real(d1(0)), real(crint_smallz(0, rho)) -! if(abs(aimag(d1(0))) .gt. 1d-15) then -! print *, ' complex d1(0)', d1(0) -! stop -! endif - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum = crint_sum + d1(i) * crint_smallz(n, rho) - -! print *, n, real(d1(i)), real(crint_smallz(n, rho)) -! if(abs(aimag(d1(i))) .gt. 1d-15) then -! print *, ' complex d1(i)', i, d1(i) -! stop -! endif - - enddo - -! print *, 'sum = ', real(crint_sum) -! if(abs(aimag(crint_sum)) .gt. 1d-15) then -! print *, ' complex crint_sum', crint_sum -! stop -! endif - - ! *** - - endif - - else - ! large z - - if(rho_mod .gt. 40.d0) then - -! print *, ' 333' -! print *, ' rho = ', rho - - rho_inv = (1.d0, 0.d0) / rho - rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv) - crint_sum = rho_tmp * d1(0) -! print *, 0, rho_tmp - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv - crint_sum = crint_sum + rho_tmp * d1(i) -! print *, n, rho_tmp - enddo - - ! *** - - else - -! print *, ' 444' -! print *, ' rho = ', rho - - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - !sq_rho = zsqrt(rho) - - - F0 = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho - crint_sum = F0 * d1(0) -! print *, 0, F0 - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - mmax = shiftr(n_pt_out, 1) - if(mmax .gt. 0) then - - allocate( Fm(mmax) ) - Fm(1:mmax) = (0.d0, 0.d0) - - do n = 0, mmax-1 - F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv - Fm(n+1) = F0 -! print *, n, F0 - enddo - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum = crint_sum + Fm(n) * d1(i) - enddo - deallocate(Fm) - endif - - ! *** - - endif - - endif - -end function crint_sum - -! --- - -complex*16 function crint_smallz(n, rho) - - BEGIN_DOC - ! Standard version of rint - END_DOC - - implicit none - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer, parameter :: kmax = 40 - double precision, parameter :: eps = 1.d-13 - - integer :: k - double precision :: delta_mod - complex*16 :: rho_k, ct, delta_k - - ct = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0) - rho_k = (1.d0, 0.d0) - crint_smallz = ct * rho_k / gamma(dble(n) + 1.5d0) - - do k = 1, kmax - - rho_k = rho_k * rho - delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0) - crint_smallz = crint_smallz + delta_k - - delta_mod = dsqrt(REAL(delta_k)*REAL(delta_k) + AIMAG(delta_k)*AIMAG(delta_k)) - if(delta_mod .lt. eps) return - enddo - - if(delta_mod > eps) then - write(*,*) ' pb in crint_smallz !' - write(*,*) ' n, rho = ', n, rho - write(*,*) ' delta_mod = ', delta_mod - stop 1 - endif - -end function crint_smallz +end ! --- diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f new file mode 100644 index 00000000..9ffcc817 --- /dev/null +++ b/src/utils/cpx_boys.irp.f @@ -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 + +! --- + diff --git a/src/utils/cpx_erf.irp.f b/src/utils/cpx_erf.irp.f index 61f81055..1c5fa61d 100644 --- a/src/utils/cpx_erf.irp.f +++ b/src/utils/cpx_erf.irp.f @@ -1,7 +1,7 @@ ! --- -complex*16 function cpx_erf(x, y) +complex*16 function cpx_erf_1(x, y) BEGIN_DOC ! @@ -25,25 +25,25 @@ complex*16 function cpx_erf(x, y) if(yabs .lt. 1.d-15) then - cpx_erf = (1.d0, 0.d0) * derf(x) + cpx_erf_1 = (1.d0, 0.d0) * derf(x) return else erf_tmp1 = (1.d0, 0.d0) * derf(x) erf_tmp2 = erf_E(x, yabs) + erf_F(x, yabs) - erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * ( erf_G(x, yabs) + erf_H(x, yabs) ) + erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * (erf_G(x, yabs) + erf_H(x, yabs)) erf_tot = erf_tmp1 + erf_tmp2 - erf_tmp3 endif if(y .gt. 0.d0) then - cpx_erf = erf_tot + cpx_erf_1 = erf_tot else - cpx_erf = CONJG(erf_tot) + cpx_erf_1 = conjg(erf_tot) endif -end function cpx_erf +end ! --- @@ -54,7 +54,7 @@ complex*16 function erf_E(x, yabs) double precision, intent(in) :: x, yabs - if( (dabs(x).gt.6.d0) .or. (x==0.d0) ) then + if((dabs(x).gt.6.d0) .or. (x==0.d0)) then erf_E = (0.d0, 0.d0) return endif @@ -70,7 +70,7 @@ complex*16 function erf_E(x, yabs) endif -end function erf_E +end ! --- @@ -109,7 +109,7 @@ double precision function erf_F(x, yabs) endif -end function erf_F +end ! --- @@ -149,7 +149,7 @@ complex*16 function erf_G(x, yabs) enddo -end function erf_G +end ! --- @@ -172,7 +172,7 @@ complex*16 function erf_H(x, yabs) endif - if( (dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0) ) then + if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then x2 = x * x ct = 0.5d0 * inv_pi @@ -186,7 +186,7 @@ complex*16 function erf_H(x, yabs) tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1 erf_H = erf_H + tmp2 - tmp_mod = dsqrt(REAL(tmp2)*REAL(tmp2) + AIMAG(tmp2)*AIMAG(tmp2)) + tmp_mod = dsqrt(real(tmp2)*real(tmp2) + aimag(tmp2)*aimag(tmp2)) if(tmp_mod .lt. 1d-15) exit enddo erf_H = ct * erf_H @@ -197,8 +197,394 @@ complex*16 function erf_H(x, yabs) endif -end function erf_H +end ! --- +complex*16 function cpx_erf_2(x, y) + + BEGIN_DOC + ! + ! compute erf(z) for z = x + i y + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, intent(in) :: x, y + + double precision :: yabs + complex*16 :: z + + yabs = dabs(y) + + if(yabs .lt. 1.d-15) then + + cpx_erf_2 = (1.d0, 0.d0) * derf(x) + return + + else + + z = x + (0.d0, 1.d0) * y + + if(x .ge. 0.d0) then + call zboysfun00(z, cpx_erf_2) + else + call zboysfun00nrp(z, cpx_erf_2) + cpx_erf_2 = cpx_erf_2 * zexp(-z) + endif + + endif + + return +end + +! --- + +subroutine zboysfun00(z, val) + + BEGIN_DOC + ! + ! Computes values of the Boys function for n=0 + ! for a complex valued argument + ! + ! Input: z --- argument, complex*16, Real(z) >= 0 + ! Output: val --- value of the Boys function n=0 + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & + 0.249999999999993161d0, & + -0.374999999999766599d0, & + 0.937499999992027020d0, & + -3.28124999972738868d0, & + 14.7656249906697030d0, & + -81.2109371803307752d0 /) + + double precision, parameter :: taylcoef(0:10) = (/ 1.0d0, & + -0.333333333333333333d0, & + 0.1d0, & + -0.238095238095238095d-01, & + 0.462962962962962963d-02, & + -0.757575757575757576d-03, & + 0.106837606837606838d-03, & + -0.132275132275132275d-04, & + 1.458916900093370682d-06, & + -1.450385222315046877d-07, & + 1.3122532963802805073d-08 /) + + double precision, parameter :: sqpio2 = 0.886226925452758014d0 + + double precision, parameter :: pp(1:22) = (/ 0.001477878263796956477d0, & + 0.013317276413725817441d0, & + 0.037063591452052541530d0, & + 0.072752512422882761543d0, & + 0.120236941228785688896d0, & + 0.179574293958937717967d0, & + 0.253534046984087292596d0, & + 0.350388652780721927513d0, & + 0.482109575931276669313d0, & + 0.663028993158374107103d0, & + 0.911814736856590885929d0, & + 1.2539502287919293d0, & + 1.7244634233573395d0, & + 2.3715248262781863d0, & + 3.2613796996078355d0, & + 4.485130169059591d0, & + 6.168062135122484d0, & + 8.48247187231787d0, & + 11.665305486296793d0, & + 16.042417132288328d0, & + 22.06192951814709d0, & + 30.340112094708307d0 /) + + double precision, parameter :: ff(1:22) = (/ 0.0866431027201416556d0, & + 0.0857720608434394764d0, & + 0.0839350436829178814d0, & + 0.0809661970413229146d0, & + 0.0769089548492978618d0, & + 0.0731552078711821626d0, & + 0.0726950035163157228d0, & + 0.0752842556089304050d0, & + 0.0770943953645196145d0, & + 0.0754250625677530441d0, & + 0.0689686192650315305d0, & + 0.05744480422143023d0, & + 0.04208199434694545d0, & + 0.025838539448223282d0, & + 0.012445024157255563d0, & + 0.004292541592599837d0, & + 0.0009354342987735969d0, & + 0.10840885466502504d-03, & + 5.271867966761674d-06, & + 7.765974039750418d-08, & + 2.2138172422680093d-10, & + 6.594161760037707d-14 /) + + complex*16, intent(in) :: z + complex*16, intent(out) :: val + + integer :: k + complex*16 :: z1, zz, y + + zz = zexp(-z) + + if(abs(z) .ge. 100.0d0) then + + ! large |z| + z1 = 1.0d0 / zsqrt(z) + y = 1.0d0 / z + val = asymcoef(7) + do k = 6, 1, -1 + val = val * y + asymcoef(k) + enddo + val = zz * val * y + z1 * sqpio2 + + else if(abs(z) .le. 0.35d0) then + + ! small |z| + val = taylcoef(10) * (1.d0, 0.d0) + do k = 9, 0, -1 + val = val * z + taylcoef(k) + enddo + + else + + ! intermediate |z| + val = sqpio2 / zsqrt(z) - 0.5d0 * zz * sum(ff(1:22)/(z+pp(1:22))) + !val = (0.d0, 0.d0) + !do k = 1, 22 + ! val += ff(k) / (z + pp(k)) + !enddo + !val = sqpio2 / zsqrt(z) - 0.5d0 * zz * val + + endif + + return +end + +! --- + +subroutine zboysfun00nrp(z, val) + + BEGIN_DOC + ! + ! Computes values of the exp(z) F(0,z) + ! (where F(0,z) is the Boys function) + ! for a complex valued argument with Real(z)<=0 + ! + ! Input: z --- argument, complex*16, !!! Real(z)<=0 !!! + ! Output: val --- value of the function !!! exp(z) F(0,z) !!!, where F(0,z) is the Boys function + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & + 0.249999999999993161d0, & + -0.374999999999766599d0, & + 0.937499999992027020d0, & + -3.28124999972738868d0, & + 14.7656249906697030d0, & + -81.2109371803307752d0 /) + + double precision, parameter :: taylcoef(0:10) = (/ 1.0d0, & + -0.333333333333333333d0, & + 0.1d0, & + -0.238095238095238095d-01, & + 0.462962962962962963d-02, & + -0.757575757575757576d-03, & + 0.106837606837606838d-03, & + -0.132275132275132275d-04, & + 1.458916900093370682d-06, & + -1.450385222315046877d-07, & + 1.3122532963802805073d-08 /) + + double precision, parameter :: tol = 1.0d-03 + double precision, parameter :: sqpio2 = 0.886226925452758014d0 ! sqrt(pi)/2 + double precision, parameter :: pi = 3.14159265358979324d0 + double precision, parameter :: etmax = 25.7903399171930621d0 + double precision, parameter :: etmax1 = 26.7903399171930621d0 + complex*16, parameter :: ima = (0.d0, 1.d0) + + double precision, parameter :: pp(1:16) = (/ 0.005299532504175031d0, & + 0.0277124884633837d0, & + 0.06718439880608407d0, & + 0.12229779582249845d0, & + 0.19106187779867811d0, & + 0.27099161117138637d0, & + 0.35919822461037054d0, & + 0.45249374508118123d0, & + 0.5475062549188188d0, & + 0.6408017753896295d0, & + 0.7290083888286136d0, & + 0.8089381222013219d0, & + 0.8777022041775016d0, & + 0.9328156011939159d0, & + 0.9722875115366163d0, & + 0.994700467495825d0 /) + + double precision, parameter :: ww(1:16) = (/ 0.013576229705876844d0, & + 0.03112676196932382d0, & + 0.04757925584124612d0, & + 0.062314485627766904d0, & + 0.07479799440828848d0, & + 0.08457825969750153d0, & + 0.09130170752246194d0, & + 0.0947253052275344d0, & + 0.0947253052275344d0, & + 0.09130170752246194d0, & + 0.08457825969750153d0, & + 0.07479799440828848d0, & + 0.062314485627766904d0, & + 0.04757925584124612d0, & + 0.03112676196932382d0, & + 0.013576229705876844d0 /) + + double precision, parameter :: qq (1:16) = (/ 0.0007243228510223928d0, & + 0.01980651726441906d0, & + 0.11641097769229371d0, & + 0.38573968881461146d0, & + 0.9414671037609641d0, & + 1.8939510935716377d0, & + 3.3275564293459383d0, & + 5.280587297262129d0, & + 7.730992222360452d0, & + 10.590207725831563d0, & + 13.706359477128965d0, & + 16.876705473663804d0, & + 19.867876155236257d0, & + 22.441333930203022d0, & + 24.380717439613566d0, & + 25.51771075067431d0 /) + + + double precision, parameter :: qq1 (1:16) = (/ 0.0007524078957852004d0,& + 0.020574499281252233d0, & + 0.12092472113522865d0, & + 0.40069643967765295d0, & + 0.9779717449089211d0, & + 1.9673875468969015d0, & + 3.4565797939091802d0, & + 5.485337886599723d0, & + 8.030755321535683d0, & + 11.000834641174064d0, & + 14.237812708111456d0, & + 17.531086359214406d0, & + 20.6382373144543d0, & + 23.31147887603379d0, & + 25.326060444703632d0, & + 26.507139770710722d0 /) + + double precision, parameter :: uu(1:16) = (/ 0.9992759394074501d0, & + 0.9803883431758104d0, & + 0.8901093330366746d0, & + 0.6799475005849274d0, & + 0.3900551639790145d0, & + 0.15047608763371934d0, & + 0.0358806749968974d0, & + 0.005089440900100864d0, & + 0.00043900830706867264d0, & + 0.000025161192619824898d0, & + 1.1153308427285078d-6, & + 4.68317018372038d-8, & + 2.3522908467181876d-9, & + 1.7941242138648815d-10, & + 2.5798173021885247d-11, & + 8.27559122014575d-12 /) + + + double precision, parameter :: uu1(1:16) = (/ 0.999247875092057d0, & + 0.979635711599488d0, & + 0.8861006617341018d0, & + 0.6698533710831932d0, & + 0.3760730980014839d0, & + 0.13982165701683388d0, & + 0.031537442321301304d0, & + 0.004147133581658446d0, & + 0.0003253024081883165d0, & + 0.000016687766678889653d0, & + 6.555359391864376d-7, & + 2.4341421258295026d-8, & + 1.0887481200652014d-9, & + 7.51542178140961d-11, & + 1.002378402152542d-11, & + 3.0767730761654096d-12 /) + + complex*16, intent(in) :: z + complex*16, intent(out) :: val + + integer :: k + complex*16 :: z1, zz, y, zsum, tmp, zt, q, p + + zz = zexp(z) + + if(abs(z) .ge. 100.0d0) then + ! large |z| + z1 = 1.0d0 / zsqrt(z) + y = 1.0d0 / z + val = asymcoef(7) + do k = 6, 1, -1 + val = val * y + asymcoef(k) + enddo + val = val * y + z1 * sqpio2 * zz + return + endif + + if(abs(z) .le. 0.35d0) then + ! small |z| + val = taylcoef(10) * (1.d0, 0.d0) + do k = 9, 0, -1 + val = val * z + taylcoef(k) + enddo + val = val * zz + return + endif + + if(abs(etmax+z) .ge. 0.5d0) then + ! intermediate |z| + zsum = (0.d0, 0.d0) + do k = 1, 16 + if(abs(z + qq(k)) .ge. tol) then + zsum = zsum + ww(k) * (zz - uu(k)) / (qq(k) + z) + else + q = z + qq(k) + p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 + zsum = zsum + ww(k) * p *zz + endif + enddo + zt = ima * sqrt(z / etmax) + tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) + val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z) + else + zsum = (0.d0, 0.d0) + do k = 1, 16 + if(abs(z + qq1(k)) .ge. tol) then + zsum = zsum + ww(k) * (zz - uu1(k)) / (qq1(k) + z) + else + q = z + qq1(k) + p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 + zsum = zsum + ww(k) * p * zz + endif + enddo + zt = ima * zsqrt(z / etmax1) + tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) + val = dsqrt(etmax1) * zsum / dsqrt(pi) + zz * tmp / zsqrt(pi*z) + endif + + return +end + +! ---