From 53041958a6ccbb18bfacb4bd623887dde33aa71d Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 29 Aug 2023 19:21:54 +0200 Subject: [PATCH] opt in v_ij_u_cst_mu_j1b_an: tested --- src/ao_many_one_e_ints/ao_erf_gauss.irp.f | 154 ++++++++++++++++++ .../grad_lapl_jmu_modif.irp.f | 131 ++++++++++++++- src/non_h_ints_mu/test_non_h_ints.irp.f | 44 ++++- 3 files changed, 317 insertions(+), 12 deletions(-) diff --git a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f index b1077161..823536cc 100644 --- a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f @@ -1245,3 +1245,157 @@ end subroutine NAI_pol_x2_mult_erf_ao ! --- +subroutine NAI_pol_012_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! ints(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + ! ints(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! ints(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! ints(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + ! ints(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! ints(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! ints(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3) + double precision, intent(out) :: ints(7) + + integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, m + integer :: power_A1(3), power_A2(3) + double precision :: Ai_center(3), Aj_center(3), alphai, alphaj, coef, coefi + double precision :: integral0, integral1, integral2 + + double precision, external :: NAI_pol_mult_erf_with1s + + ASSERT(beta .ge. 0.d0) + if(beta .lt. 1d-10) then + call NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) + return + endif + + ints = 0.d0 + + power_Ai(1:3) = ao_power(i_ao,1:3) + power_Aj(1:3) = ao_power(j_ao,1:3) + + Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3) + Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3) + + n_pt_in = n_pt_max_integrals + + do i = 1, ao_prim_num(i_ao) + alphai = ao_expo_ordered_transp (i,i_ao) + coefi = ao_coef_normalized_ordered_transp(i,i_ao) + + do j = 1, ao_prim_num(j_ao) + alphaj = ao_expo_ordered_transp (j,j_ao) + coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao) + + integral0 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in) + ints(1) += coef * integral0 + + do m = 1, 3 + + power_A1 = power_Ai + power_A1(m) += 1 + integral1 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A1, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in) + ints(1+m) += coef * (integral1 + Ai_center(m)*integral0) + + power_A2 = power_Ai + power_A2(m) += 2 + integral2 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A2, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in) + ints(4+m) += coef * (integral2 + Ai_center(m) * (2.d0*integral1 + Ai_center(m)*integral0)) + enddo + enddo + enddo + +end subroutine NAI_pol_012_mult_erf_ao_with1s + +! --- + +subroutine NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! int(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + ! int(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! int(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! int(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + ! int(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! int(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! int(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: mu_in, C_center(3) + double precision, intent(out) :: ints(7) + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, m + integer :: power_A1(3), power_A2(3) + double precision :: A_center(3), B_center(3), alpha, beta, coef + double precision :: integral0, integral1, integral2 + + double precision :: NAI_pol_mult_erf + + ints = 0.d0 + + num_A = ao_nucl(i_ao) + power_A(1:3) = ao_power(i_ao,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_nucl(j_ao) + power_B(1:3) = ao_power(j_ao,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + n_pt_in = n_pt_max_integrals + + do i = 1, ao_prim_num(i_ao) + alpha = ao_expo_ordered_transp(i,i_ao) + + do j = 1, ao_prim_num(j_ao) + beta = ao_expo_ordered_transp(j,j_ao) + coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) + + integral0 = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in) + ints(1) += coef * integral0 + + do m = 1, 3 + + power_A1 = power_A + power_A1(m) += 1 + integral1 = NAI_pol_mult_erf(A_center, B_center, power_A1, power_B, alpha, beta, C_center, n_pt_in, mu_in) + + ints(1+m) += coef * (integral1 + A_center(m)*integral0) + + power_A2 = power_A + power_A2(m) += 2 + integral2 = NAI_pol_mult_erf(A_center, B_center, power_A2, power_B, alpha, beta, C_center, n_pt_in, mu_in) + + ints(4+m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0)) + enddo + enddo + enddo + +end subroutine NAI_pol_012_mult_erf_ao + +! --- + diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index 9af3f9a9..24b33eb5 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -299,15 +299,12 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)] +BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC ! ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) ! - ! TODO - ! one subroutine for all integrals - ! END_DOC include 'constants.include.F' @@ -325,7 +322,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s - print*, ' providing v_ij_u_cst_mu_j1b_an ...' + print*, ' providing v_ij_u_cst_mu_j1b_an_old ...' call wall_time(wall0) provide mu_erf final_grid_points j1b_pen @@ -333,7 +330,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin ct = inv_sq_pi_2 / mu_erf - v_ij_u_cst_mu_j1b_an = 0.d0 + v_ij_u_cst_mu_j1b_an_old = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & @@ -342,7 +339,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & !$OMP final_grid_points, mu_erf, ct, & !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & - !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an) + !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an_old) !$OMP DO do ipoint = 1, n_points_final_grid @@ -413,6 +410,125 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin ! --- + v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = tmp + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = v_ij_u_cst_mu_j1b_an_old(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for v_ij_u_cst_mu_j1b_an_old', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) + ! + END_DOC + + include 'constants.include.F' + + implicit none + integer :: i, j, ipoint, i_1s + double precision :: r(3), r1_2 + double precision :: int_o + double precision :: int_c(7), int_e(7) + double precision :: coef, beta, B_center(3) + double precision :: tmp, ct + double precision :: wall0, wall1 + + double precision, external :: overlap_gauss_r12_ao_with1s + double precision, external :: NAI_pol_mult_erf_ao_with1s + + print*, ' providing v_ij_u_cst_mu_j1b_an ...' + call wall_time(wall0) + + provide mu_erf final_grid_points j1b_pen + PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent + + ct = inv_sq_pi_2 / mu_erf + + v_ij_u_cst_mu_j1b_an = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & + !$OMP r1_2, tmp, int_c, int_e, int_o) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & + !$OMP final_grid_points, mu_erf, ct, & + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & + !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an) + !$OMP DO + do ipoint = 1, n_points_final_grid + + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + r1_2 = 0.5d0 * (r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) + + do i = 1, ao_num + do j = i, ao_num + + ! --- + + coef = List_all_comb_b2_coef (1) + beta = List_all_comb_b2_expo (1) + B_center(1) = List_all_comb_b2_cent(1,1) + B_center(2) = List_all_comb_b2_cent(2,1) + B_center(3) = List_all_comb_b2_cent(3,1) + + call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c) + call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e) + + int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j) + + tmp = coef & + * ( r1_2 * (int_c(1) - int_e(1)) & + - r(1) * (int_c(2) - int_e(2)) - r(2) * (int_c(3) - int_e(3)) - r(3) * (int_c(4) - int_e(4)) & + + 0.5d0 * (int_c(5) + int_c(6) + int_c(7) - int_e(5) - int_e(6) - int_e(7)) & + - ct * int_o & + ) + + ! --- + + do i_1s = 2, List_all_comb_b2_size + + coef = List_all_comb_b2_coef (i_1s) + if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 + beta = List_all_comb_b2_expo (i_1s) + B_center(1) = List_all_comb_b2_cent(1,i_1s) + B_center(2) = List_all_comb_b2_cent(2,i_1s) + B_center(3) = List_all_comb_b2_cent(3,i_1s) + + call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c) + call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e) + + int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j) + + tmp = tmp + coef & + * ( r1_2 * (int_c(1) - int_e(1)) & + - r(1) * (int_c(2) - int_e(2)) - r(2) * (int_c(3) - int_e(3)) - r(3) * (int_c(4) - int_e(4)) & + + 0.5d0 * (int_c(5) + int_c(6) + int_c(7) - int_e(5) - int_e(6) - int_e(7)) & + - ct * int_o & + ) + + enddo + + ! --- + v_ij_u_cst_mu_j1b_an(j,i,ipoint) = tmp enddo enddo @@ -434,4 +550,3 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin END_PROVIDER ! --- - diff --git a/src/non_h_ints_mu/test_non_h_ints.irp.f b/src/non_h_ints_mu/test_non_h_ints.irp.f index aff53c2d..e7718a40 100644 --- a/src/non_h_ints_mu/test_non_h_ints.irp.f +++ b/src/non_h_ints_mu/test_non_h_ints.irp.f @@ -14,7 +14,9 @@ program test_non_h !call routine_grad_squared() !call routine_fit() - call test_ipp() + !call test_ipp() + + call test_v_ij_u_cst_mu_j1b_an() end ! --- @@ -545,9 +547,43 @@ end subroutine grad1_aos_ik_grad1_esquare ! --- - - - +subroutine test_v_ij_u_cst_mu_j1b_an() + + implicit none + integer :: i, j, ipoint + double precision :: I_old, I_new + double precision :: norm, accu, thr, diff + + PROVIDE v_ij_u_cst_mu_j1b_an_old v_ij_u_cst_mu_j1b_an + + thr = 1d-12 + norm = 0.d0 + accu = 0.d0 + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + + I_old = v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) + I_new = v_ij_u_cst_mu_j1b_an (j,i,ipoint) + + diff = dabs(I_new-I_old) + if(diff .gt. thr) then + print *, ' problem on:', j, i, ipoint + print *, ' old value :', I_old + print *, ' new value :', I_new + stop + endif + + accu += diff + norm += dabs(I_old) + enddo + enddo + enddo + + print*, ' accuracy(%) = ', 100.d0 * accu / norm + + return +end subroutine test_v_ij_u_cst_mu_j1b_an()