remove Gauss_Prod when expo = 0

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 ...'
@@ -49,7 +50,24 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
           cycle
         endif
 
-       do i_1s = 1, List_comb_thr_b3_size(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_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)
@@ -59,26 +77,22 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
   
           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
-  
 !           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
+        enddo
 
       enddo
     enddo
   enddo
-     enddo
-   enddo
-
  !$OMP END DO
  !$OMP END PARALLEL
 
@@ -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
           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
+
         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

src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f

@@ -192,10 +192,12 @@ 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, external :: overlap_gauss_r12_ao_with1s
   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
+
   print*, ' providing v_ij_u_cst_mu_j1b_test ...'
 
   dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
@@ -216,7 +218,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
  !$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, external :: overlap_gauss_r12_ao_with1s
   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
+
   dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
 
   provide mu_erf final_grid_points j1b_pen
@@ -309,7 +327,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
  !$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

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'
+  !print *, ' coeff, expo & cent of list b2'
-  do i = 1, List_all_comb_b2_size
+  !do i = 1, List_all_comb_b2_size
-    print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i)
+  !  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)
+  !  print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i)
-  enddo
+  !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'
+  !print *, ' coeff, expo & cent of list b3'
-  do i = 1, List_all_comb_b3_size
+  !do i = 1, List_all_comb_b3_size
-    print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i)
+  !  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)
+  !  print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i)
-  enddo
+  !enddo
 
 END_PROVIDER