diff --git a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f index 9294b139..7240f4a5 100644 --- a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f @@ -12,16 +12,19 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] END_DOC implicit none - integer :: power_A(3), power_B(3) - integer :: i, j, k, l, m, n, ii, jj - double precision :: c, Z, C_center(3) - double precision :: phiA, KA2 - double precision :: phiB, KB2 - complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) - complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) - complex*16 :: C1, C2, I1, I2 - complex*16 :: NAI_pol_mult_cgtos + integer :: power_A(3), power_B(3) + integer :: i, j, k, l, m, n, ii, jj + double precision :: c, Z, C_center(3) + double precision :: phiA, KA2 + double precision :: phiB, KB2 + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) + complex*16 :: C1, C2, I1, I2 + + complex*16, external :: NAI_pol_mult_cgtos + + ao_integrals_n_e_cgtos = 0.d0 @@ -140,7 +143,6 @@ complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, a dist_AC += abs(Ae_center(i) - C_center(i) * (1.d0, 0.d0)) enddo - if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13) .or. use_pw) then continue @@ -158,24 +160,27 @@ complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, a p_inv = (1.d0, 0.d0) / p rho = alpha * beta * p_inv - dist = (0.d0, 0.d0) - dist_integral = (0.d0, 0.d0) - do i = 1, 3 - P_center(i) = (alpha * Ae_center(i) + beta * Be_center(i)) * p_inv - dist += (Ae_center(i) - Be_center(i)) * (Ae_center(i) - Be_center(i)) - dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) - enddo + dist = (Ae_center(1) - Be_center(1)) * (Ae_center(1) - Be_center(1)) & + + (Ae_center(2) - Be_center(2)) * (Ae_center(2) - Be_center(2)) & + + (Ae_center(3) - Be_center(3)) * (Ae_center(3) - Be_center(3)) const_factor = dist * rho - const = p * dist_integral - - if(abs(const_factor) > 80.d0) then + if(real(const_factor) > 80.d0) then NAI_pol_mult_cgtos = (0.d0, 0.d0) return endif + P_center(1) = (alpha * Ae_center(1) + beta * Be_center(1)) * p_inv + P_center(2) = (alpha * Ae_center(2) + beta * Be_center(2)) * p_inv + P_center(3) = (alpha * Ae_center(3) + beta * Be_center(3)) * p_inv + + dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) & + + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) & + + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3)) + + const = p * dist_integral factor = zexp(-const_factor) - coeff = dtwo_pi * factor * p_inv + coeff = dtwo_pi * factor * p_inv 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 @@ -214,12 +219,12 @@ subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & double precision, intent(in) :: C_center(3) complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3) integer, intent(out) :: n_pt_out - complex*16, intent(out) :: d(0:n_pt_in) + complex*16, intent(inout) :: d(0:n_pt_in) integer :: a_x, b_x, a_y, b_y, a_z, b_z integer :: n_pt1, n_pt2, n_pt3, dim, i, n_pt_tmp complex*16 :: p, P_center(3), rho, p_inv, p_inv_2 - complex*16 :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2) + complex*16 :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2) complex*16 :: d1(0:n_pt_in), d2(0:n_pt_in), d3(0:n_pt_in) ASSERT (n_pt_in > 1) @@ -228,9 +233,9 @@ subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & p_inv = (1.d0, 0.d0) / p p_inv_2 = 0.5d0 * p_inv - do i = 1, 3 - P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv - enddo + P_center(1) = (alpha * A_center(1) + beta * B_center(1)) * p_inv + P_center(2) = (alpha * A_center(2) + beta * B_center(2)) * p_inv + P_center(3) = (alpha * A_center(3) + beta * B_center(3)) * p_inv do i = 0, n_pt_in d(i) = (0.d0, 0.d0) @@ -257,6 +262,7 @@ subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & a_x = power_A(1) b_x = power_B(1) + call I_x1_pol_mult_one_e_cgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in) if(n_pt1 < 0) then @@ -281,6 +287,7 @@ subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & a_y = power_A(2) b_y = power_B(2) + call I_x1_pol_mult_one_e_cgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in) if(n_pt2 < 0) then @@ -305,6 +312,7 @@ subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & a_z = power_A(3) b_z = power_B(3) + call I_x1_pol_mult_one_e_cgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in) if(n_pt3 < 0) then @@ -319,11 +327,9 @@ subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & n_pt_tmp = 0 call multiply_cpoly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp) - do i = 0, n_pt_tmp - d1(i) = (0.d0, 0.d0) - enddo n_pt_out = 0 + d1(0:n_pt_tmp) = (0.d0, 0.d0) call multiply_cpoly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out) do i = 0, n_pt_out d(i) = d1(i) @@ -354,13 +360,13 @@ recursive subroutine I_x1_pol_mult_one_e_cgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt dim = n_pt_in - if( (a==0) .and. (c==0)) then + if((a == 0) .and. (c == 0)) then nd = 0 d(0) = (1.d0, 0.d0) return - elseif( (c < 0) .or. (nd < 0) ) then + elseif((c < 0) .or. (nd < 0)) then nd = -1 return diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f index 3f436743..339156db 100644 --- a/src/ao_one_e_ints/screening.irp.f +++ b/src/ao_one_e_ints/screening.irp.f @@ -3,7 +3,7 @@ logical function ao_one_e_integral_zero(i,k) integer, intent(in) :: i,k ao_one_e_integral_zero = .False. - if (.not.((io_ao_integrals_overlap/='None').or.is_periodic)) then + if (.not.((io_ao_integrals_overlap/='None').or.is_periodic.or.use_cgtos)) then if (ao_overlap_abs(i,k) < ao_one_e_integrals_threshold) then ao_one_e_integral_zero = .True. return diff --git a/src/ao_two_e_ints/deb_2eint_cgtos.irp.f b/src/ao_two_e_ints/deb_2eint_cgtos.irp.f deleted file mode 100644 index b871ff5a..00000000 --- a/src/ao_two_e_ints/deb_2eint_cgtos.irp.f +++ /dev/null @@ -1,153 +0,0 @@ - -! --- - -subroutine deb_ao_2eint_cgtos(i, j, k, l) - - BEGIN_DOC - ! integral of the AO basis or (ij|kl) - ! i(r1) j(r1) 1/r12 k(r2) l(r2) - END_DOC - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: 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_q1(3), iorder_q2(3) - complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3) - complex*16 :: expo1, expo2, expo3, expo4 - complex*16 :: P1_center(3), pp1 - complex*16 :: P2_center(3), pp2 - complex*16 :: Q1_center(3), qq1 - complex*16 :: Q2_center(3), qq2 - - - - dim1 = n_pt_max_integrals - - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - - if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then - - !print*, ao_prim_num(i), ao_prim_num(j), ao_prim_num(k), ao_prim_num(l) - - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) * (1.d0, 0.d0) - J_center(p) = nucl_coord(num_j,p) * (1.d0, 0.d0) - K_center(p) = nucl_coord(num_k,p) * (1.d0, 0.d0) - L_center(p) = nucl_coord(num_l,p) * (1.d0, 0.d0) - enddo - - do p = 1, ao_prim_num(i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - !print*, "expo1 = ", expo1 - !print*, "center1 = ", I_center - - do q = 1, ao_prim_num(j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - !print*, "expo2 = ", expo2 - !print*, "center2 = ", J_center - - pp1 = expo1 + expo2 - P1_center(1:3) = (expo1 * I_center(1:3) + expo2 * J_center(1:3)) / pp1 - iorder_p1(1:3) = I_power(1:3) + J_power(1:3) - - pp2 = conjg(expo1) + expo2 - P2_center(1:3) = (conjg(expo1) * I_center(1:3) + expo2 * J_center(1:3)) / pp2 - iorder_p2(1:3) = I_power(1:3) + J_power(1:3) - - do r = 1, ao_prim_num(k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - !print*, "expo3 = ", expo3 - !print*, "center3 = ", K_center - - do s = 1, ao_prim_num(l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - !print*, "expo4 = ", expo4 - !print*, "center4 = ", L_center - - qq1 = expo3 + expo4 - Q1_center(1:3) = (expo3 * K_center(1:3) + expo4 * L_center(1:3)) / qq1 - iorder_q1(1:3) = K_power(1:3) + L_power(1:3) - - qq2 = conjg(expo3) + expo4 - Q2_center(1:3) = (conjg(expo3) * K_center(1:3) + expo4 * L_center(1:3)) / qq2 - iorder_q2(1:3) = K_power(1:3) + L_power(1:3) - - call deb_cboys(P1_center, pp1, iorder_p1, Q1_center, qq1, iorder_q1) - call deb_cboys(P1_center, pp1, iorder_p1, Q2_center, qq2, iorder_q2) - call deb_cboys(P2_center, pp2, iorder_p2, Q1_center, qq1, iorder_q1) - call deb_cboys(P2_center, pp2, iorder_p2, Q2_center, qq2, iorder_q2) - call deb_cboys(conjg(P2_center), conjg(pp2), iorder_p2, Q1_center, qq1, iorder_q1) - call deb_cboys(conjg(P2_center), conjg(pp2), iorder_p2, Q2_center, qq2, iorder_q2) - call deb_cboys(conjg(P1_center), conjg(pp1), iorder_p1, Q1_center, qq1, iorder_q1) - call deb_cboys(conjg(P1_center), conjg(pp1), iorder_p1, Q2_center, qq2, iorder_q2) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif ! same centers - - return -end - -! --- - -subroutine deb_cboys(P_center, p, iorder_p, Q_center, q, iorder_q) - - - implicit none - include 'utils/constants.include.F' - - integer, intent(in) :: iorder_p(3), iorder_q(3) - complex*16, intent(in) :: P_center(3), p - complex*16, intent(in) :: Q_center(3), q - - integer :: iorder, n - complex*16 :: dist, rho - complex*16 :: int1, int2 - - complex*16, external :: crint_2 - - - dist = (P_center(1) - Q_center(1)) * (P_center(1) - Q_center(1)) & - + (P_center(2) - Q_center(2)) * (P_center(2) - Q_center(2)) & - + (P_center(3) - Q_center(3)) * (P_center(3) - Q_center(3)) - rho = dist * p * q / (p + q) - - if(abs(rho) .lt. 1d-15) return - - iorder = 2*iorder_p(1)+2*iorder_q(1) + 2*iorder_p(2)+2*iorder_q(2) + 2*iorder_p(3)+2*iorder_q(3) - n = shiftr(iorder, 1) - - !write(33,*) n, real(rho), aimag(rho) - !print*, n, real(rho), aimag(rho) - - int1 = crint_2(n, rho) - call crint_quad_12(n, rho, 1000000, int2) - - if(abs(int1 - int2) .gt. 1d-5) then - print*, ' important error found: ' - print*, p!, P_center - print*, q!, Q_center - print*, dist - print*, " n, tho = ", n, real(rho), aimag(rho) - print*, real(int1), real(int2), dabs(real(int1-int2)) - print*, aimag(int1), aimag(int2), dabs(aimag(int1-int2)) - stop - endif - -end - -! --- - diff --git a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f index 3cdb6b73..a7521b78 100644 --- a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f +++ b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f @@ -882,15 +882,13 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ complex*16, intent(in) :: Q_new(0:max_dim,3), Q_center(3), fact_q, q, q_inv integer :: i, j, nx, ny, nz, n_Ix, n_Iy, n_Iz, iorder, n_pt_tmp, n_pt_out - double precision :: tmp_mod - double precision :: ppq_re, ppq_im, ppq_mod, sq_ppq_re, sq_ppq_im complex*16 :: pq, pq_inv, pq_inv_2, p01_1, p01_2, p10_1, p10_2, ppq, sq_ppq complex*16 :: rho, dist, const complex*16 :: accu, tmp_p, tmp_q 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, external :: crint_sum !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol @@ -919,15 +917,13 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ n_Ix = 0 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)) - if(tmp_mod < thresh) cycle + tmp_p = P_new(i,1) + if(zabs(tmp_p) < 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)) - if(tmp_mod < thresh) cycle + tmp_q = tmp_p * Q_new(j,1) + if(zabs(tmp_q) < thresh) cycle !DIR$ FORCEINLINE call give_cpolynom_mult_center_x(P_center(1), Q_center(1), i, j, p, q, iorder, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, dx, nx) @@ -950,15 +946,13 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ n_Iy = 0 do i = 0, iorder_p(2) - tmp_p = P_new(i,2) - tmp_mod = dsqrt(REAL(tmp_p)*REAL(tmp_p) + AIMAG(tmp_p)*AIMAG(tmp_p)) - if(tmp_mod < thresh) cycle + tmp_p = P_new(i,2) + if(zabs(tmp_p) < thresh) cycle do j = 0, iorder_q(2) - tmp_q = tmp_p * Q_new(j,2) - tmp_mod = dsqrt(REAL(tmp_q)*REAL(tmp_q) + AIMAG(tmp_q)*AIMAG(tmp_q)) - if(tmp_mod < thresh) cycle + tmp_q = tmp_p * Q_new(j,2) + if(zabs(tmp_q) < thresh) cycle !DIR$ FORCEINLINE call give_cpolynom_mult_center_x(P_center(2), Q_center(2), i, j, p, q, iorder, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, dy, ny) @@ -982,15 +976,13 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ n_Iz = 0 do i = 0, iorder_p(3) - tmp_p = P_new(i,3) - tmp_mod = dsqrt(REAL(tmp_p)*REAL(tmp_p) + AIMAG(tmp_p)*AIMAG(tmp_p)) - if(tmp_mod < thresh) cycle + tmp_p = P_new(i,3) + if(zabs(tmp_p) < thresh) cycle do j = 0, iorder_q(3) - tmp_q = tmp_p * Q_new(j,3) - tmp_mod = dsqrt(REAL(tmp_q)*REAL(tmp_q) + AIMAG(tmp_q)*AIMAG(tmp_q)) - if(tmp_mod < thresh) cycle + tmp_q = tmp_p * Q_new(j,3) + if(zabs(tmp_q) < thresh) cycle !DIR$ FORCEINLINE call give_cpolynom_mult_center_x(P_center(3), Q_center(3), i, j, p, q, iorder, pq_inv, pq_inv_2, p10_1, p01_1, p10_2, p01_2, dz, nz) @@ -1028,6 +1020,9 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ !DIR$ FORCEINLINE call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) + if(n_pt_out == -1) then + return + endif accu = crint_sum(n_pt_out, const, d1) @@ -1056,7 +1051,6 @@ complex*16 function ERI_cgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_y, integer :: a_x_2, b_x_2, c_x_2, d_x_2, a_y_2, b_y_2, c_y_2, d_y_2, a_z_2, b_z_2, c_z_2, d_z_2 integer :: i, j, k, l, n_pt integer :: nx, ny, nz - double precision :: ppq_re, ppq_im, ppq_mod, sq_ppq_re, sq_ppq_im complex*16 :: p, q, ppq, sq_ppq, coeff, I_f ERI_cgtos = (0.d0, 0.d0) @@ -1494,7 +1488,7 @@ recursive subroutine I_x1_pol_mult_a1_cgtos(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt complex*16 :: Y(0:max_dim) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y - if( (c < 0) .or. (nd < 0) ) then + if((c < 0) .or. (nd < 0)) then nd = -1 return endif diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index ec048ab4..3734d4a0 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -46,6 +46,7 @@ double precision function ao_two_e_integral(i, j, k, l) if(use_cgtos) then ao_two_e_integral = ao_two_e_integral_cgtos(i, j, k, l) + return else if (use_only_lr) then diff --git a/src/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f deleted file mode 100644 index 99b5b36a..00000000 --- a/src/hartree_fock/deb_ao_2e_int.irp.f +++ /dev/null @@ -1,637 +0,0 @@ - -program deb_ao_2e_int - - implicit none - - !call check_ao_two_e_integral_cgtos() - !call check_crint1() - !call check_crint2() - !call check_crint3() - !call check_crint4() - call check_crint5() - !call check_crint6() - -end - -! --- - -subroutine check_ao_two_e_integral_cgtos() - - implicit none - - integer :: i, j, k, l - double precision :: acc, nrm, dif - double precision :: tmp1, tmp2 - double precision :: pw, pw0 - double precision :: t1, t2, tt - - double precision, external :: ao_two_e_integral - double precision, external :: ao_two_e_integral_cgtos - - acc = 0.d0 - nrm = 0.d0 - - pw0 = dble(ao_num**3) - pw = 0.d0 - tt = 0.d0 - do i = 1, ao_num - call wall_time(t1) - do j = 1, ao_num - do k = 1, ao_num - do l = 1, ao_num - - call deb_ao_2eint_cgtos(i, j, k, l) - - !tmp1 = ao_two_e_integral (i, j, k, l) - !tmp2 = ao_two_e_integral_cgtos(i, j, k, l) - - !print*, i, j, k, l - - !dif = abs(tmp1 - tmp2) - !!if(dif .gt. 1d-10) then - ! print*, ' error on:', i, j, k, l - ! print*, tmp1, tmp2, dif - ! !stop - !!endif - !acc += dif - !nrm += abs(tmp1) - enddo - enddo - enddo - call wall_time(t2) - tt += t2 - t1 - print*, " % done = ", 100.d0 * dble(i) / ao_num - print*, ' ellepsed time (sec) =', tt - 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 - - 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) - 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)) - - 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 - -! --- - -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 = 1 - - 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 - - n = 0 - !rho = (-6.83897018210218d0, -7.24479852507338d0) - rho = (-9.83206247355480d0, 0.445269582329036d0) - - 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 - -! --- - -subroutine check_crint4() - - implicit none - - integer :: i_test, n_test - integer :: i, seed_size, clock_time - double precision :: xr(1), x, shift - double precision :: yr(1), y - double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im - double precision :: t1, t2, t_int1, t_int2 - complex*16 :: rho - complex*16 :: int1, int2, int3 - integer, allocatable :: seed(:) - - - - 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 = 100 - shift = 15.d0 - - acc_re = 0.d0 - nrm_re = 0.d0 - acc_im = 0.d0 - nrm_im = 0.d0 - do i_test = 1, n_test - - call random_number(xr) - call random_number(yr) - - x = 1.d0 * (2.d0 * shift * xr(1) - shift) - y = 1.d0 * (2.d0 * shift * yr(1) - shift) - - rho = x + (0.d0, 1.d0) * y - - call wall_time(t1) - call zboysfun00_1(rho, int1) - call wall_time(t2) - t_int1 = t_int1 + t2 - t1 - - call wall_time(t1) - call zboysfun00_2(rho, int2) - call wall_time(t2) - t_int2 = t_int2 + t2 - t1 - - dif_re = dabs(real(int1) - real(int2)) - dif_im = dabs(aimag(int1) - aimag(int2)) - if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then - print*, ' important error found: ' - print*, " rho = ", x, y - print*, real(int1), real(int2), dif_re - print*, aimag(int1), aimag(int2), dif_im - call crint_quad_12(0, rho, 10000000, int3) - if(zabs(int1 - int3) .lt. zabs(int2 - int3)) then - print*, ' implementation 2 seems to be wrong' - else - print*, ' implementation 1 seems to be wrong' - print*, ' quad 10000000:', real(int3), aimag(int3) - call crint_quad_12(0, rho, 100000000, int3) - print*, ' quad 100000000:', real(int3), aimag(int3) - endif - !print*, ' quad:', real(int3), aimag(int3) - !stop - endif - - if((real(int1) /= real(int1)) .or. (aimag(int1) /= aimag(int1)) .or. & - (real(int2) /= real(int2)) .or. (aimag(int2) /= aimag(int2)) ) then - cycle - else - acc_re += dif_re - acc_im += dif_im - nrm_re += dabs(real(int1)) - nrm_im += dabs(aimag(int1)) - 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*, "zerf_1 wall time (sec) = ", t_int1 - print*, "zerf_2 wall time (sec) = ", t_int2 - - - deallocate(seed) - -end - -! --- - -subroutine check_crint5() - - implicit none - - integer :: i_test, n_test - integer :: i, seed_size, clock_time - integer :: n - double precision :: xr(1), yr(1), nr(1), x, shift, y - double precision :: dif1_re, dif1_im, acc1_re, acc1_im - double precision :: dif2_re, dif2_im, acc2_re, acc2_im - double precision :: nrm_re, nrm_im - double precision :: t1, t2, t_int1, t_int2 - complex*16 :: rho - complex*16 :: int1, int2, int_ref - integer, allocatable :: seed(:) - - complex*16, external :: crint_1, 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 = 100 - - acc1_re = 0.d0 - acc1_im = 0.d0 - acc2_re = 0.d0 - acc2_im = 0.d0 - nrm_re = 0.d0 - nrm_im = 0.d0 - do i_test = 1, n_test - - call random_number(xr) - call random_number(yr) - call random_number(nr) - - x = 1.d+1 * (30.d0 * xr(1) - 15.d0) - y = 1.d+1 * (30.d0 * yr(1) - 15.d0) - n = int(16.d0 * nr(1)) - - rho = x + (0.d0, 1.d0) * y - - call wall_time(t1) - int1 = crint_1(n, rho) - call wall_time(t2) - t_int1 = t_int1 + t2 - t1 - - call wall_time(t1) - int2 = crint_2(n, rho) - call wall_time(t2) - t_int2 = t_int2 + t2 - t1 - - call crint_quad_12(n, rho, 10000000, int_ref) - - dif1_re = dabs(real(int1) - real(int_ref)) - dif1_im = dabs(aimag(int1) - aimag(int_ref)) - - dif2_re = dabs(real(int2) - real(int_ref)) - dif2_im = dabs(aimag(int2) - aimag(int_ref)) - - if((dif2_re .gt. 1d-7) .or. (dif2_im .gt. 1d-7)) then - print*, ' important error found: ' - print*, " n, rho = ", n, x, y - print*, real(int1), real(int2), real(int_ref) - print*, aimag(int1), aimag(int2), aimag(int_ref) - !stop - endif - - acc1_re += dif1_re - acc1_im += dif1_im - - acc2_re += dif2_re - acc2_im += dif2_im - - nrm_re += dabs(real(int_ref)) - nrm_im += dabs(aimag(int_ref)) - enddo - - print*, "accuracy on boys_1 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc1_im / (nrm_im + 1d-15) - print*, "accuracy on boys_2 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc2_im / (nrm_im + 1d-15) - - print*, "boys_1 wall time (sec) = ", t_int1 - print*, "boys_2 wall time (sec) = ", t_int2 - - - deallocate(seed) - -end - -! --- - -subroutine check_crint6() - - implicit none - - integer :: i_test, n_test - integer :: i, seed_size, clock_time - integer :: n - double precision :: xr(1), yr(1), nr(1), x, shift, y - double precision :: dif_re, dif_im, acc_re, acc_im - double precision :: nrm_re, nrm_im - double precision :: t1, t2, t_int1, t_int2 - complex*16 :: rho - complex*16 :: int1, int2, int3 - integer, allocatable :: seed(:) - - complex*16, external :: crint_1, 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 = 100 - - acc_re = 0.d0 - acc_im = 0.d0 - nrm_re = 0.d0 - nrm_im = 0.d0 - do i_test = 1, n_test - - call random_number(xr) - call random_number(yr) - call random_number(nr) - - x = 1.d0 * (30.d0 * xr(1) - 15.d0) - y = 1.d0 * (30.d0 * yr(1) - 15.d0) - n = int(16.d0 * nr(1)) - - rho = x + (0.d0, 1.d0) * y - - call wall_time(t1) - int1 = crint_1(n, rho) - call wall_time(t2) - t_int1 = t_int1 + t2 - t1 - - call wall_time(t1) - int2 = crint_2(n, rho) - call wall_time(t2) - t_int2 = t_int2 + t2 - t1 - - dif_re = dabs(real(int1) - real(int2)) - dif_im = dabs(aimag(int1) - aimag(int2)) - - if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then - print*, ' important error found: ' - print*, " n, rho = ", n, x, y - print*, real(int1), real(int2), dif_re - print*, aimag(int1), aimag(int2), dif_im - call crint_quad_12(n, rho, 100000000, int3) - print*, ' quad 100000000:', real(int3), aimag(int3) - !print*, ' quad 100000000:', dabs(real(int1) - real(int3)), dabs(aimag(int1) - aimag(int3)) - !stop - endif - - acc_re += dif_re - acc_im += dif_im - nrm_re += dabs(real(int1)) - nrm_im += dabs(aimag(int1)) - enddo - - print*, "diff (%):", 100.d0 * acc_re / (nrm_re + 1d-15), 100.d0 * acc_im / (nrm_im + 1d-15) - - print*, "boys_1 wall time (sec) = ", t_int1 - print*, "boys_2 wall time (sec) = ", t_int2 - - - deallocate(seed) - -end - -! --- - diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index 83616969..f4a0b7de 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -21,7 +21,7 @@ complex*16 function overlap_cgaussian_x(Ae_center, Be_center, alpha, beta, power integer :: i, iorder_p complex*16 :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p - complex*16 :: Fc_integral + complex*16, external :: Fc_integral call give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_p, iorder_p, & @@ -42,6 +42,7 @@ complex*16 function overlap_cgaussian_x(Ae_center, Be_center, alpha, beta, power overlap_cgaussian_x = overlap_cgaussian_x * fact_p + return end ! --- @@ -92,7 +93,7 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow F_integral_tab(i) = Fc_integral(i, inv_sq_p) enddo - ab = alpha * beta * inv_sq_p + ab = alpha * beta * inv_sq_p * inv_sq_p arg = ab * (Ae_center(1) - Be_center(1)) & * (Ae_center(1) - Be_center(1)) diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index 075d9bbc..b6d3ac98 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -164,6 +164,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) END_DOC implicit none + complex*16, intent(in) :: a, b, xa(3), xb(3) complex*16, intent(out) :: p, k, xp(3) @@ -196,6 +197,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 + return end ! --- @@ -254,15 +256,10 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) complex*16, intent(inout) :: d(0:nb+nc) integer, intent(out) :: nd - integer :: ndtmp, ib, ic + integer :: ib, ic - if(ior(nc, nb) >= 0) then ! True if nc>=0 and nb>=0 - continue - else - return - endif + if(ior(nc, nb) < 0) return ! False if nc>=0 and nb>=0 - ndtmp = nb + nc do ic = 0, nc d(ic) = d(ic) + c(ic) * b(0) @@ -275,9 +272,8 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) enddo enddo - do nd = ndtmp, 0, -1 - if(abs(d(nd)) .lt. 1.d-15) cycle - exit + do nd = nb + nc, 0, -1 + if(d(nd) /= (0.d0, 0.d0)) exit enddo end @@ -421,11 +417,12 @@ complex*16 function Fc_integral(n, inv_sq_p) implicit none include 'constants.include.F' - integer, intent(in) :: n - complex*16, intent(in) :: inv_sq_p + integer, intent(in) :: n + complex*16, intent(in) :: inv_sq_p + complex*16 :: inv_sq_p2, inv_sq_p3, inv_sq_p4 ! (n)! - double precision :: fact + double precision, external :: fact if(n < 0) then Fc_integral = (0.d0, 0.d0) @@ -438,13 +435,29 @@ complex*16 function Fc_integral(n, inv_sq_p) return endif - if(n == 0) then + select case(n) + case(0) Fc_integral = sqpi * inv_sq_p - return - endif - - Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1)) + case(2) + Fc_integral = 0.5d0 * sqpi * inv_sq_p * inv_sq_p * inv_sq_p + case(4) + inv_sq_p2 = inv_sq_p * inv_sq_p + Fc_integral = 0.75d0 * sqpi * inv_sq_p * inv_sq_p2 * inv_sq_p2 + case(6) + inv_sq_p3 = inv_sq_p * inv_sq_p * inv_sq_p + Fc_integral = 1.875d0 * sqpi * inv_sq_p * inv_sq_p3 * inv_sq_p3 + case(8) + inv_sq_p3 = inv_sq_p * inv_sq_p * inv_sq_p + Fc_integral = 6.5625d0 * sqpi * inv_sq_p3 * inv_sq_p3 * inv_sq_p3 + case(10) + inv_sq_p2 = inv_sq_p * inv_sq_p + inv_sq_p4 = inv_sq_p2 * inv_sq_p2 + Fc_integral = 29.53125d0 * sqpi * inv_sq_p * inv_sq_p2 * inv_sq_p4 * inv_sq_p4 + case default + Fc_integral = 2.d0 * sqpi * (0.5d0 * inv_sq_p)**(n+1) * fact(n) / fact(shiftr(n, 1)) + end select + return end ! ---