From e436e22f2a51fa56fed0cfebcaf46b1c6e7cac5e Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 14:28:29 +0100 Subject: [PATCH] remove Gauss_Prod when expo = 0 --- src/ao_many_one_e_ints/grad2_jmu_manu.irp.f | 143 ++++++++++++------ .../grad_lapl_jmu_manu.irp.f | 79 ++++++---- src/ao_many_one_e_ints/listj1b.irp.f | 20 +-- 3 files changed, 156 insertions(+), 86 deletions(-) diff --git a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f index f01ed5ba..8e253d75 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f @@ -18,6 +18,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n double precision :: int_gauss, dsqpi_3_2, int_j1b double precision :: factor_ij_1s, beta_ij, center_ij_1s(3), sq_pi_3_2 double precision, allocatable :: int_fit_v(:) + double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao_with1s print*, ' providing int2_grad1u2_grad2u2_j1b2_test ...' @@ -39,48 +40,61 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, & !$OMP ao_overlap_abs,sq_pi_3_2) !$OMP DO SCHEDULE(dynamic) - 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) - do i = 1, ao_num - do j = i, ao_num - if(ao_overlap_abs(j,i) .lt. 1.d-12) then - cycle - endif - - do i_1s = 1, List_comb_thr_b3_size(j,i) + 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) + do i = 1, ao_num + do j = i, ao_num + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + cycle + endif - coef = List_comb_thr_b3_coef (i_1s,j,i) - beta = List_comb_thr_b3_expo (i_1s,j,i) - int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) - B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) - B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) - B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b3_j1b(1,j,i) + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_1_erf_x_2(i_fit) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) + if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle + int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j) + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- - do i_fit = 1, ng_fit_jast + do i_1s = 2, List_comb_thr_b3_size(j,i) + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) - expo_fit = expo_gauss_1_erf_x_2(i_fit) - !DIR$ FORCEINLINE - call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) - coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_1_erf_x_2(i_fit) + !DIR$ FORCEINLINE + call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef ! if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version - if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle - + if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle ! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & ! expo_fit, i, j, int_fit_v, n_points_final_grid) - int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss - - enddo + int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss + enddo enddo - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 1, ao_num @@ -239,9 +253,27 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ do i = 1, ao_num do j = i, ao_num - tmp = 0.d0 - do i_1s = 1, List_comb_thr_b3_size(j,i) + + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b3_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_j_mu_x_2(i_fit) + coef_fit = coef_gauss_j_mu_x_2(i_fit) + if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle + int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) + tmp += coef_fit * int_fit + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b3_size(j,i) coef = List_comb_thr_b3_coef (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i) @@ -252,23 +284,15 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) do i_fit = 1, ng_fit_jast - expo_fit = expo_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit) !DIR$ FORCEINLINE call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) ! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle - - ! --- - - int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - - tmp += coef * coef_fit * int_fit + int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + tmp += coef * coef_fit * int_fit enddo - - ! --- - enddo int2_u2_j1b2_test(j,i,ipoint) = tmp @@ -451,13 +475,34 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p do ipoint = 1, n_points_final_grid do i = 1, ao_num do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle + + if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10) cycle + r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) r(3) = final_grid_points(3,ipoint) tmp = 0.d0 - do i_1s = 1, List_comb_thr_b3_size(j,i) + + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b3_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_j_mu_1_erf(i_fit) + if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle + coef_fit = coef_gauss_j_mu_1_erf(i_fit) + int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r) + tmp += coef_fit * int_fit + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b3_size(j,i) coef = List_comb_thr_b3_coef (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i) @@ -469,9 +514,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & + (B_center(2) - r(2)) * (B_center(2) - r(2)) & + (B_center(3) - r(3)) * (B_center(3) - r(3)) - do i_fit = 1, ng_fit_jast - expo_fit = expo_gauss_j_mu_1_erf(i_fit) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.1.d-15)cycle diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f index a6a55810..5c9f81e9 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f @@ -192,9 +192,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po double precision :: coef, beta, B_center(3) double precision :: tmp double precision :: wall0, wall1 + double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot + double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b + double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao_with1s - double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b print*, ' providing v_ij_u_cst_mu_j1b_test ...' @@ -208,15 +210,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & - !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & - !$OMP SHARED (n_points_final_grid, ao_num, & - !$OMP final_grid_points, ng_fit_jast, & + !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP final_grid_points, ng_fit_jast, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & - !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & - !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & + !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP DO - !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -227,8 +228,26 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle tmp = 0.d0 - do i_1s = 1, List_comb_thr_b2_size(j,i) + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b2_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_j_mu_x(i_fit) + coef_fit = coef_gauss_j_mu_x(i_fit) + if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle + int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) + tmp += coef_fit * int_fit + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b2_size(j,i) coef = List_comb_thr_b2_coef (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) @@ -236,18 +255,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) - do i_fit = 1, ng_fit_jast - expo_fit = expo_gauss_j_mu_x(i_fit) coef_fit = coef_gauss_j_mu_x(i_fit) coeftot = coef * coef_fit if(dabs(coeftot).lt.1.d-15)cycle - double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp += coef * coef_fit * int_fit enddo enddo @@ -288,9 +303,12 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, double precision :: coef, beta, B_center(3) double precision :: tmp double precision :: wall0, wall1 + double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot + double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b + double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao_with1s - double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b + dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) provide mu_erf final_grid_points j1b_pen @@ -299,17 +317,16 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, v_ij_u_cst_mu_j1b_ng_1_test = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & - !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & - !$OMP SHARED (n_points_final_grid, ao_num, & - !$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, & - !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & - !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & - !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & + !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, & + !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & + !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP DO - !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -320,8 +337,22 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle tmp = 0.d0 - do i_1s = 1, List_comb_thr_b2_size(j,i) + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b2_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + expo_fit = expo_good_j_mu_1gauss + int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) + tmp += int_fit + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b2_size(j,i) coef = List_comb_thr_b2_coef (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) @@ -329,18 +360,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) - ! do i_fit = 1, ng_fit_jast - expo_fit = expo_good_j_mu_1gauss coef_fit = 1.d0 coeftot = coef * coef_fit if(dabs(coeftot).lt.1.d-15)cycle - double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp += coef * coef_fit * int_fit ! enddo enddo diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index e27bf723..4698cb27 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -102,11 +102,11 @@ END_PROVIDER List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) enddo - print *, ' coeff, expo & cent of list b2' - do i = 1, List_all_comb_b2_size - print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i) - print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i) - enddo + !print *, ' coeff, expo & cent of list b2' + !do i = 1, List_all_comb_b2_size + ! print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i) + ! print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i) + !enddo END_PROVIDER @@ -225,11 +225,11 @@ END_PROVIDER List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i)) enddo - print *, ' coeff, expo & cent of list b3' - do i = 1, List_all_comb_b3_size - print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i) - print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i) - enddo + !print *, ' coeff, expo & cent of list b3' + !do i = 1, List_all_comb_b3_size + ! print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i) + ! print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i) + !enddo END_PROVIDER