added ao_many_one_e_ints/ bi_ortho_mos/

external/qp2-dependencies

@@ -1 +1 @@
-Subproject commit 242151e03d1d6bf042387226431d82d35845686a
+Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8

src/ao_many_one_e_ints/NEED

@@ -0,0 +1,5 @@
+ao_one_e_ints
+ao_two_e_ints
+becke_numerical_grid
+mo_one_e_ints
+dft_utils_in_r

src/ao_many_one_e_ints/README.rst

@@ -0,0 +1,25 @@
+==================
+ao_many_one_e_ints
+==================
+
+This module contains A LOT of one-electron integrals of the type 
+A_ij( r ) = \int dr' phi_i(r') w(r,r') phi_j(r') 
+where r is a point in real space. 
+
++) ao_gaus_gauss.irp.f: w(r,r') is a exp(-(r-r')^2)  , and can be multiplied by x/y/z
++) ao_erf_gauss.irp.f : w(r,r') is a exp(-(r-r')^2) erf(mu * |r-r'|)/|r-r'| , and can be multiplied by x/y/z
++) ao_erf_gauss_grad.irp.f: w(r,r') is a exp(-(r-r')^2) erf(mu * |r-r'|)/|r-r'| , and can be multiplied by x/y/z, but evaluated with also one gradient of an AO function. 
+
+Fit of a Slater function and corresponding integrals
+----------------------------------------------------
+The file fit_slat_gauss.irp.f contains many useful providers/routines to fit a Slater function with 20 gaussian. 
++) coef_fit_slat_gauss : coefficients of the gaussians to fit e^(-x)
++) expo_fit_slat_gauss : exponents of the gaussians to fit e^(-x)
+
+Integrals involving Slater functions : stg_gauss_int.irp.f
+
+Taylor expansion of full correlation factor
+-------------------------------------------
+In taylor_exp.irp.f you might find interesting integrals of the type 
+\int dr' exp( e^{-alpha |r-r|' - beta |r-r'|^2}) phi_i(r') phi_j(r') 
+evaluated as a Taylor expansion of the exponential. 

src/ao_many_one_e_ints/ao_erf_gauss_grad.irp.f

@@ -0,0 +1,150 @@
+subroutine phi_j_erf_mu_r_dxyz_phi(i,j,mu_in, C_center, dxyz_ints)
+ implicit none
+ BEGIN_DOC
+! dxyz_ints(1/2/3) = int dr phi_i(r) [erf(mu  |r - C|)/|r-C|]  d/d(x/y/z) phi_i(r)
+ END_DOC
+ integer, intent(in) :: i,j
+ double precision, intent(in) :: mu_in, C_center(3)
+ double precision, intent(out):: dxyz_ints(3)
+ integer :: num_A,power_A(3), num_b, power_B(3),power_B_tmp(3)
+ double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf,coef,thr
+ integer :: n_pt_in,l,m,mm
+ thr = 1.d-12
+ dxyz_ints = 0.d0
+ if(ao_overlap_abs(j,i).lt.thr)then
+  return
+ endif
+
+ n_pt_in = n_pt_max_integrals
+ ! j 
+ num_A = ao_nucl(j)
+ power_A(1:3)= ao_power(j,1:3)
+ A_center(1:3) = nucl_coord(num_A,1:3)
+ ! i 
+ num_B = ao_nucl(i)
+ power_B(1:3)= ao_power(i,1:3)
+ B_center(1:3) = nucl_coord(num_B,1:3)
+
+ do l=1,ao_prim_num(j)
+  alpha = ao_expo_ordered_transp(l,j)
+  do m=1,ao_prim_num(i)
+    beta = ao_expo_ordered_transp(m,i)
+    coef = ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i) 
+    if(dabs(coef).lt.thr)cycle
+    do mm = 1, 3
+     ! (d/dx phi_i ) * phi_j 
+     ! d/dx * (x - B_x)^b_x exp(-beta * (x -B_x)^2)= [b_x * (x - B_x)^(b_x - 1) - 2 beta * (x - B_x)^(b_x + 1)] exp(-beta * (x -B_x)^2)
+     !
+     ! first contribution :: b_x (x - B_x)^(b_x-1) :: integral with b_x=>b_x-1 multiplied by b_x
+     power_B_tmp = power_B
+     power_B_tmp(mm) += -1
+     contrib = NAI_pol_mult_erf(A_center,B_center,power_A,power_B_tmp,alpha,beta,C_center,n_pt_in,mu_in)  
+     dxyz_ints(mm) += contrib * dble(power_B(mm)) * coef 
+                                                  
+     ! second contribution ::  - 2 beta * (x - B_x)^(b_x + 1) :: integral with b_x=> b_x+1 multiplied by -2 * beta
+     power_B_tmp = power_B
+     power_B_tmp(mm) += 1
+     contrib = NAI_pol_mult_erf(A_center,B_center,power_A,power_B_tmp,alpha,beta,C_center,n_pt_in,mu_in)  
+     dxyz_ints(mm) += contrib * (-2.d0 * beta )  * coef 
+                                                 
+    enddo
+  enddo
+ enddo
+end
+
+
+
+
+subroutine phi_j_erf_mu_r_dxyz_phi_bis(i,j,mu_in, C_center, dxyz_ints)
+ implicit none
+ BEGIN_DOC
+! dxyz_ints(1/2/3) = int dr phi_j(r) [erf(mu  |r - C|)/|r-C|]  d/d(x/y/z) phi_i(r)
+ END_DOC
+ integer, intent(in) :: i,j
+ double precision, intent(in) :: mu_in, C_center(3)
+ double precision, intent(out):: dxyz_ints(3)
+ integer :: num_A,power_A(3), num_b, power_B(3),power_B_tmp(3)
+ double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf
+ double precision :: thr, coef
+ integer :: n_pt_in,l,m,mm,kk
+ thr = 1.d-12
+ dxyz_ints = 0.d0
+ if(ao_overlap_abs(j,i).lt.thr)then
+  return
+ endif
+
+ n_pt_in = n_pt_max_integrals
+ ! j == A 
+ num_A = ao_nucl(j)
+ power_A(1:3)= ao_power(j,1:3)
+ A_center(1:3) = nucl_coord(num_A,1:3)
+ ! i == B
+ num_B = ao_nucl(i)
+ power_B(1:3)= ao_power(i,1:3)
+ B_center(1:3) = nucl_coord(num_B,1:3)
+
+ dxyz_ints = 0.d0
+ do l=1,ao_prim_num(j)
+  alpha = ao_expo_ordered_transp(l,j)
+  do m=1,ao_prim_num(i)
+    beta = ao_expo_ordered_transp(m,i)
+    do kk = 1, 2 ! loop over the extra terms induced by the d/dx/y/z * AO(i)
+     do mm = 1, 3
+      power_B_tmp = power_B
+      power_B_tmp(mm) = power_ord_grad_transp(kk,mm,i)
+      coef = ao_coef_normalized_ordered_transp(l,j) * ao_coef_ord_grad_transp(kk,mm,m,i) 
+      if(dabs(coef).lt.thr)cycle
+      contrib = NAI_pol_mult_erf(A_center,B_center,power_A,power_B_tmp,alpha,beta,C_center,n_pt_in,mu_in)  
+      dxyz_ints(mm) += contrib * coef 
+     enddo
+    enddo
+  enddo
+ enddo
+end
+
+subroutine phi_j_erf_mu_r_xyz_dxyz_phi(i,j,mu_in, C_center, dxyz_ints)
+ implicit none
+ BEGIN_DOC
+! dxyz_ints(1/2/3) = int dr phi_j(r) x/y/z [erf(mu  |r - C|)/|r-C|]  d/d(x/y/z) phi_i(r)
+ END_DOC
+ integer, intent(in) :: i,j
+ double precision, intent(in) :: mu_in, C_center(3)
+ double precision, intent(out):: dxyz_ints(3)
+ integer :: num_A,power_A(3), num_b, power_B(3),power_B_tmp(3)
+ double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf
+ double precision :: thr, coef
+ integer :: n_pt_in,l,m,mm,kk
+ thr = 1.d-12
+ dxyz_ints = 0.d0
+ if(ao_overlap_abs(j,i).lt.thr)then
+  return
+ endif
+
+ n_pt_in = n_pt_max_integrals
+ ! j == A 
+ num_A = ao_nucl(j)
+ power_A(1:3)= ao_power(j,1:3)
+ A_center(1:3) = nucl_coord(num_A,1:3)
+ ! i == B
+ num_B = ao_nucl(i)
+ power_B(1:3)= ao_power(i,1:3)
+ B_center(1:3) = nucl_coord(num_B,1:3)
+
+ dxyz_ints = 0.d0
+ do l=1,ao_prim_num(j)
+  alpha = ao_expo_ordered_transp(l,j)
+  do m=1,ao_prim_num(i)
+    beta = ao_expo_ordered_transp(m,i)
+    do kk = 1, 4 ! loop over the extra terms induced by the x/y/z * d dx/y/z AO(i)
+     do mm = 1, 3
+      power_B_tmp = power_B
+      power_B_tmp(mm) = power_ord_xyz_grad_transp(kk,mm,i)
+      coef = ao_coef_normalized_ordered_transp(l,j) * ao_coef_ord_xyz_grad_transp(kk,mm,m,i) 
+      if(dabs(coef).lt.thr)cycle
+      contrib = NAI_pol_mult_erf(A_center,B_center,power_A,power_B_tmp,alpha,beta,C_center,n_pt_in,mu_in)  
+      dxyz_ints(mm) += contrib * coef 
+     enddo
+    enddo
+  enddo
+ enddo
+end

src/ao_many_one_e_ints/ao_gaus_gauss.irp.f

@@ -0,0 +1,426 @@
+! ---
+
+subroutine overlap_gauss_xyz_r12_ao(D_center,delta,i,j,gauss_ints)
+
+ implicit none
+ BEGIN_DOC
+! gauss_ints(m) = \int dr AO_i(r) AO_j(r) x/y/z e^{-delta |r-D_center|^2}
+!
+! with m == 1 ==> x, m == 2 ==> y, m == 3 ==> z
+ END_DOC
+ integer, intent(in) :: i,j
+ double precision, intent(in)  :: D_center(3), delta
+ double precision, intent(out) :: gauss_ints(3)
+
+ integer :: num_a,num_b,power_A(3), power_B(3),l,k,m
+ double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,gauss_ints_tmp(3)
+ gauss_ints = 0.d0
+ if(ao_overlap_abs(j,i).lt.1.d-12)then
+  return
+ endif
+ num_A = ao_nucl(i)
+ power_A(1:3)= ao_power(i,1:3)
+ A_center(1:3) = nucl_coord(num_A,1:3)
+ num_B = ao_nucl(j)
+ power_B(1:3)= ao_power(j,1:3)
+ B_center(1:3) = nucl_coord(num_B,1:3)
+ do l=1,ao_prim_num(i)
+  alpha = ao_expo_ordered_transp(l,i)
+  do k=1,ao_prim_num(j)
+   beta = ao_expo_ordered_transp(k,j)
+   call overlap_gauss_xyz_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,gauss_ints_tmp)
+   do m = 1, 3
+    gauss_ints(m) += gauss_ints_tmp(m) *  ao_coef_normalized_ordered_transp(l,i)             &
+                                       *  ao_coef_normalized_ordered_transp(k,j)
+   enddo
+  enddo
+ enddo
+
+end
+
+
+
+double precision function overlap_gauss_xyz_r12_ao_specific(D_center,delta,i,j,mx)
+ implicit none
+ BEGIN_DOC
+! \int dr AO_i(r) AO_j(r) x/y/z e^{-delta |r-D_center|^2}
+!
+! with mx == 1 ==> x, mx == 2 ==> y, mx == 3 ==> z
+ END_DOC
+ integer, intent(in) :: i,j,mx
+ double precision, intent(in)  :: D_center(3), delta
+
+ integer :: num_a,num_b,power_A(3), power_B(3),l,k
+ double precision :: gauss_int
+ double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta
+ double precision :: overlap_gauss_xyz_r12_specific
+ overlap_gauss_xyz_r12_ao_specific = 0.d0
+ if(ao_overlap_abs(j,i).lt.1.d-12)then
+  return
+ endif
+ num_A = ao_nucl(i)
+ power_A(1:3)= ao_power(i,1:3)
+ A_center(1:3) = nucl_coord(num_A,1:3)
+ num_B = ao_nucl(j)
+ power_B(1:3)= ao_power(j,1:3)
+ B_center(1:3) = nucl_coord(num_B,1:3)
+ do l=1,ao_prim_num(i)
+  alpha = ao_expo_ordered_transp(l,i)
+  do k=1,ao_prim_num(j)
+   beta = ao_expo_ordered_transp(k,j)
+   gauss_int = overlap_gauss_xyz_r12_specific(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,mx)
+   overlap_gauss_xyz_r12_ao_specific = gauss_int *  ao_coef_normalized_ordered_transp(l,i)             &
+                                                 *  ao_coef_normalized_ordered_transp(k,j)
+  enddo
+ enddo
+end
+
+
+subroutine overlap_gauss_r12_all_ao(D_center,delta,aos_ints)
+ implicit none
+ double precision, intent(in) :: D_center(3), delta
+ double precision, intent(out):: aos_ints(ao_num,ao_num)
+
+ integer :: num_a,num_b,power_A(3), power_B(3),l,k,i,j
+ double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,analytical_j
+ aos_ints = 0.d0
+ do i = 1, ao_num
+  do j = 1, ao_num
+   if(ao_overlap_abs(j,i).lt.1.d-12)cycle
+   num_A = ao_nucl(i)
+   power_A(1:3)= ao_power(i,1:3)
+   A_center(1:3) = nucl_coord(num_A,1:3)
+   num_B = ao_nucl(j)
+   power_B(1:3)= ao_power(j,1:3)
+   B_center(1:3) = nucl_coord(num_B,1:3)
+   do l=1,ao_prim_num(i)
+    alpha = ao_expo_ordered_transp(l,i)
+    do k=1,ao_prim_num(j)
+     beta = ao_expo_ordered_transp(k,j)
+     analytical_j = overlap_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta)
+     aos_ints(j,i) += analytical_j *  ao_coef_normalized_ordered_transp(l,i)             &
+                                   *  ao_coef_normalized_ordered_transp(k,j)
+    enddo
+   enddo
+  enddo
+ enddo
+end
+
+! ---
+
+! TODO :: PUT CYCLES IN LOOPS
+double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
+
+  BEGIN_DOC
+  ! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2}
+  END_DOC
+
+  implicit none
+  integer,          intent(in) :: i, j
+  double precision, intent(in) :: D_center(3), delta
+
+  integer                      :: power_A(3), power_B(3), l, k
+  double precision             :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j
+
+  double precision, external   :: overlap_gauss_r12
+
+  overlap_gauss_r12_ao = 0.d0
+
+  if(ao_overlap_abs(j,i).lt.1.d-12) then
+    return
+  endif
+
+  power_A(1:3) = ao_power(i,1:3)
+  power_B(1:3) = ao_power(j,1:3)
+
+  A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
+  B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
+
+  do l = 1, ao_prim_num(i)
+    alpha = ao_expo_ordered_transp           (l,i)
+    coef1 = ao_coef_normalized_ordered_transp(l,i)
+
+    do k = 1, ao_prim_num(j)
+      beta = ao_expo_ordered_transp(k,j)
+      coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
+
+      if(dabs(coef) .lt. 1d-12) cycle
+
+      analytical_j = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
+
+      overlap_gauss_r12_ao += coef * analytical_j
+    enddo
+  enddo
+
+end function overlap_gauss_r12_ao
+
+! --
+
+double precision function overlap_abs_gauss_r12_ao(D_center, delta, i, j)
+
+  BEGIN_DOC
+  ! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2}
+  END_DOC
+
+  implicit none
+  integer,          intent(in) :: i, j
+  double precision, intent(in) :: D_center(3), delta
+
+  integer                      :: power_A(3), power_B(3), l, k
+  double precision             :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j
+
+  double precision, external   :: overlap_abs_gauss_r12
+
+  overlap_abs_gauss_r12_ao = 0.d0
+
+  if(ao_overlap_abs(j,i).lt.1.d-12) then
+    return
+  endif
+
+  power_A(1:3) = ao_power(i,1:3)
+  power_B(1:3) = ao_power(j,1:3)
+
+  A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
+  B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
+
+  do l = 1, ao_prim_num(i)
+    alpha = ao_expo_ordered_transp           (l,i)
+    coef1 = ao_coef_normalized_ordered_transp(l,i)
+
+    do k = 1, ao_prim_num(j)
+      beta = ao_expo_ordered_transp(k,j)
+      coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
+
+      if(dabs(coef) .lt. 1d-12) cycle
+
+      analytical_j = overlap_abs_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
+
+      overlap_abs_gauss_r12_ao += dabs(coef * analytical_j)
+    enddo
+  enddo
+
+end function overlap_gauss_r12_ao
+
+! --
+
+subroutine overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
+
+  BEGIN_DOC
+  !
+  ! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2}
+  !
+  ! n_points: nb of integrals <= min(LD_D, LD_resv)
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, j, LD_D, LD_resv, n_points
+  double precision, intent(in)  :: D_center(LD_D,3), delta
+  double precision, intent(out) :: resv(LD_resv)
+
+  integer                       :: ipoint
+  integer                       :: power_A(3), power_B(3), l, k
+  double precision              :: A_center(3), B_center(3), alpha, beta, coef, coef1 
+  double precision, allocatable :: analytical_j(:)
+
+  resv(:) = 0.d0
+  if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+    return
+  endif
+
+  power_A(1:3) = ao_power(i,1:3)
+  power_B(1:3) = ao_power(j,1:3)
+
+  A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
+  B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
+
+  allocate(analytical_j(n_points))
+
+  do l = 1, ao_prim_num(i)
+    alpha = ao_expo_ordered_transp           (l,i)
+    coef1 = ao_coef_normalized_ordered_transp(l,i)
+
+    do k = 1, ao_prim_num(j)
+      beta = ao_expo_ordered_transp(k,j)
+      coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
+
+      if(dabs(coef) .lt. 1d-12) cycle
+
+      call overlap_gauss_r12_v(D_center, LD_D, delta, A_center, B_center, power_A, power_B, alpha, beta, analytical_j, n_points, n_points)
+
+      do ipoint = 1, n_points
+        resv(ipoint) = resv(ipoint) + coef * analytical_j(ipoint)
+      enddo
+
+    enddo
+  enddo
+
+  deallocate(analytical_j)
+
+end subroutine overlap_gauss_r12_ao_v
+
+! ---
+
+double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, delta, i, j)
+
+  BEGIN_DOC
+  !
+  ! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2}
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in) :: i, j
+  double precision, intent(in) :: B_center(3), beta, D_center(3), delta
+
+  integer                      :: power_A1(3), power_A2(3), l, k
+  double precision             :: A1_center(3), A2_center(3), alpha1, alpha2, coef1, coef12, analytical_j
+  double precision             :: G_center(3), gama, fact_g, gama_inv
+
+  double precision, external   :: overlap_gauss_r12, overlap_gauss_r12_ao
+
+  if(beta .lt. 1d-10) then
+    overlap_gauss_r12_ao_with1s = overlap_gauss_r12_ao(D_center, delta, i, j)
+    return
+  endif
+
+  overlap_gauss_r12_ao_with1s = 0.d0
+
+  if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+    return
+  endif
+
+  ! e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} = fact_g e^{-gama |r - G|^2}
+
+  gama        = beta + delta
+  gama_inv    = 1.d0 / gama
+  G_center(1) = (beta * B_center(1) + delta * D_center(1)) * gama_inv
+  G_center(2) = (beta * B_center(2) + delta * D_center(2)) * gama_inv
+  G_center(3) = (beta * B_center(3) + delta * D_center(3)) * gama_inv
+  fact_g      = beta * delta * gama_inv * ( (B_center(1) - D_center(1)) * (B_center(1) - D_center(1)) &
+                                          + (B_center(2) - D_center(2)) * (B_center(2) - D_center(2)) &
+                                          + (B_center(3) - D_center(3)) * (B_center(3) - D_center(3)) )
+  if(fact_g .gt. 10d0) return
+  fact_g = dexp(-fact_g)
+
+  ! ---
+
+  power_A1(1:3) = ao_power(i,1:3)
+  power_A2(1:3) = ao_power(j,1:3)
+
+  A1_center(1:3) = nucl_coord(ao_nucl(i),1:3)
+  A2_center(1:3) = nucl_coord(ao_nucl(j),1:3)
+
+  do l = 1, ao_prim_num(i)
+    alpha1 = ao_expo_ordered_transp                    (l,i)
+    coef1  = fact_g * ao_coef_normalized_ordered_transp(l,i)
+    if(dabs(coef1) .lt. 1d-12) cycle
+
+    do k = 1, ao_prim_num(j)
+      alpha2 = ao_expo_ordered_transp                   (k,j)
+      coef12 = coef1 * ao_coef_normalized_ordered_transp(k,j)
+      if(dabs(coef12) .lt. 1d-12) cycle
+
+      analytical_j = overlap_gauss_r12(G_center, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2)
+
+      overlap_gauss_r12_ao_with1s += coef12 * analytical_j
+    enddo
+  enddo
+
+end function overlap_gauss_r12_ao_with1s
+
+! ---
+
+subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
+
+  BEGIN_DOC
+  !
+  ! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2}
+  ! using an array of D_centers.
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, j, n_points, LD_D, LD_resv
+  double precision, intent(in)  :: B_center(3), beta, D_center(LD_D,3), delta
+  double precision, intent(out) :: resv(LD_resv)
+
+  integer                       :: ipoint
+  integer                       :: power_A1(3), power_A2(3), l, k
+  double precision              :: A1_center(3), A2_center(3), alpha1, alpha2, coef1
+  double precision              :: coef12, coef12f
+  double precision              :: gama, gama_inv
+  double precision              :: bg, dg, bdg
+  double precision, allocatable :: fact_g(:), G_center(:,:), analytical_j(:)
+
+  if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+    return
+  endif
+
+  ASSERT(beta .gt. 0.d0)
+
+  if(beta .lt. 1d-10) then
+    call overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_points)
+    return
+  endif
+
+  resv(:) = 0.d0
+
+  ! e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} = fact_g e^{-gama |r - G|^2}
+
+  gama     = beta + delta
+  gama_inv = 1.d0 / gama
+
+  power_A1(1:3) = ao_power(i,1:3)
+  power_A2(1:3) = ao_power(j,1:3)
+
+  A1_center(1:3) = nucl_coord(ao_nucl(i),1:3)
+  A2_center(1:3) = nucl_coord(ao_nucl(j),1:3)
+
+  allocate(fact_g(n_points), G_center(n_points,3), analytical_j(n_points))
+
+  bg  = beta  * gama_inv
+  dg  = delta * gama_inv
+  bdg = bg * delta 
+
+  do ipoint = 1, n_points
+
+    G_center(ipoint,1) = bg * B_center(1) + dg * D_center(ipoint,1)
+    G_center(ipoint,2) = bg * B_center(2) + dg * D_center(ipoint,2)
+    G_center(ipoint,3) = bg * B_center(3) + dg * D_center(ipoint,3)
+    fact_g(ipoint) = bdg * ( (B_center(1) - D_center(ipoint,1)) * (B_center(1) - D_center(ipoint,1)) &
+                           + (B_center(2) - D_center(ipoint,2)) * (B_center(2) - D_center(ipoint,2)) &
+                           + (B_center(3) - D_center(ipoint,3)) * (B_center(3) - D_center(ipoint,3)) )
+
+    if(fact_g(ipoint) < 10d0) then
+      fact_g(ipoint) = dexp(-fact_g(ipoint))
+    else
+      fact_g(ipoint) = 0.d0
+    endif
+
+  enddo
+
+  do l = 1, ao_prim_num(i)
+    alpha1 = ao_expo_ordered_transp           (l,i)
+    coef1  = ao_coef_normalized_ordered_transp(l,i)
+
+    do k = 1, ao_prim_num(j)
+      alpha2 = ao_expo_ordered_transp                   (k,j)
+      coef12 = coef1 * ao_coef_normalized_ordered_transp(k,j)
+      if(dabs(coef12) .lt. 1d-12) cycle
+
+      call overlap_gauss_r12_v(G_center, n_points, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, analytical_j, n_points, n_points)
+
+      do ipoint = 1, n_points
+        coef12f = coef12 * fact_g(ipoint)
+        resv(ipoint) += coef12f * analytical_j(ipoint)
+      enddo
+    enddo
+  enddo
+
+  deallocate(fact_g, G_center, analytical_j)
+
+end subroutine overlap_gauss_r12_ao_with1s_v
+
+! ---
+

src/ao_many_one_e_ints/fit_slat_gauss.irp.f

@@ -0,0 +1,94 @@
+ BEGIN_PROVIDER [integer, n_max_fit_slat]
+ implicit none
+ BEGIN_DOC
+! number of gaussian to fit exp(-x)
+!
+! I took 20 gaussians from the program bassto.f
+ END_DOC
+ n_max_fit_slat = 20
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, coef_fit_slat_gauss, (n_max_fit_slat)]
+&BEGIN_PROVIDER [double precision, expo_fit_slat_gauss, (n_max_fit_slat)]
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+ ! fit the exp(-x) as 
+ !
+ ! \sum_{i = 1, n_max_fit_slat} coef_fit_slat_gauss(i) * exp(-expo_fit_slat_gauss(i) * x**2)
+ !
+ ! The coefficient are taken from the program bassto.f
+ END_DOC
+
+
+      expo_fit_slat_gauss(01)=30573.77073000000
+      coef_fit_slat_gauss(01)=0.00338925525
+      expo_fit_slat_gauss(02)=5608.45238100000
+      coef_fit_slat_gauss(02)=0.00536433869
+      expo_fit_slat_gauss(03)=1570.95673400000
+      coef_fit_slat_gauss(03)=0.00818702846
+      expo_fit_slat_gauss(04)=541.39785110000
+      coef_fit_slat_gauss(04)=0.01202047655
+      expo_fit_slat_gauss(05)=212.43469630000
+      coef_fit_slat_gauss(05)=0.01711289568
+      expo_fit_slat_gauss(06)=91.31444574000
+      coef_fit_slat_gauss(06)=0.02376001022
+      expo_fit_slat_gauss(07)=42.04087246000
+      coef_fit_slat_gauss(07)=0.03229121736
+      expo_fit_slat_gauss(08)=20.43200443000
+      coef_fit_slat_gauss(08)=0.04303646818
+      expo_fit_slat_gauss(09)=10.37775161000
+      coef_fit_slat_gauss(09)=0.05624657578
+      expo_fit_slat_gauss(10)=5.46880754500
+      coef_fit_slat_gauss(10)=0.07192311571
+      expo_fit_slat_gauss(11)=2.97373529200
+      coef_fit_slat_gauss(11)=0.08949389001
+      expo_fit_slat_gauss(12)=1.66144190200
+      coef_fit_slat_gauss(12)=0.10727599240
+      expo_fit_slat_gauss(13)=0.95052560820
+      coef_fit_slat_gauss(13)=0.12178961750
+      expo_fit_slat_gauss(14)=0.55528683970
+      coef_fit_slat_gauss(14)=0.12740141870
+      expo_fit_slat_gauss(15)=0.33043360020
+      coef_fit_slat_gauss(15)=0.11759168160
+      expo_fit_slat_gauss(16)=0.19982303230
+      coef_fit_slat_gauss(16)=0.08953504394
+      expo_fit_slat_gauss(17)=0.12246840760
+      coef_fit_slat_gauss(17)=0.05066721317
+      expo_fit_slat_gauss(18)=0.07575825322
+      coef_fit_slat_gauss(18)=0.01806363869
+      expo_fit_slat_gauss(19)=0.04690146243
+      coef_fit_slat_gauss(19)=0.00305632563
+      expo_fit_slat_gauss(20)=0.02834749861
+      coef_fit_slat_gauss(20)=0.00013317513
+
+
+
+END_PROVIDER 
+
+double precision function slater_fit_gam(x,gam)
+ implicit none
+ double precision, intent(in) :: x,gam
+ BEGIN_DOC
+! fit of the function exp(-gam * x) with gaussian functions 
+ END_DOC
+ integer :: i
+ slater_fit_gam = 0.d0
+ do i = 1, n_max_fit_slat
+  slater_fit_gam += coef_fit_slat_gauss(i) * dexp(-expo_fit_slat_gauss(i) * gam * gam * x * x)
+ enddo
+end
+
+subroutine expo_fit_slater_gam(gam,expos)
+ implicit none
+ BEGIN_DOC
+! returns the array of the exponents of the gaussians to fit exp(-gam*x)
+ END_DOC
+ double precision, intent(in)  :: gam
+ double precision, intent(out) :: expos(n_max_fit_slat)
+ integer :: i
+ do i = 1, n_max_fit_slat
+  expos(i) = expo_fit_slat_gauss(i) * gam * gam
+ enddo
+end
+

src/ao_many_one_e_ints/grad2_jmu_manu.irp.f

@@ -0,0 +1,517 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), expo_fit, coef_fit
+  double precision              :: coef, beta, B_center(3)
+  double precision              :: tmp
+  double precision              :: wall0, wall1
+  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_with1s
+
+  print*, ' providing int2_grad1u2_grad2u2_j1b2_test ...'
+
+  sq_pi_3_2 = (dacos(-1.d0))**(1.5d0)
+
+  provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef
+  call wall_time(wall0)
+
+  int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                                              &
+     !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,                                     &
+     !$OMP          coef_fit, expo_fit, int_fit_v, tmp,int_gauss,int_j1b,factor_ij_1s,beta_ij,center_ij_1s) &
+     !$OMP SHARED  (n_points_final_grid, ao_num, final_grid_points,List_comb_thr_b3_size,                   &
+     !$OMP          final_grid_points_transp, ng_fit_jast,                                                  &
+     !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,                                             &
+     !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,                                           &
+     !$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)
+
+         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)
+  
+         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
+
+   !$OMP END DO
+   !$OMP END PARALLEL
+
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      do j = 1, i-1
+        int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao_num, n_points_final_grid)]
+!
+!  BEGIN_DOC
+!  !
+!  ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
+!  !
+!  END_DOC
+!
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), expo_fit, coef_fit
+  double precision              :: coef, beta, B_center(3)
+  double precision              :: tmp
+  double precision              :: wall0, wall1
+
+  double precision, allocatable :: int_fit_v(:),big_array(:,:,:)
+  double precision, external    :: overlap_gauss_r12_ao_with1s
+
+  print*, ' providing int2_grad1u2_grad2u2_j1b2_test_v ...'
+
+  provide mu_erf final_grid_points_transp j1b_pen
+  call wall_time(wall0)
+
+ double precision :: int_j1b
+ big_array(:,:,:) = 0.d0
+ allocate(big_array(n_points_final_grid,ao_num, ao_num))
+ !$OMP PARALLEL DEFAULT (NONE)                                       &
+     !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,&
+     !$OMP          coef_fit, expo_fit, int_fit_v, tmp,int_j1b)                &
+     !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b3_size,&
+     !$OMP          final_grid_points_transp, ng_fit_jast,               &
+     !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,      &
+     !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,    &
+     !$OMP          List_comb_thr_b3_cent, big_array,&
+     !$OMP          ao_abs_comb_b3_j1b,ao_overlap_abs)
+!
+ allocate(int_fit_v(n_points_final_grid))
+ !$OMP DO SCHEDULE(dynamic)
+ 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)
+
+         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)
+!         if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle
+         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)
+
+       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) * coef
+
+         call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, size(final_grid_points_transp,1),&
+               expo_fit, i, j, int_fit_v, size(int_fit_v,1),n_points_final_grid)
+
+         do ipoint = 1, n_points_final_grid
+           big_array(ipoint,j,i) += coef_fit * int_fit_v(ipoint)
+         enddo
+
+       enddo
+
+     enddo
+   enddo
+ enddo
+ !$OMP END DO
+ deallocate(int_fit_v)
+ !$OMP END PARALLEL
+ do i = 1, ao_num
+   do j = i, ao_num
+    do ipoint = 1, n_points_final_grid
+     int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,j,i)
+    enddo
+   enddo
+  enddo
+
+  do ipoint = 1, n_points_final_grid
+    do i = 2, ao_num
+      do j = 1, i-1
+        int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,i,j)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_v', wall1 - wall0
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), int_fit, expo_fit, coef_fit
+  double precision              :: coef, beta, B_center(3), tmp
+  double precision              :: wall0, wall1,int_j1b
+
+  double precision, external    :: overlap_gauss_r12_ao
+  double precision, external    :: overlap_gauss_r12_ao_with1s
+  double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),sq_pi_3_2
+
+  print*, ' providing int2_u2_j1b2_test ...'
+
+  sq_pi_3_2 = (dacos(-1.d0))**(1.5d0)
+
+  provide mu_erf final_grid_points j1b_pen
+  call wall_time(wall0)
+
+  int2_u2_j1b2_test = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, tmp, int_j1b,factor_ij_1s,beta_ij,center_ij_1s)          & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b3_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2,           &
+ !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2,       & 
+ !$OMP          List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b)
+ !$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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+
+        tmp = 0.d0
+        do i_1s = 1, 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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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)
+
+          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
+      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
+        int2_u2_j1b2_test(j,i,ipoint) = int2_u2_j1b2_test(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_u2_j1b2_test', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n_points_final_grid, 3)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint, i_1s, i_fit
+  double precision :: r(3), int_fit(3), expo_fit, coef_fit
+  double precision :: coef, beta, B_center(3), dist
+  double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp
+  double precision :: tmp_x, tmp_y, tmp_z, int_j1b
+  double precision :: wall0, wall1, sq_pi_3_2,sq_alpha
+
+  print*, ' providing int2_u_grad1u_x_j1b2_test ...'
+
+  sq_pi_3_2 = dacos(-1.D0)**(1.d0)
+  provide mu_erf final_grid_points j1b_pen
+  call wall_time(wall0)
+
+  int2_u_grad1u_x_j1b2_test = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, alpha_1s, dist,        &
+ !$OMP          alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp,     & 
+ !$OMP          tmp_x, tmp_y, tmp_z,int_j1b,sq_alpha)                        & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b3_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
+ !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,       & 
+ !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2)
+ !$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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp_x = 0.d0
+        tmp_y = 0.d0
+        tmp_z = 0.d0
+        do i_1s = 1, 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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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)
+          do i_fit = 1, ng_fit_jast
+    
+            expo_fit = expo_gauss_j_mu_1_erf(i_fit)
+            coef_fit = coef_gauss_j_mu_1_erf(i_fit)
+    
+            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)) 
+
+            alpha_1s     = beta + expo_fit
+            alpha_1s_inv = 1.d0 / alpha_1s 
+
+            centr_1s(1)  = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
+            centr_1s(2)  = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
+            centr_1s(3)  = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
+
+            expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist 
+            coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
+            sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
+!            if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle ! old version
+            if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. 1d-10) cycle
+            
+            call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
+
+            tmp_x += coef_tmp * int_fit(1)
+            tmp_y += coef_tmp * int_fit(2)
+            tmp_z += coef_tmp * int_fit(3)
+          enddo
+
+          ! ---
+
+        enddo
+
+        int2_u_grad1u_x_j1b2_test(j,i,ipoint,1) = tmp_x
+        int2_u_grad1u_x_j1b2_test(j,i,ipoint,2) = tmp_y
+        int2_u_grad1u_x_j1b2_test(j,i,ipoint,3) = tmp_z
+      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
+        int2_u_grad1u_x_j1b2_test(j,i,ipoint,1) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,1)
+        int2_u_grad1u_x_j1b2_test(j,i,ipoint,2) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,2)
+        int2_u_grad1u_x_j1b2_test(j,i,ipoint,3) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,3)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_u_grad1u_x_j1b2_test', wall1 - wall0
+
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu]
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), int_fit, expo_fit, coef_fit, coef_tmp
+  double precision              :: coef, beta, B_center(3), dist
+  double precision              :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
+  double precision              :: wall0, wall1
+  double precision, external    :: NAI_pol_mult_erf_ao_with1s
+  double precision :: j12_mu_r12,int_j1b
+  double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
+  double precision :: beta_ij,center_ij_1s(3),factor_ij_1s
+
+  print*, ' providing int2_u_grad1u_j1b2_test ...'
+
+  dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
+
+  provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent
+  call wall_time(wall0)
+
+
+  int2_u_grad1u_j1b2_test = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist,   &
+ !$OMP          beta_ij,center_ij_1s,factor_ij_1s,               &
+ !$OMP          int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp)     &
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
+ !$OMP          final_grid_points, ng_fit_jast,                  &
+ !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
+ !$OMP          ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2,   &
+ !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,  ao_abs_comb_b3_j1b,     &
+ !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test)
+ !$OMP DO
+  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)
+
+          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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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)
+          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
+            coef_fit = coef_gauss_j_mu_1_erf(i_fit)
+
+            alpha_1s     = beta + expo_fit
+            alpha_1s_inv = 1.d0 / alpha_1s
+            centr_1s(1)  = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
+            centr_1s(2)  = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
+            centr_1s(3)  = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
+
+            expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
+            if(expo_coef_1s .gt. 20.d0) cycle
+            coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
+            if(dabs(coef_tmp) .lt. 1d-08) cycle
+
+            int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s,  1.d+9, r)
+
+            tmp += coef_tmp * int_fit
+          enddo
+        enddo
+
+        int2_u_grad1u_j1b2_test(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
+        int2_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0
+
+END_PROVIDER
+
+! ---

src/ao_many_one_e_ints/grad2_jmu_modif.irp.f

@@ -0,0 +1,420 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), int_fit, expo_fit, coef_fit
+  double precision              :: coef, beta, B_center(3)
+  double precision              :: tmp
+  double precision              :: wall0, wall1
+
+  double precision, external    :: overlap_gauss_r12_ao
+  double precision, external    :: overlap_gauss_r12_ao_with1s
+
+  print*, ' providing int2_grad1u2_grad2u2_j1b2 ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  int2_grad1u2_grad2u2_j1b2 = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, tmp)                   & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,         &
+ !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,       & 
+ !$OMP          List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2)
+ !$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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp = 0.d0
+        do i_fit = 1, ng_fit_jast
+
+          expo_fit = expo_gauss_1_erf_x_2(i_fit)
+          coef_fit = coef_gauss_1_erf_x_2(i_fit)
+
+          ! ---
+
+          int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
+          tmp += -0.25d0 * coef_fit * int_fit
+!          if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle
+
+          ! ---
+
+          do i_1s = 2, List_all_comb_b3_size
+
+            coef        = List_all_comb_b3_coef  (i_1s)
+            beta        = List_all_comb_b3_expo  (i_1s)
+            B_center(1) = List_all_comb_b3_cent(1,i_1s)
+            B_center(2) = List_all_comb_b3_cent(2,i_1s)
+            B_center(3) = List_all_comb_b3_cent(3,i_1s)
+
+            int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
+
+            tmp += -0.25d0 * coef * coef_fit * int_fit
+          enddo
+
+          ! ---
+
+        enddo
+
+        int2_grad1u2_grad2u2_j1b2(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
+        int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_grad1u2_grad2u2_j1b2 =', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), int_fit, expo_fit, coef_fit
+  double precision              :: coef, beta, B_center(3), tmp
+  double precision              :: wall0, wall1
+
+  double precision, external    :: overlap_gauss_r12_ao
+  double precision, external    :: overlap_gauss_r12_ao_with1s
+
+  print*, ' providing int2_u2_j1b2 ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  int2_u2_j1b2 = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, tmp)                   & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2,           &
+ !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,       & 
+ !$OMP          List_all_comb_b3_cent, int2_u2_j1b2)
+ !$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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp = 0.d0
+        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)
+
+          ! ---
+
+          int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
+          tmp += coef_fit * int_fit
+!          if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle
+
+          ! ---
+
+          do i_1s = 2, List_all_comb_b3_size
+         
+            coef        = List_all_comb_b3_coef  (i_1s)
+            beta        = List_all_comb_b3_expo  (i_1s)
+            B_center(1) = List_all_comb_b3_cent(1,i_1s)
+            B_center(2) = List_all_comb_b3_cent(2,i_1s)
+            B_center(3) = List_all_comb_b3_cent(3,i_1s)
+
+            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(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
+        int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_u2_j1b2', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_points_final_grid, 3)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint, i_1s, i_fit
+  double precision :: r(3), int_fit(3), expo_fit, coef_fit
+  double precision :: coef, beta, B_center(3), dist
+  double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp
+  double precision :: tmp_x, tmp_y, tmp_z
+  double precision :: wall0, wall1
+
+  print*, ' providing int2_u_grad1u_x_j1b2 ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  int2_u_grad1u_x_j1b2 = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, alpha_1s, dist,        &
+ !$OMP          alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp,     & 
+ !$OMP          tmp_x, tmp_y, tmp_z)                                & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
+ !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,       & 
+ !$OMP          List_all_comb_b3_cent, int2_u_grad1u_x_j1b2)
+ !$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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp_x = 0.d0
+        tmp_y = 0.d0
+        tmp_z = 0.d0
+        do i_fit = 1, ng_fit_jast
+
+          expo_fit = expo_gauss_j_mu_1_erf(i_fit)
+          coef_fit = coef_gauss_j_mu_1_erf(i_fit)
+
+          ! ---
+
+          call NAI_pol_x_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r, int_fit)
+          tmp_x += coef_fit * int_fit(1)
+          tmp_y += coef_fit * int_fit(2)
+          tmp_z += coef_fit * int_fit(3)
+!          if( dabs(coef_fit)*(dabs(int_fit(1)) + dabs(int_fit(2)) + dabs(int_fit(3))) .lt. 3d-10 ) cycle
+
+          ! ---
+
+          do i_1s = 2, List_all_comb_b3_size
+          
+            coef        = List_all_comb_b3_coef  (i_1s)
+            beta        = List_all_comb_b3_expo  (i_1s)
+            B_center(1) = List_all_comb_b3_cent(1,i_1s)
+            B_center(2) = List_all_comb_b3_cent(2,i_1s)
+            B_center(3) = List_all_comb_b3_cent(3,i_1s)
+            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)) 
+
+            alpha_1s     = beta + expo_fit
+            alpha_1s_inv = 1.d0 / alpha_1s 
+
+            centr_1s(1)  = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
+            centr_1s(2)  = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
+            centr_1s(3)  = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
+
+            expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist 
+            coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
+!            if(dabs(coef_tmp) .lt. 1d-12) cycle
+            
+            call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
+
+            tmp_x += coef_tmp * int_fit(1)
+            tmp_y += coef_tmp * int_fit(2)
+            tmp_z += coef_tmp * int_fit(3)
+          enddo
+
+          ! ---
+
+        enddo
+
+        int2_u_grad1u_x_j1b2(j,i,ipoint,1) = tmp_x
+        int2_u_grad1u_x_j1b2(j,i,ipoint,2) = tmp_y
+        int2_u_grad1u_x_j1b2(j,i,ipoint,3) = tmp_z
+      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
+        int2_u_grad1u_x_j1b2(j,i,ipoint,1) = int2_u_grad1u_x_j1b2(i,j,ipoint,1)
+        int2_u_grad1u_x_j1b2(j,i,ipoint,2) = int2_u_grad1u_x_j1b2(i,j,ipoint,2)
+        int2_u_grad1u_x_j1b2(j,i,ipoint,3) = int2_u_grad1u_x_j1b2(i,j,ipoint,3)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_u_grad1u_x_j1b2 = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu]
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_1s, i_fit
+  double precision              :: r(3), int_fit, expo_fit, coef_fit, coef_tmp
+  double precision              :: coef, beta, B_center(3), dist
+  double precision              :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
+  double precision              :: wall0, wall1
+  double precision, external    :: NAI_pol_mult_erf_ao_with1s
+
+  print*, ' providing int2_u_grad1u_j1b2 ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  int2_u_grad1u_j1b2 = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist,   &
+ !$OMP          alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp)     & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
+ !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,       & 
+ !$OMP          List_all_comb_b3_cent, int2_u_grad1u_j1b2)
+ !$OMP DO
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      do j = i, ao_num
+        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_fit = 1, ng_fit_jast
+
+          expo_fit = expo_gauss_j_mu_1_erf(i_fit)
+          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)
+!          if(dabs(coef_fit)*dabs(int_fit) .lt. 1d-12) cycle
+
+          tmp += coef_fit * int_fit
+
+          ! ---
+
+          do i_1s = 2, List_all_comb_b3_size
+
+            coef        = List_all_comb_b3_coef  (i_1s)
+            beta        = List_all_comb_b3_expo  (i_1s)
+            B_center(1) = List_all_comb_b3_cent(1,i_1s)
+            B_center(2) = List_all_comb_b3_cent(2,i_1s)
+            B_center(3) = List_all_comb_b3_cent(3,i_1s)
+            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))
+
+            alpha_1s     = beta + expo_fit
+            alpha_1s_inv = 1.d0 / alpha_1s 
+            centr_1s(1)  = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
+            centr_1s(2)  = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
+            centr_1s(3)  = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
+
+            expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
+            if(expo_coef_1s .gt. 80.d0) cycle
+            coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
+            if(dabs(coef_tmp) .lt. 1d-12) cycle
+
+            int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s,  1.d+9, r)
+
+            tmp += coef_tmp * int_fit
+          enddo
+
+          ! ---
+
+        enddo
+
+        int2_u_grad1u_j1b2(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
+        int2_u_grad1u_j1b2(j,i,ipoint) = int2_u_grad1u_j1b2(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_u_grad1u_j1b2', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+

src/ao_many_one_e_ints/grad2_jmu_modif_vect.irp.f

@@ -0,0 +1,453 @@
+!
+!! ---
+!
+!BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
+!
+!  BEGIN_DOC
+!  !
+!  ! -\frac{1}{4} int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
+!  !
+!  END_DOC
+!
+!  implicit none
+!  integer                       :: i, j, ipoint, i_1s, i_fit
+!  integer                       :: i_mask_grid
+!  double precision              :: r(3), expo_fit, coef_fit
+!  double precision              :: coef, beta, B_center(3)
+!  double precision              :: wall0, wall1
+!
+!  integer,          allocatable :: n_mask_grid(:)
+!  double precision, allocatable :: r_mask_grid(:,:)
+!  double precision, allocatable :: int_fit_v(:)
+!
+!  print*, ' providing int2_grad1u2_grad2u2_j1b2'
+!
+!  provide mu_erf final_grid_points_transp j1b_pen
+!  call wall_time(wall0)
+!
+!  int2_grad1u2_grad2u2_j1b2(:,:,:) = 0.d0
+!
+! !$OMP PARALLEL DEFAULT (NONE)                                     &
+! !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,&
+! !$OMP          coef_fit, expo_fit, int_fit_v, n_mask_grid,        &
+! !$OMP          i_mask_grid, r_mask_grid)                          &
+! !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size,&
+! !$OMP          final_grid_points_transp, n_max_fit_slat,          &
+! !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,        &
+! !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,      &
+! !$OMP          List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2,  &
+! !$OMP          ao_overlap_abs)
+!
+!  allocate(int_fit_v(n_points_final_grid))
+!  allocate(n_mask_grid(n_points_final_grid))
+!  allocate(r_mask_grid(n_points_final_grid,3))
+!
+! !$OMP DO SCHEDULE(dynamic)
+!  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_fit = 1, n_max_fit_slat
+!
+!        expo_fit = expo_gauss_1_erf_x_2(i_fit)
+!        coef_fit = coef_gauss_1_erf_x_2(i_fit) * (-0.25d0)
+!
+!        ! ---
+!
+!        call overlap_gauss_r12_ao_v(final_grid_points_transp, n_points_final_grid, expo_fit, i, j, int_fit_v, n_points_final_grid, n_points_final_grid)
+!
+!        i_mask_grid = 0    ! dim
+!        n_mask_grid = 0    ! ind
+!        r_mask_grid = 0.d0 ! val
+!        do ipoint = 1, n_points_final_grid
+!
+!          int2_grad1u2_grad2u2_j1b2(j,i,ipoint) += coef_fit * int_fit_v(ipoint)
+!
+!          if(dabs(int_fit_v(ipoint)) .gt. 1d-10) then
+!            i_mask_grid += 1
+!            n_mask_grid(i_mask_grid  ) = ipoint
+!            r_mask_grid(i_mask_grid,1) = final_grid_points_transp(ipoint,1)
+!            r_mask_grid(i_mask_grid,2) = final_grid_points_transp(ipoint,2)
+!            r_mask_grid(i_mask_grid,3) = final_grid_points_transp(ipoint,3)
+!          endif
+!
+!        enddo
+!
+!        if(i_mask_grid .eq. 0) cycle
+!
+!        ! ---
+!
+!        do i_1s = 2, List_all_comb_b3_size
+!
+!          coef        = List_all_comb_b3_coef  (i_1s) * coef_fit
+!          beta        = List_all_comb_b3_expo  (i_1s)
+!          B_center(1) = List_all_comb_b3_cent(1,i_1s)
+!          B_center(2) = List_all_comb_b3_cent(2,i_1s)
+!          B_center(3) = List_all_comb_b3_cent(3,i_1s)
+!
+!          call overlap_gauss_r12_ao_with1s_v(B_center, beta, r_mask_grid, n_points_final_grid, expo_fit, i, j, int_fit_v, n_points_final_grid, i_mask_grid)
+!
+!          do ipoint = 1, i_mask_grid
+!            int2_grad1u2_grad2u2_j1b2(j,i,n_mask_grid(ipoint)) += coef * int_fit_v(ipoint)
+!          enddo
+!
+!        enddo
+!
+!        ! ---
+!
+!      enddo
+!    enddo
+!  enddo
+! !$OMP END DO
+!
+!  deallocate(n_mask_grid)
+!  deallocate(r_mask_grid)
+!  deallocate(int_fit_v)
+!
+! !$OMP END PARALLEL
+!
+!  do ipoint = 1, n_points_final_grid
+!    do i = 2, ao_num
+!      do j = 1, i-1
+!        int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
+!      enddo
+!    enddo
+!  enddo
+!
+!  call wall_time(wall1)
+!  print*, ' wall time for int2_grad1u2_grad2u2_j1b2', wall1 - wall0
+!
+!END_PROVIDER
+!
+!! ---
+!
+!BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
+!
+!  BEGIN_DOC
+!  !
+!  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2
+!  !
+!  END_DOC
+!
+!  implicit none
+!  integer                       :: i, j, ipoint, i_1s, i_fit
+!  integer                       :: i_mask_grid
+!  double precision              :: r(3), expo_fit, coef_fit
+!  double precision              :: coef, beta, B_center(3), tmp
+!  double precision              :: wall0, wall1
+!
+!  integer,          allocatable :: n_mask_grid(:)
+!  double precision, allocatable :: r_mask_grid(:,:)
+!  double precision, allocatable :: int_fit_v(:)
+!
+!  print*, ' providing int2_u2_j1b2'
+!
+!  provide mu_erf final_grid_points_transp j1b_pen
+!  call wall_time(wall0)
+!
+!  int2_u2_j1b2(:,:,:) = 0.d0
+!
+! !$OMP PARALLEL DEFAULT (NONE)                                      &
+! !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+! !$OMP          coef_fit, expo_fit, int_fit_v,                      &
+! !$OMP          i_mask_grid, n_mask_grid, r_mask_grid )             &
+! !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size, &
+! !$OMP          final_grid_points_transp, n_max_fit_slat,           &
+! !$OMP          expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2,           &
+! !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,       &
+! !$OMP          List_all_comb_b3_cent, int2_u2_j1b2)
+!
+!  allocate(n_mask_grid(n_points_final_grid))
+!  allocate(r_mask_grid(n_points_final_grid,3))
+!  allocate(int_fit_v(n_points_final_grid))
+!
+! !$OMP DO SCHEDULE(dynamic)
+!  do i = 1, ao_num
+!    do j = i, ao_num
+!
+!      do i_fit = 1, n_max_fit_slat
+!
+!        expo_fit = expo_gauss_j_mu_x_2(i_fit)
+!        coef_fit = coef_gauss_j_mu_x_2(i_fit)
+!
+!        ! ---
+!
+!        call overlap_gauss_r12_ao_v(final_grid_points_transp, n_points_final_grid, expo_fit, i, j, int_fit_v, n_points_final_grid, n_points_final_grid)
+!
+!        i_mask_grid = 0    ! dim
+!        n_mask_grid = 0    ! ind
+!        r_mask_grid = 0.d0 ! val
+!
+!        do ipoint = 1, n_points_final_grid
+!          int2_u2_j1b2(j,i,ipoint) += coef_fit * int_fit_v(ipoint)
+!
+!          if(dabs(int_fit_v(ipoint)) .gt. 1d-10) then
+!            i_mask_grid += 1
+!            n_mask_grid(i_mask_grid  ) = ipoint
+!            r_mask_grid(i_mask_grid,1) = final_grid_points_transp(ipoint,1)
+!            r_mask_grid(i_mask_grid,2) = final_grid_points_transp(ipoint,2)
+!            r_mask_grid(i_mask_grid,3) = final_grid_points_transp(ipoint,3)
+!          endif
+!        enddo
+!
+!        if(i_mask_grid .eq. 0) cycle
+!
+!        ! ---
+!
+!        do i_1s = 2, List_all_comb_b3_size
+!
+!          coef        = List_all_comb_b3_coef  (i_1s) * coef_fit
+!          beta        = List_all_comb_b3_expo  (i_1s)
+!          B_center(1) = List_all_comb_b3_cent(1,i_1s)
+!          B_center(2) = List_all_comb_b3_cent(2,i_1s)
+!          B_center(3) = List_all_comb_b3_cent(3,i_1s)
+!
+!          call overlap_gauss_r12_ao_with1s_v(B_center, beta, r_mask_grid, n_points_final_grid, expo_fit, i, j, int_fit_v, n_points_final_grid, i_mask_grid)
+!
+!          do ipoint = 1, i_mask_grid
+!            int2_u2_j1b2(j,i,n_mask_grid(ipoint)) += coef * int_fit_v(ipoint)
+!          enddo
+!
+!        enddo
+!
+!        ! ---
+!
+!      enddo
+!    enddo
+!  enddo
+! !$OMP END DO
+!
+!  deallocate(n_mask_grid)
+!  deallocate(r_mask_grid)
+!  deallocate(int_fit_v)
+!
+! !$OMP END PARALLEL
+!
+!  do ipoint = 1, n_points_final_grid
+!    do i = 2, ao_num
+!      do j = 1, i-1
+!        int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint)
+!      enddo
+!    enddo
+!  enddo
+!
+!  call wall_time(wall1)
+!  print*, ' wall time for int2_u2_j1b2', wall1 - wall0
+!
+!END_PROVIDER
+!
+!! ---
+!
+!BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_points_final_grid, 3)]
+!
+!  BEGIN_DOC
+!  !
+!  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
+!  !
+!  END_DOC
+!
+!  implicit none
+!
+!  integer                       :: i, j, ipoint, i_1s, i_fit
+!  integer                       :: i_mask_grid1, i_mask_grid2, i_mask_grid3, i_mask_grid(3)
+!  double precision              :: x, y, z, expo_fit, coef_fit
+!  double precision              :: coef, beta, B_center(3)
+!  double precision              :: alpha_1s, alpha_1s_inv, expo_coef_1s
+!  double precision              :: wall0, wall1
+!
+!  integer,          allocatable :: n_mask_grid(:,:)
+!  double precision, allocatable :: r_mask_grid(:,:,:)
+!  double precision, allocatable :: int_fit_v(:,:), dist(:,:), centr_1s(:,:,:)
+!
+!  print*, ' providing int2_u_grad1u_x_j1b2'
+!
+!  provide mu_erf final_grid_points_transp j1b_pen
+!  call wall_time(wall0)
+!
+!  int2_u_grad1u_x_j1b2(:,:,:,:) = 0.d0
+!
+! !$OMP PARALLEL DEFAULT (NONE)                                          &
+! !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, x, y, z, coef, beta,         &
+! !$OMP          coef_fit, expo_fit, int_fit_v, alpha_1s, dist, B_center,&
+! !$OMP          alpha_1s_inv, centr_1s, expo_coef_1s,                   &
+! !$OMP          i_mask_grid1, i_mask_grid2, i_mask_grid3, i_mask_grid,  &
+! !$OMP          n_mask_grid, r_mask_grid)                               &
+! !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b3_size,     &
+! !$OMP          final_grid_points_transp, n_max_fit_slat,               &
+! !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,           &
+! !$OMP          List_all_comb_b3_coef, List_all_comb_b3_expo,           &
+! !$OMP          List_all_comb_b3_cent, int2_u_grad1u_x_j1b2)
+!
+!  allocate(dist(n_points_final_grid,3))
+!  allocate(centr_1s(n_points_final_grid,3,3))
+!  allocate(n_mask_grid(n_points_final_grid,3))
+!  allocate(r_mask_grid(n_points_final_grid,3,3))
+!  allocate(int_fit_v(n_points_final_grid,3))
+!
+! !$OMP DO SCHEDULE(dynamic)
+!  do i = 1, ao_num
+!    do j = i, ao_num
+!      do i_fit = 1, n_max_fit_slat
+!
+!        expo_fit = expo_gauss_j_mu_1_erf(i_fit)
+!        coef_fit = coef_gauss_j_mu_1_erf(i_fit)
+!
+!        ! ---
+!
+!        call NAI_pol_x_mult_erf_ao_with1s_v0(i, j, expo_fit, final_grid_points_transp, n_points_final_grid, 1.d+9, final_grid_points_transp, n_points_final_grid, int_fit_v, n_points_final_grid, n_points_final_grid)
+!
+!        i_mask_grid1 = 0    ! dim
+!        i_mask_grid2 = 0    ! dim
+!        i_mask_grid3 = 0    ! dim
+!        n_mask_grid  = 0    ! ind
+!        r_mask_grid  = 0.d0 ! val
+!        do ipoint = 1, n_points_final_grid
+!
+!          ! ---
+!
+!          int2_u_grad1u_x_j1b2(j,i,ipoint,1) += coef_fit * int_fit_v(ipoint,1)
+!
+!          if(dabs(int_fit_v(ipoint,1)) .gt. 1d-10) then
+!            i_mask_grid1 += 1
+!            n_mask_grid(i_mask_grid1,  1) = ipoint
+!            r_mask_grid(i_mask_grid1,1,1) = final_grid_points_transp(ipoint,1)
+!            r_mask_grid(i_mask_grid1,2,1) = final_grid_points_transp(ipoint,2)
+!            r_mask_grid(i_mask_grid1,3,1) = final_grid_points_transp(ipoint,3)
+!          endif
+!
+!          ! ---
+!
+!          int2_u_grad1u_x_j1b2(j,i,ipoint,2) += coef_fit * int_fit_v(ipoint,2)
+!
+!          if(dabs(int_fit_v(ipoint,2)) .gt. 1d-10) then
+!            i_mask_grid2 += 1
+!            n_mask_grid(i_mask_grid2,  2) = ipoint
+!            r_mask_grid(i_mask_grid2,1,2) = final_grid_points_transp(ipoint,1)
+!            r_mask_grid(i_mask_grid2,2,2) = final_grid_points_transp(ipoint,2)
+!            r_mask_grid(i_mask_grid2,3,2) = final_grid_points_transp(ipoint,3)
+!          endif
+!
+!          ! ---
+!
+!          int2_u_grad1u_x_j1b2(j,i,ipoint,3) += coef_fit * int_fit_v(ipoint,3)
+!
+!          if(dabs(int_fit_v(ipoint,3)) .gt. 1d-10) then
+!            i_mask_grid3 += 1
+!            n_mask_grid(i_mask_grid3,  3) = ipoint
+!            r_mask_grid(i_mask_grid3,1,3) = final_grid_points_transp(ipoint,1)
+!            r_mask_grid(i_mask_grid3,2,3) = final_grid_points_transp(ipoint,2)
+!            r_mask_grid(i_mask_grid3,3,3) = final_grid_points_transp(ipoint,3)
+!          endif
+!
+!          ! ---
+!
+!        enddo
+!
+!        if((i_mask_grid1+i_mask_grid2+i_mask_grid3) .eq. 0) cycle
+!
+!        i_mask_grid(1) = i_mask_grid1
+!        i_mask_grid(2) = i_mask_grid2
+!        i_mask_grid(3) = i_mask_grid3
+!
+!        ! ---
+!
+!        do i_1s = 2, List_all_comb_b3_size
+!
+!          coef        = List_all_comb_b3_coef  (i_1s) * coef_fit
+!          beta        = List_all_comb_b3_expo  (i_1s)
+!          B_center(1) = List_all_comb_b3_cent(1,i_1s)
+!          B_center(2) = List_all_comb_b3_cent(2,i_1s)
+!          B_center(3) = List_all_comb_b3_cent(3,i_1s)
+!
+!          alpha_1s     = beta + expo_fit
+!          alpha_1s_inv = 1.d0 / alpha_1s
+!          expo_coef_1s = beta * expo_fit * alpha_1s_inv 
+!
+!          do ipoint = 1, i_mask_grid1
+!
+!            x = r_mask_grid(ipoint,1,1)
+!            y = r_mask_grid(ipoint,2,1)
+!            z = r_mask_grid(ipoint,3,1)
+!
+!            centr_1s(ipoint,1,1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * x)
+!            centr_1s(ipoint,2,1) = alpha_1s_inv * (beta * B_center(2) + expo_fit * y)
+!            centr_1s(ipoint,3,1) = alpha_1s_inv * (beta * B_center(3) + expo_fit * z)
+!
+!            dist(ipoint,1) = (B_center(1) - x) * (B_center(1) - x) + (B_center(2) - y) * (B_center(2) - y) + (B_center(3) - z) * (B_center(3) - z)
+!          enddo
+!
+!          do ipoint = 1, i_mask_grid2
+!
+!            x = r_mask_grid(ipoint,1,2)
+!            y = r_mask_grid(ipoint,2,2)
+!            z = r_mask_grid(ipoint,3,2)
+!
+!            centr_1s(ipoint,1,2) = alpha_1s_inv * (beta * B_center(1) + expo_fit * x)
+!            centr_1s(ipoint,2,2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * y)
+!            centr_1s(ipoint,3,2) = alpha_1s_inv * (beta * B_center(3) + expo_fit * z)
+!
+!            dist(ipoint,2) = (B_center(1) - x) * (B_center(1) - x) + (B_center(2) - y) * (B_center(2) - y) + (B_center(3) - z) * (B_center(3) - z)
+!          enddo
+!
+!          do ipoint = 1, i_mask_grid3
+!
+!            x = r_mask_grid(ipoint,1,3)
+!            y = r_mask_grid(ipoint,2,3)
+!            z = r_mask_grid(ipoint,3,3)
+!
+!            centr_1s(ipoint,1,3) = alpha_1s_inv * (beta * B_center(1) + expo_fit * x)
+!            centr_1s(ipoint,2,3) = alpha_1s_inv * (beta * B_center(2) + expo_fit * y)
+!            centr_1s(ipoint,3,3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * z)
+!
+!            dist(ipoint,3) = (B_center(1) - x) * (B_center(1) - x) + (B_center(2) - y) * (B_center(2) - y) + (B_center(3) - z) * (B_center(3) - z)
+!          enddo
+!
+!          call NAI_pol_x_mult_erf_ao_with1s_v(i, j, alpha_1s, centr_1s, n_points_final_grid, 1.d+9, r_mask_grid, n_points_final_grid, int_fit_v, n_points_final_grid, i_mask_grid)
+!
+!          do ipoint = 1, i_mask_grid1
+!            int2_u_grad1u_x_j1b2(j,i,n_mask_grid(ipoint,1),1) += coef * dexp(-expo_coef_1s * dist(ipoint,1)) * int_fit_v(ipoint,1)
+!          enddo
+!
+!          do ipoint = 1, i_mask_grid2
+!            int2_u_grad1u_x_j1b2(j,i,n_mask_grid(ipoint,2),2) += coef * dexp(-expo_coef_1s * dist(ipoint,2)) * int_fit_v(ipoint,2)
+!          enddo
+!
+!          do ipoint = 1, i_mask_grid3
+!            int2_u_grad1u_x_j1b2(j,i,n_mask_grid(ipoint,3),3) += coef * dexp(-expo_coef_1s * dist(ipoint,3)) * int_fit_v(ipoint,3)
+!          enddo
+!
+!        enddo
+!
+!        ! ---
+!
+!      enddo
+!    enddo
+!  enddo
+! !$OMP END DO
+!
+!  deallocate(dist)
+!  deallocate(centr_1s)
+!  deallocate(n_mask_grid)
+!  deallocate(r_mask_grid)
+!  deallocate(int_fit_v)
+!
+! !$OMP END PARALLEL
+!
+!  do ipoint = 1, n_points_final_grid
+!    do i = 2, ao_num
+!      do j = 1, i-1
+!        int2_u_grad1u_x_j1b2(j,i,ipoint,1) = int2_u_grad1u_x_j1b2(i,j,ipoint,1)
+!        int2_u_grad1u_x_j1b2(j,i,ipoint,2) = int2_u_grad1u_x_j1b2(i,j,ipoint,2)
+!        int2_u_grad1u_x_j1b2(j,i,ipoint,3) = int2_u_grad1u_x_j1b2(i,j,ipoint,3)
+!      enddo
+!    enddo
+!  enddo
+!
+!  call wall_time(wall1)
+!  print*, ' wall time for int2_u_grad1u_x_j1b2 =', wall1 - wall0
+!
+!END_PROVIDER
+!

src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f

@@ -0,0 +1,369 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R|
+  !
+  END_DOC
+
+  implicit none
+  integer                    :: i, j, ipoint, i_1s
+  double precision           :: r(3), int_mu, int_coulomb
+  double precision           :: coef, beta, B_center(3)
+  double precision           :: tmp,int_j1b
+  double precision           :: wall0, wall1
+  double precision, external :: NAI_pol_mult_erf_ao_with1s
+  double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
+
+  print*, ' providing v_ij_erf_rk_cst_mu_j1b_test ...'
+
+  dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
+  provide mu_erf final_grid_points j1b_pen
+  call wall_time(wall0)
+
+  v_ij_erf_rk_cst_mu_j1b_test = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                         &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_j1b)& 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
+ !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b,  &
+ !$OMP          v_ij_erf_rk_cst_mu_j1b_test, mu_erf,                                   &
+ !$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)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, 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)
+
+          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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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)
+          ! TODO :: cycle on the 1 - erf(mur12)
+          int_mu      = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
+          int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r)
+
+          tmp += coef * (int_mu - int_coulomb)
+        enddo
+
+        v_ij_erf_rk_cst_mu_j1b_test(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_erf_rk_cst_mu_j1b_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print*, ' wall time for v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid, 3)]
+
+  BEGIN_DOC
+  ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint, i_1s
+  double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
+  double precision :: tmp_x, tmp_y, tmp_z
+  double precision :: wall0, wall1
+  double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b,factor_ij_1s,beta_ij,center_ij_1s
+
+  print*, ' providing x_v_ij_erf_rk_cst_mu_j1b_test ...'
+
+  dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
+
+  provide expo_erfc_mu_gauss ao_prod_sigma ao_prod_center
+  call wall_time(wall0)
+
+  x_v_ij_erf_rk_cst_mu_j1b_test = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                        &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb,      & 
+ !$OMP          int_j1b, tmp_x, tmp_y, tmp_z,factor_ij_1s,beta_ij,center_ij_1s)       & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
+ !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,  &
+ !$OMP          x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b,         &
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma)
+! !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
+ !$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)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
+
+        tmp_x = 0.d0
+        tmp_y = 0.d0
+        tmp_z = 0.d0
+        do i_1s = 1, 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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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)
+
+!          if(ao_prod_center(1,j,i).ne.10000.d0)then
+!           ! approximate 1 - erf(mu r12) by a gaussian * 10
+!           !DIR$ FORCEINLINE
+!           call gaussian_product(expo_erfc_mu_gauss,r,     &
+!                ao_prod_sigma(j,i),ao_prod_center(1,j,i),  & 
+!                factor_ij_1s,beta_ij,center_ij_1s)
+!           if(dabs(coef * factor_ij_1s*int_j1b*10.d0 * dsqpi_3_2 * beta_ij**(-1.5d0)).lt.1.d-10)cycle 
+!          endif
+          call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints        )
+          call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r, ints_coulomb)
+
+          tmp_x += coef * (ints(1) - ints_coulomb(1))
+          tmp_y += coef * (ints(2) - ints_coulomb(2))
+          tmp_z += coef * (ints(3) - ints_coulomb(3))
+        enddo
+
+        x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = tmp_x
+        x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = tmp_y
+        x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = tmp_z
+      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
+        x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1)
+        x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2)
+        x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+! TODO analytically
+BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (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
+
+  implicit none
+  integer                    :: i, j, ipoint, i_1s, i_fit
+  double precision           :: r(3), int_fit, expo_fit, coef_fit
+  double precision           :: coef, beta, B_center(3)
+  double precision           :: tmp
+  double precision           :: wall0, wall1
+
+  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 ...'
+
+  dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
+
+  provide mu_erf final_grid_points j1b_pen
+  call wall_time(wall0)
+
+  v_ij_u_cst_mu_j1b_test = 0.d0
+
+ !$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          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          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)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, 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)
+
+          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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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
+
+        v_ij_u_cst_mu_j1b_test(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_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_test(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print*, ' wall time for v_ij_u_cst_mu_j1b_test', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) with u(mu,r12) \approx 1/2 mu e^{-2.5 * mu (r12)^2}
+  !
+  END_DOC
+
+  implicit none
+  integer                    :: i, j, ipoint, i_1s
+  double precision           :: r(3), int_fit, expo_fit, coef_fit
+  double precision           :: coef, beta, B_center(3)
+  double precision           :: tmp
+  double precision           :: wall0, wall1
+
+  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
+  call wall_time(wall0)
+
+  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          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          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)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, 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)
+
+          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)
+          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          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
+
+        v_ij_u_cst_mu_j1b_ng_1_test(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_ng_1_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_ng_1_test(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print*, ' wall time for v_ij_u_cst_mu_j1b_ng_1_test', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+

src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f

@@ -0,0 +1,300 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R|
+  !
+  END_DOC
+
+  implicit none
+  integer                    :: i, j, ipoint, i_1s
+  double precision           :: r(3), int_mu, int_coulomb
+  double precision           :: coef, beta, B_center(3)
+  double precision           :: tmp
+  double precision           :: wall0, wall1
+  double precision, external :: NAI_pol_mult_erf_ao_with1s
+
+  print *, ' providing v_ij_erf_rk_cst_mu_j1b ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  v_ij_erf_rk_cst_mu_j1b = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                         &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, &
+ !$OMP          List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent,   &
+ !$OMP          v_ij_erf_rk_cst_mu_j1b, mu_erf)
+ !$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)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp = 0.d0
+
+        ! ---
+
+        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)
+
+        int_mu      = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
+        int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r)
+!        if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle
+
+        tmp += coef * (int_mu - int_coulomb)
+
+        ! ---
+
+        do i_1s = 2, List_all_comb_b2_size
+
+          coef        = List_all_comb_b2_coef  (i_1s)
+          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)
+
+          int_mu      = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
+          int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r)
+
+          tmp += coef * (int_mu - int_coulomb)
+        enddo
+
+        ! ---
+
+        v_ij_erf_rk_cst_mu_j1b(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_erf_rk_cst_mu_j1b(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print*, ' wall time for v_ij_erf_rk_cst_mu_j1b', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid, 3)]
+
+  BEGIN_DOC
+  ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint, i_1s
+  double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
+  double precision :: tmp_x, tmp_y, tmp_z
+  double precision :: wall0, wall1
+
+  print*, ' providing x_v_ij_erf_rk_cst_mu_j1b ...'
+  call wall_time(wall0)
+
+  x_v_ij_erf_rk_cst_mu_j1b = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                        &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb,      & 
+ !$OMP          tmp_x, tmp_y, tmp_z)                                                  & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,&
+ !$OMP          List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent,  &
+ !$OMP          x_v_ij_erf_rk_cst_mu_j1b, mu_erf)
+ !$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)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp_x = 0.d0
+        tmp_y = 0.d0
+        tmp_z = 0.d0
+
+        ! ---
+
+        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_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints        )
+        call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r, ints_coulomb)
+
+!        if( dabs(coef)*(dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle
+
+        tmp_x += coef * (ints(1) - ints_coulomb(1))
+        tmp_y += coef * (ints(2) - ints_coulomb(2))
+        tmp_z += coef * (ints(3) - ints_coulomb(3))
+
+        ! ---
+
+        do i_1s = 2, List_all_comb_b2_size
+
+          coef        = List_all_comb_b2_coef  (i_1s)
+          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_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints        )
+          call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r, ints_coulomb)
+
+          tmp_x += coef * (ints(1) - ints_coulomb(1))
+          tmp_y += coef * (ints(2) - ints_coulomb(2))
+          tmp_z += coef * (ints(3) - ints_coulomb(3))
+        enddo
+
+        ! ---
+
+        x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = tmp_x
+        x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = tmp_y
+        x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = tmp_z
+      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
+        x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1)
+        x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2)
+        x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b =', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+! TODO analytically
+BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (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
+
+  implicit none
+  integer                    :: i, j, ipoint, i_1s, i_fit
+  double precision           :: r(3), int_fit, expo_fit, coef_fit
+  double precision           :: coef, beta, B_center(3)
+  double precision           :: tmp
+  double precision           :: wall0, wall1
+
+  double precision, external :: overlap_gauss_r12_ao_with1s
+
+  print*, ' providing v_ij_u_cst_mu_j1b ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  v_ij_u_cst_mu_j1b = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
+ !$OMP          coef_fit, expo_fit, int_fit, tmp)                   & 
+ !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b2_size, & 
+ !$OMP          final_grid_points, ng_fit_jast,                     &
+ !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,               &
+ !$OMP          List_all_comb_b2_coef, List_all_comb_b2_expo,       & 
+ !$OMP          List_all_comb_b2_cent, v_ij_u_cst_mu_j1b)
+ !$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)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp = 0.d0
+        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)
+
+          ! ---
+
+          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)
+
+          int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
+!          if(dabs(int_fit*coef) .lt. 1d-12) cycle
+
+          tmp += coef * coef_fit * int_fit
+
+          ! ---
+
+          do i_1s = 2, List_all_comb_b2_size
+
+            coef        = List_all_comb_b2_coef  (i_1s)
+            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)
+
+            int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
+
+            tmp += coef * coef_fit * int_fit
+          enddo
+
+          ! ---
+
+        enddo
+
+        v_ij_u_cst_mu_j1b(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(j,i,ipoint) = v_ij_u_cst_mu_j1b(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print*, ' wall time for v_ij_u_cst_mu_j1b', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+

src/ao_many_one_e_ints/grad_related_ints.irp.f

@@ -0,0 +1,437 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1) / |r - R|
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: r(3)
+  double precision :: int_mu, int_coulomb
+  double precision :: wall0, wall1
+
+  double precision :: NAI_pol_mult_erf_ao
+
+  print*, ' providing v_ij_erf_rk_cst_mu ...'
+
+  provide mu_erf final_grid_points 
+  call wall_time(wall0)
+
+  v_ij_erf_rk_cst_mu = 0.d0
+
+ !$OMP PARALLEL                                       &
+ !$OMP DEFAULT (NONE)                                 &
+ !$OMP PRIVATE (i, j, ipoint, r, int_mu, int_coulomb) & 
+ !$OMP SHARED  (ao_num, n_points_final_grid, v_ij_erf_rk_cst_mu, final_grid_points, mu_erf)
+ !$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
+
+        int_mu      = NAI_pol_mult_erf_ao(i, j, mu_erf, r)
+        int_coulomb = NAI_pol_mult_erf_ao(i, j,  1.d+9, r)
+
+        v_ij_erf_rk_cst_mu(j,i,ipoint) = int_mu - int_coulomb
+      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_erf_rk_cst_mu(j,i,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print*, ' wall time for v_ij_erf_rk_cst_mu = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_transp, (n_points_final_grid, ao_num, ao_num)]
+
+  BEGIN_DOC
+  ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: r(3)
+  double precision :: int_mu, int_coulomb
+  double precision :: wall0, wall1
+  double precision :: NAI_pol_mult_erf_ao
+
+  print *, ' providing v_ij_erf_rk_cst_mu_transp ...'
+
+  provide mu_erf final_grid_points 
+  call wall_time(wall0)
+
+ !$OMP PARALLEL                                  &
+ !$OMP DEFAULT (NONE)                            &
+ !$OMP PRIVATE (i,j,ipoint,r,int_mu,int_coulomb) & 
+ !$OMP SHARED (ao_num,n_points_final_grid,v_ij_erf_rk_cst_mu_transp,final_grid_points,mu_erf)
+ !$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
+        int_mu      = NAI_pol_mult_erf_ao(i, j, mu_erf, r)
+        int_coulomb = NAI_pol_mult_erf_ao(i, j,  1.d+9, r)
+
+        v_ij_erf_rk_cst_mu_transp(ipoint,j,i) = int_mu - int_coulomb
+      enddo
+    enddo
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  do i = 2, ao_num
+    do j = 1, i-1
+      do ipoint = 1, n_points_final_grid
+        v_ij_erf_rk_cst_mu_transp(ipoint,j,i) = v_ij_erf_rk_cst_mu_transp(ipoint,i,j)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print *, ' wall time for v_ij_erf_rk_cst_mu_transp = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: r(3), ints(3), ints_coulomb(3)
+  double precision :: wall0, wall1
+
+  print*, ' providing x_v_ij_erf_rk_cst_mu_tmp ...'
+
+  call wall_time(wall0)
+
+ !$OMP PARALLEL                                 &
+ !$OMP DEFAULT (NONE)                           &
+ !$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) & 
+ !$OMP SHARED (ao_num,n_points_final_grid,x_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf)
+ !$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
+
+        call NAI_pol_x_mult_erf_ao(i, j, mu_erf, r, ints        )
+        call NAI_pol_x_mult_erf_ao(i, j, 1.d+9 , r, ints_coulomb)
+
+        x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = ints(1) - ints_coulomb(1)
+        x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = ints(2) - ints_coulomb(2)
+        x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = ints(3) - ints_coulomb(3)
+      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
+        x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint)
+        x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint)
+        x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print *, ' wall time for x_v_ij_erf_rk_cst_mu_tmp = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points_final_grid, 3)]
+
+  BEGIN_DOC
+  ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  print *, ' providing x_v_ij_erf_rk_cst_mu ...'
+
+  call wall_time(wall0)
+
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      do j = 1, ao_num
+        x_v_ij_erf_rk_cst_mu(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint)
+        x_v_ij_erf_rk_cst_mu(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint)
+        x_v_ij_erf_rk_cst_mu(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print *, ' wall time for x_v_ij_erf_rk_cst_mu = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_transp, (ao_num, ao_num,3,n_points_final_grid)]
+
+  BEGIN_DOC
+  ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  print *, ' providing x_v_ij_erf_rk_cst_mu_transp ...'
+
+  call wall_time(wall0)
+
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      do j = 1, ao_num
+        x_v_ij_erf_rk_cst_mu_transp(j,i,1,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint)
+        x_v_ij_erf_rk_cst_mu_transp(j,i,2,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint)
+        x_v_ij_erf_rk_cst_mu_transp(j,i,3,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print *, ' wall time for x_v_ij_erf_rk_cst_mu_transp = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_transp_bis, (n_points_final_grid, ao_num, ao_num, 3)]
+
+  BEGIN_DOC
+  ! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  print *, ' providing x_v_ij_erf_rk_cst_mu_transp_bis ...'
+
+  call wall_time(wall0)
+
+  do i = 1, ao_num
+    do j = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+        x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,1) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint)
+        x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,2) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint)
+        x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,3) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print *, ' wall time for x_v_ij_erf_rk_cst_mu_transp_bis = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu_tmp, (3, n_points_final_grid, ao_num, ao_num)]
+
+  BEGIN_DOC
+  ! d_dx_v_ij_erf_rk_cst_mu_tmp(m,R,j,i) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) 
+  !
+  ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
+  END_DOC
+
+ implicit none
+ integer          :: i, j, ipoint
+ double precision :: r(3), ints(3), ints_coulomb(3)
+ double precision :: wall0, wall1
+
+  print *, ' providing d_dx_v_ij_erf_rk_cst_mu_tmp ...'
+
+  call wall_time(wall0)
+
+ !$OMP PARALLEL                                 &
+ !$OMP DEFAULT (NONE)                           &
+ !$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) & 
+ !$OMP SHARED (ao_num,n_points_final_grid,d_dx_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf)
+ !$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 = 1, ao_num
+        call phi_j_erf_mu_r_dxyz_phi(j, i, mu_erf, r, ints)
+        call phi_j_erf_mu_r_dxyz_phi(j, i,  1.d+9, r, ints_coulomb)
+
+        d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) = ints(1) - ints_coulomb(1)
+        d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) = ints(2) - ints_coulomb(2)
+        d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) = ints(3) - ints_coulomb(3)
+      enddo
+    enddo
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  call wall_time(wall1)
+  print *, ' wall time for d_dx_v_ij_erf_rk_cst_mu_tmp = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid, ao_num, ao_num, 3)]
+
+  BEGIN_DOC
+  !
+  ! d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) 
+  !
+  ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  print *, ' providing d_dx_v_ij_erf_rk_cst_mu ...'
+
+  call wall_time(wall0)
+  do i = 1, ao_num
+    do j = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+        d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,1) = d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i)
+        d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,2) = d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i)
+        d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,3) = d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print *, ' wall time for d_dx_v_ij_erf_rk_cst_mu = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu_tmp, (3, n_points_final_grid, ao_num, ao_num)]
+
+  BEGIN_DOC
+  !
+  ! x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,j,i,R) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) 
+  !
+  ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: r(3), ints(3), ints_coulomb(3)
+  double precision :: wall0, wall1
+
+  print *, ' providing x_d_dx_v_ij_erf_rk_cst_mu_tmp ...'
+
+  call wall_time(wall0)
+
+ !$OMP PARALLEL                                 &
+ !$OMP DEFAULT (NONE)                           &
+ !$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) & 
+ !$OMP SHARED (ao_num,n_points_final_grid,x_d_dx_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf)
+ !$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 = 1, ao_num
+        call phi_j_erf_mu_r_xyz_dxyz_phi(j, i, mu_erf, r, ints)
+        call phi_j_erf_mu_r_xyz_dxyz_phi(j, i,  1.d+9, r, ints_coulomb)
+
+        x_d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) = ints(1) - ints_coulomb(1)
+        x_d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) = ints(2) - ints_coulomb(2)
+        x_d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) = ints(3) - ints_coulomb(3)
+      enddo
+    enddo
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  call wall_time(wall1)
+  print *, ' wall time for x_d_dx_v_ij_erf_rk_cst_mu_tmp = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid,ao_num, ao_num,3)]
+
+  BEGIN_DOC
+  !
+  ! x_d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r) 
+  !
+  ! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  print *, ' providing x_d_dx_v_ij_erf_rk_cst_mu ...'
+
+  call wall_time(wall0)
+
+  do i = 1, ao_num
+    do j = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+        x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,1) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i)
+        x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,2) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i)
+        x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,3) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i)
+      enddo
+    enddo
+  enddo
+ 
+  call wall_time(wall1)
+  print *, ' wall time for x_d_dx_v_ij_erf_rk_cst_mu = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+

src/ao_many_one_e_ints/list_grid.irp.f

@@ -0,0 +1,59 @@
+ BEGIN_PROVIDER [ integer, n_pts_grid_ao_prod, (ao_num, ao_num)]
+&BEGIN_PROVIDER [ integer, max_n_pts_grid_ao_prod]
+ implicit none
+ integer :: i,j,ipoint
+ double precision :: overlap, r(3),thr, overlap_abs_gauss_r12_ao,overlap_gauss_r12_ao
+ double precision :: sigma,dist,center_ij(3),fact_gauss, alpha, center(3)
+ n_pts_grid_ao_prod = 0
+ thr = 1.d-11
+ print*,' expo_good_j_mu_1gauss = ',expo_good_j_mu_1gauss
+ !$OMP PARALLEL DEFAULT (NONE)                                      &
+ !$OMP PRIVATE (ipoint, i, j, r, overlap, fact_gauss, alpha, center,dist,sigma,center_ij) &
+ !$OMP SHARED  (n_points_final_grid, ao_num, thr, ao_overlap_abs_grid,n_pts_grid_ao_prod,expo_good_j_mu_1gauss,&
+ !$OMP          final_grid_points,ao_prod_center,ao_prod_sigma,ao_nucl)
+ !$OMP DO
+ do i = 1, ao_num
+! do i = 3,3
+  do j = 1, ao_num
+! do i = 22,22
+!  do j = 9,9
+   center_ij(1:3) = ao_prod_center(1:3,j,i)
+   sigma = ao_prod_sigma(j,i)
+   sigma *= sigma
+   sigma = 0.5d0 /sigma
+!   if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
+   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)
+    dist  = (center_ij(1) - r(1))*(center_ij(1) - r(1))
+    dist += (center_ij(2) - r(2))*(center_ij(2) - r(2))
+    dist += (center_ij(3) - r(3))*(center_ij(3) - r(3))
+    dist = dsqrt(dist)
+    call gaussian_product(sigma, center_ij, expo_good_j_mu_1gauss, r, fact_gauss, alpha, center)
+!    print*,''
+!    print*,j,i,ao_overlap_abs_grid(j,i),ao_overlap_abs(j,i)
+!    print*,r
+!    print*,dist,sigma
+!    print*,fact_gauss
+    if( fact_gauss*ao_overlap_abs_grid(j,i).lt.1.d-11)cycle
+    if(ao_nucl(i) == ao_nucl(j))then
+     overlap = overlap_abs_gauss_r12_ao(r, expo_good_j_mu_1gauss, i, j)
+    else
+     overlap = overlap_gauss_r12_ao(r, expo_good_j_mu_1gauss, i, j)
+    endif
+!    print*,overlap
+    if(dabs(overlap).lt.thr)cycle
+    n_pts_grid_ao_prod(j,i) += 1
+   enddo
+  enddo
+ enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+ integer :: list(ao_num)
+ do i = 1, ao_num
+  list(i) = maxval(n_pts_grid_ao_prod(:,i))
+ enddo
+ max_n_pts_grid_ao_prod = maxval(list)
+END_PROVIDER

src/ao_many_one_e_ints/listj1b.irp.f

@@ -0,0 +1,237 @@
+
+! ---
+
+BEGIN_PROVIDER [ integer, List_all_comb_b2_size]
+
+  implicit none
+
+  List_all_comb_b2_size = 2**nucl_num
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)]
+
+  implicit none
+  integer :: i, j
+
+  if(nucl_num .gt. 32) then
+    print *, ' nucl_num = ', nucl_num, '> 32'
+    stop
+  endif
+
+  List_all_comb_b2 = 0
+
+  do i = 0, List_all_comb_b2_size-1
+    do j = 0, nucl_num-1
+      if (btest(i,j)) then
+        List_all_comb_b2(j+1,i+1) = 1
+      endif
+    enddo
+  enddo
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, List_all_comb_b2_coef, (   List_all_comb_b2_size)]
+&BEGIN_PROVIDER [ double precision, List_all_comb_b2_expo, (   List_all_comb_b2_size)]
+&BEGIN_PROVIDER [ double precision, List_all_comb_b2_cent, (3, List_all_comb_b2_size)]
+
+  implicit none
+  integer          :: i, j, k, phase
+  double precision :: tmp_alphaj, tmp_alphak
+  double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z
+
+  provide j1b_pen
+
+  List_all_comb_b2_coef = 0.d0
+  List_all_comb_b2_expo = 0.d0
+  List_all_comb_b2_cent = 0.d0
+
+  do i = 1, List_all_comb_b2_size
+
+    tmp_cent_x = 0.d0
+    tmp_cent_y = 0.d0
+    tmp_cent_z = 0.d0
+    do j = 1, nucl_num
+      tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
+      List_all_comb_b2_expo(i) += tmp_alphaj
+      tmp_cent_x               += tmp_alphaj * nucl_coord(j,1)
+      tmp_cent_y               += tmp_alphaj * nucl_coord(j,2)
+      tmp_cent_z               += tmp_alphaj * nucl_coord(j,3)
+    enddo
+
+    if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle
+
+    List_all_comb_b2_cent(1,i) = tmp_cent_x / List_all_comb_b2_expo(i) 
+    List_all_comb_b2_cent(2,i) = tmp_cent_y / List_all_comb_b2_expo(i)
+    List_all_comb_b2_cent(3,i) = tmp_cent_z / List_all_comb_b2_expo(i)
+  enddo
+
+  ! ---
+
+  do i = 1, List_all_comb_b2_size
+
+    do j = 2, nucl_num, 1
+      tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
+      do k = 1, j-1, 1
+        tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k)
+
+        List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+                                                              + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+                                                              + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
+      enddo
+    enddo
+
+    if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle
+
+    List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i)
+  enddo
+
+  ! ---
+
+  do i = 1, List_all_comb_b2_size
+
+    phase = 0
+    do j = 1, nucl_num
+      phase += List_all_comb_b2(j,i)
+    enddo
+
+    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
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ integer, List_all_comb_b3_size]
+
+  implicit none
+
+  List_all_comb_b3_size = 3**nucl_num
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)]
+
+  implicit none
+  integer              :: i, j, ii, jj
+  integer, allocatable :: M(:,:), p(:)
+
+  if(nucl_num .gt. 32) then
+    print *, ' nucl_num = ', nucl_num, '> 32'
+    stop
+  endif
+
+  List_all_comb_b3(:,:)                     = 0
+  List_all_comb_b3(:,List_all_comb_b3_size) = 2
+
+  allocate(p(nucl_num))
+  p = 0
+
+  do i = 2, List_all_comb_b3_size-1
+    do j = 1, nucl_num
+
+      ii = 0
+      do jj = 1, j-1, 1
+        ii = ii + p(jj) * 3**(jj-1)
+      enddo
+      p(j) = modulo(i-1-ii, 3**j) / 3**(j-1)
+
+      List_all_comb_b3(j,i) = p(j)
+    enddo
+  enddo
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, List_all_comb_b3_coef, (   List_all_comb_b3_size)]
+&BEGIN_PROVIDER [ double precision, List_all_comb_b3_expo, (   List_all_comb_b3_size)]
+&BEGIN_PROVIDER [ double precision, List_all_comb_b3_cent, (3, List_all_comb_b3_size)]
+
+  implicit none
+  integer          :: i, j, k, phase
+  double precision :: tmp_alphaj, tmp_alphak, facto
+
+  provide j1b_pen
+
+  List_all_comb_b3_coef = 0.d0
+  List_all_comb_b3_expo = 0.d0
+  List_all_comb_b3_cent = 0.d0
+
+  do i = 1, List_all_comb_b3_size
+
+    do j = 1, nucl_num
+      tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j)
+      List_all_comb_b3_expo(i)   += tmp_alphaj
+      List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1)
+      List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2)
+      List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3)
+
+    enddo
+
+    if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle
+    ASSERT(List_all_comb_b3_expo(i) .gt. 0d0)
+
+    List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i) 
+    List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i)
+    List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i)
+  enddo
+
+  ! ---
+
+  do i = 1, List_all_comb_b3_size
+
+    do j = 2, nucl_num, 1
+      tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j)
+      do k = 1, j-1, 1
+        tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k)
+
+        List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+                                                              + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+                                                              + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
+      enddo
+    enddo
+
+    if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle
+
+    List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i)
+  enddo
+
+  ! ---
+
+  do i = 1, List_all_comb_b3_size
+
+    facto = 1.d0
+    phase = 0
+    do j = 1, nucl_num
+      tmp_alphaj = dble(List_all_comb_b3(j,i)) 
+
+      facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj))
+      phase += List_all_comb_b3(j,i)
+    enddo
+
+    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
+
+END_PROVIDER
+
+! ---
+

src/ao_many_one_e_ints/listj1b_sorted.irp.f

@@ -0,0 +1,191 @@
+
+ BEGIN_PROVIDER [ integer, List_comb_thr_b2_size, (ao_num, ao_num)]
+&BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size]
+ implicit none
+ integer :: i_1s,i,j,ipoint
+ double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: r(3),weight,dist
+ thr = 1.d-15
+ List_comb_thr_b2_size = 0
+ do i = 1, ao_num
+  do j = i, ao_num
+   do i_1s = 1, List_all_comb_b2_size
+     coef        = List_all_comb_b2_coef  (i_1s)
+     if(dabs(coef).lt.1.d-15)cycle
+     beta        = List_all_comb_b2_expo  (i_1s)
+     beta = max(beta,1.d-12)
+     center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
+     int_j1b = 0.d0
+     do ipoint = 1, n_points_extra_final_grid
+      r(1:3) = final_grid_points_extra(1:3,ipoint)
+      weight = final_weight_at_r_vector_extra(ipoint)
+      dist  = ( center(1) - r(1) )*( center(1) - r(1) )
+      dist += ( center(2) - r(2) )*( center(2) - r(2) )
+      dist += ( center(3) - r(3) )*( center(3) - r(3) )
+      int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
+     enddo
+     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+      List_comb_thr_b2_size(j,i) += 1
+     endif
+   enddo
+  enddo 
+ enddo
+ do i = 1, ao_num
+  do j = 1, i-1
+    List_comb_thr_b2_size(j,i) = List_comb_thr_b2_size(i,j)
+  enddo
+ enddo
+ integer :: list(ao_num)
+ do i = 1, ao_num
+  list(i) = maxval(List_comb_thr_b2_size(:,i))
+ enddo
+ max_List_comb_thr_b2_size = maxval(list) 
+ 
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, (   max_List_comb_thr_b2_size,ao_num, ao_num )]
+&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, (   max_List_comb_thr_b2_size,ao_num, ao_num )]
+&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3, max_List_comb_thr_b2_size,ao_num, ao_num )]
+&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)]
+ implicit none
+ integer :: i_1s,i,j,ipoint,icount
+ double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: r(3),weight,dist
+ thr = 1.d-15
+ ao_abs_comb_b2_j1b = 10000000.d0
+ do i = 1, ao_num
+  do j = i, ao_num
+   icount = 0
+   do i_1s = 1, List_all_comb_b2_size
+     coef        = List_all_comb_b2_coef  (i_1s)
+     if(dabs(coef).lt.1.d-12)cycle
+     beta        = List_all_comb_b2_expo  (i_1s)
+     center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
+     int_j1b = 0.d0
+     do ipoint = 1, n_points_extra_final_grid
+      r(1:3) = final_grid_points_extra(1:3,ipoint)
+      weight = final_weight_at_r_vector_extra(ipoint)
+      dist  = ( center(1) - r(1) )*( center(1) - r(1) )
+      dist += ( center(2) - r(2) )*( center(2) - r(2) )
+      dist += ( center(3) - r(3) )*( center(3) - r(3) )
+      int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
+     enddo
+     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+      icount += 1
+      List_comb_thr_b2_coef(icount,j,i) = coef
+      List_comb_thr_b2_expo(icount,j,i) = beta
+      List_comb_thr_b2_cent(1:3,icount,j,i) = center(1:3)
+      ao_abs_comb_b2_j1b(icount,j,i) = int_j1b
+     endif
+   enddo
+  enddo 
+ enddo
+
+ do i = 1, ao_num
+  do j = 1, i-1
+    do icount = 1, List_comb_thr_b2_size(j,i)
+     List_comb_thr_b2_coef(icount,j,i) = List_comb_thr_b2_coef(icount,i,j)
+     List_comb_thr_b2_expo(icount,j,i) = List_comb_thr_b2_expo(icount,i,j)
+     List_comb_thr_b2_cent(1:3,icount,j,i) = List_comb_thr_b2_cent(1:3,icount,i,j)
+    enddo
+  enddo
+ enddo
+ 
+END_PROVIDER 
+
+
+ BEGIN_PROVIDER [ integer, List_comb_thr_b3_size, (ao_num, ao_num)]
+&BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size]
+ implicit none
+ integer :: i_1s,i,j,ipoint
+ double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: r(3),weight,dist
+ thr = 1.d-15
+ List_comb_thr_b3_size = 0
+ do i = 1, ao_num
+  do j = 1, ao_num
+   do i_1s = 1, List_all_comb_b3_size
+     coef        = List_all_comb_b3_coef  (i_1s)
+     beta        = List_all_comb_b3_expo  (i_1s)
+     center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
+     if(dabs(coef).lt.thr)cycle
+     int_j1b = 0.d0
+     do ipoint = 1, n_points_extra_final_grid
+      r(1:3) = final_grid_points_extra(1:3,ipoint)
+      weight = final_weight_at_r_vector_extra(ipoint)
+      dist  = ( center(1) - r(1) )*( center(1) - r(1) )
+      dist += ( center(2) - r(2) )*( center(2) - r(2) )
+      dist += ( center(3) - r(3) )*( center(3) - r(3) )
+      int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
+     enddo
+     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+      List_comb_thr_b3_size(j,i) += 1
+     endif
+   enddo
+  enddo 
+ enddo
+! do i = 1, ao_num
+!  do j = 1, i-1
+!    List_comb_thr_b3_size(j,i) = List_comb_thr_b3_size(i,j)
+!  enddo
+! enddo
+ integer :: list(ao_num)
+ do i = 1, ao_num
+  list(i) = maxval(List_comb_thr_b3_size(:,i))
+ enddo
+ max_List_comb_thr_b3_size = maxval(list) 
+ print*,'max_List_comb_thr_b3_size =  ',max_List_comb_thr_b3_size
+ 
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision, List_comb_thr_b3_coef, (   max_List_comb_thr_b3_size,ao_num, ao_num )]
+&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_expo, (   max_List_comb_thr_b3_size,ao_num, ao_num )]
+&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_cent, (3, max_List_comb_thr_b3_size,ao_num, ao_num )]
+&BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)]
+ implicit none
+ integer :: i_1s,i,j,ipoint,icount
+ double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: r(3),weight,dist
+ thr = 1.d-15
+ ao_abs_comb_b3_j1b = 10000000.d0
+ do i = 1, ao_num
+  do j = 1, ao_num
+   icount = 0
+   do i_1s = 1, List_all_comb_b3_size
+     coef        = List_all_comb_b3_coef  (i_1s)
+     beta        = List_all_comb_b3_expo  (i_1s)
+     beta = max(beta,1.d-12)
+     center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
+     if(dabs(coef).lt.thr)cycle
+     int_j1b = 0.d0
+     do ipoint = 1, n_points_extra_final_grid
+      r(1:3) = final_grid_points_extra(1:3,ipoint)
+      weight = final_weight_at_r_vector_extra(ipoint)
+      dist  = ( center(1) - r(1) )*( center(1) - r(1) )
+      dist += ( center(2) - r(2) )*( center(2) - r(2) )
+      dist += ( center(3) - r(3) )*( center(3) - r(3) )
+      int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
+     enddo
+     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+      icount += 1
+      List_comb_thr_b3_coef(icount,j,i) = coef
+      List_comb_thr_b3_expo(icount,j,i) = beta
+      List_comb_thr_b3_cent(1:3,icount,j,i) = center(1:3)
+      ao_abs_comb_b3_j1b(icount,j,i) = int_j1b
+     endif
+   enddo
+  enddo 
+ enddo
+
+! do i = 1, ao_num
+!  do j = 1, i-1
+!    do icount = 1, List_comb_thr_b3_size(j,i)
+!     List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j)
+!     List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j)
+!     List_comb_thr_b3_cent(1:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j)
+!    enddo
+!  enddo
+! enddo
+ 
+END_PROVIDER 
+

src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f

@@ -0,0 +1,195 @@
+double precision function NAI_pol_mult_erf_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu)
+  BEGIN_DOC
+  ! Computes the following integral R^3 :
+  !
+  ! .. math::
+  ! 
+  !   \int dr  (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !   \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$ exp(-delta (r - D)^2 ).
+  !
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), delta  ! pure gaussian "D" 
+ double precision, intent(in)    :: C_center(3),mu      ! coulomb center "C" and "mu" in the erf(mu*x)/x function
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+
+ double precision  :: NAI_pol_mult_erf
+ ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+ double precision  :: A_new(0:max_dim,3)! new polynom 
+ double precision  :: A_center_new(3)   ! new center
+ integer           :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+ double precision  :: alpha_new         ! new exponent
+ double precision  :: fact_a_new        ! constant factor
+ double precision  :: accu,coefx,coefy,coefz,coefxy,coefxyz,thr
+ integer           :: d(3),i,lx,ly,lz,iorder_tmp(3)
+ thr = 1.d-10
+ d = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+
+ ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order                                
+ call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new , & 
+                                      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+ ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+ accu = 0.d0
+ do lx = 0, iorder_a_new(1)
+  coefx = A_new(lx,1)
+  if(dabs(coefx).lt.thr)cycle
+  iorder_tmp(1) = lx
+  do ly = 0, iorder_a_new(2)
+   coefy = A_new(ly,2)
+   coefxy = coefx * coefy 
+   if(dabs(coefxy).lt.thr)cycle
+   iorder_tmp(2) = ly
+   do lz = 0, iorder_a_new(3)
+    coefz = A_new(lz,3)
+    coefxyz = coefxy * coefz 
+    if(dabs(coefxyz).lt.thr)cycle
+    iorder_tmp(3) = lz
+    accu += coefxyz * NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B,alpha_new,beta,C_center,n_pt_max_integrals,mu)
+   enddo
+  enddo
+ enddo
+ NAI_pol_mult_erf_gauss_r12 = fact_a_new * accu 
+end
+
+subroutine erfc_mu_gauss_xyz(D_center,delta,mu,A_center,B_center,power_A,power_B,alpha,beta,n_pt_in,xyz_ints)
+  BEGIN_DOC
+  ! Computes the following integral :
+  !
+  ! .. math::
+  ! 
+  !   \int dr exp(-delta (r - D)^2 ) x/y/z * (1 - erf(mu |r-r'|))/ |r-r'| * (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  !   xyz_ints(1) = x , xyz_ints(2) = y, xyz_ints(3) = z, xyz_ints(4) = x^0 
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), delta,mu  ! pure gaussian "D" and mu parameter
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3),n_pt_in
+ double precision, intent(out)   :: xyz_ints(4)
+
+ double precision  :: NAI_pol_mult_erf
+ ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+ double precision  :: A_new(0:max_dim,3)! new polynom 
+ double precision  :: A_center_new(3)   ! new center
+ integer           :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+ double precision  :: alpha_new         ! new exponent
+ double precision  :: fact_a_new        ! constant factor
+ double precision  :: accu,coefx,coefy,coefz,coefxy,coefxyz,thr,contrib,contrib_inf,mu_inf
+ integer           :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1,mm
+ integer           :: power_B_tmp(3)
+ dim1=100
+ mu_inf = 1.d+10
+ thr = 1.d-10
+ d = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+
+ ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order                                
+ call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new , & 
+                                      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+ ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+ xyz_ints = 0.d0
+ do lx = 0, iorder_a_new(1)
+  coefx = A_new(lx,1)
+  if(dabs(coefx).lt.thr)cycle
+  iorder_tmp(1) = lx
+  do ly = 0, iorder_a_new(2)
+   coefy = A_new(ly,2)
+   coefxy = coefx * coefy 
+   if(dabs(coefxy).lt.thr)cycle
+   iorder_tmp(2) = ly
+   do lz = 0, iorder_a_new(3)
+    coefz = A_new(lz,3)
+    coefxyz = coefxy * coefz 
+    if(dabs(coefxyz).lt.thr)cycle
+    iorder_tmp(3) = lz
+     power_B_tmp = power_B
+     contrib = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B_tmp,alpha_new,beta,D_center,n_pt_in,mu)  
+     contrib_inf = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B_tmp,alpha_new,beta,D_center,n_pt_in,mu_inf)  
+     xyz_ints(4) += (contrib_inf - contrib) * coefxyz ! usual term with no x/y/z 
+                                      
+     do mm = 1, 3 
+      ! (x phi_i ) * phi_j 
+      ! x * (x - B_x)^b_x = B_x (x - B_x)^b_x + 1 * (x - B_x)^{b_x+1}
+      
+      !
+      ! first contribution :: B_x (x - B_x)^b_x :: usual integral multiplied by B_x
+      power_B_tmp = power_B
+      contrib_inf = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B_tmp,alpha_new,beta,D_center,n_pt_in,mu_inf)  
+      contrib = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B_tmp,alpha_new,beta,D_center,n_pt_in,mu)  
+      xyz_ints(mm) += (contrib_inf - contrib) * B_center(mm) * coefxyz 
+                                                       
+      !
+      ! second contribution :: (x - B_x)^(b_x+1) :: integral with b_x=>b_x+1 
+      power_B_tmp(mm) += 1
+      contrib = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B_tmp,alpha_new,beta,D_center,n_pt_in,mu)  
+      contrib_inf = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B_tmp,alpha_new,beta,D_center,n_pt_in,mu_inf)  
+      xyz_ints(mm) += (contrib_inf -  contrib) * coefxyz     
+     enddo
+   enddo
+  enddo
+ enddo
+ xyz_ints *= fact_a_new 
+end
+
+
+double precision function erf_mu_gauss(D_center,delta,mu,A_center,B_center,power_A,power_B,alpha,beta,n_pt_in)
+  BEGIN_DOC
+  ! Computes the following integral :
+  !
+  ! .. math::
+  ! 
+  !   \int dr exp(-delta (r - D)^2 ) erf(mu*|r-r'|)/ |r-r'| * (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), delta,mu  ! pure gaussian "D" and mu parameter
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3),n_pt_in
+
+ double precision  :: NAI_pol_mult_erf
+ ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+ double precision  :: A_new(0:max_dim,3)! new polynom 
+ double precision  :: A_center_new(3)   ! new center
+ integer           :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+ double precision  :: alpha_new         ! new exponent
+ double precision  :: fact_a_new        ! constant factor
+ double precision  :: accu,coefx,coefy,coefz,coefxy,coefxyz,thr,contrib,contrib_inf,mu_inf
+ integer           :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1,mm
+ dim1=100
+ mu_inf = 1.d+10
+ thr = 1.d-10
+ d = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+
+ ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order                                
+ call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new , & 
+                                      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+ ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+ erf_mu_gauss = 0.d0
+ do lx = 0, iorder_a_new(1)
+  coefx = A_new(lx,1)
+  if(dabs(coefx).lt.thr)cycle
+  iorder_tmp(1) = lx
+  do ly = 0, iorder_a_new(2)
+   coefy = A_new(ly,2)
+   coefxy = coefx * coefy 
+   if(dabs(coefxy).lt.thr)cycle
+   iorder_tmp(2) = ly
+   do lz = 0, iorder_a_new(3)
+    coefz = A_new(lz,3)
+    coefxyz = coefxy * coefz 
+    if(dabs(coefxyz).lt.thr)cycle
+    iorder_tmp(3) = lz
+    contrib = NAI_pol_mult_erf(A_center_new,B_center,iorder_tmp,power_B,alpha_new,beta,D_center,n_pt_in,mu)
+    erf_mu_gauss += contrib * coefxyz     
+   enddo
+  enddo
+ enddo
+ erf_mu_gauss *= fact_a_new 
+end
+

src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f

@@ -0,0 +1,340 @@
+! ---
+
+double precision function overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
+
+  BEGIN_DOC
+  !
+  ! Computes the following integral :
+  !
+  ! .. math                      ::
+  !
+  !   \int dr exp(-delta (r - D)^2 ) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  END_DOC
+
+  include 'constants.include.F'
+
+  implicit none
+  double precision, intent(in) :: D_center(3), delta  ! pure gaussian "D"
+  double precision, intent(in) :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+  integer, intent(in)          :: power_A(3),power_B(3)
+
+  double precision             :: overlap_x,overlap_y,overlap_z,overlap
+  ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+  double precision             :: A_new(0:max_dim,3)! new polynom
+  double precision             :: A_center_new(3)   ! new center
+  integer                      :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+  double precision             :: alpha_new         ! new exponent
+  double precision             :: fact_a_new        ! constant factor
+  double precision             :: accu, coefx, coefy, coefz, coefxy, coefxyz, thr
+  integer                      :: d(3), i, lx, ly, lz, iorder_tmp(3), dim1
+
+  dim1 = 100
+  thr  = 1.d-10
+  d(:) = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+  overlap_gauss_r12 = 0.d0
+
+  ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order
+  call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new ,&
+      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+  if(fact_a_new.lt.thr)return
+  ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+  accu = 0.d0
+  do lx = 0, iorder_a_new(1)
+    coefx = A_new(lx,1)*fact_a_new
+    if(dabs(coefx).lt.thr)cycle
+    iorder_tmp(1) = lx
+
+    do ly = 0, iorder_a_new(2)
+      coefy  = A_new(ly,2)
+      coefxy = coefx * coefy
+      if(dabs(coefxy) .lt. thr) cycle
+      iorder_tmp(2) = ly
+
+      do lz = 0, iorder_a_new(3)
+        coefz   = A_new(lz,3)
+        coefxyz = coefxy * coefz
+        if(dabs(coefxyz) .lt. thr) cycle
+        iorder_tmp(3) = lz
+
+        call overlap_gaussian_xyz( A_center_new, B_center, alpha_new, beta, iorder_tmp, power_B &
+                                 , overlap_x, overlap_y, overlap_z, overlap, dim1)
+
+        accu += coefxyz * overlap
+      enddo
+    enddo
+  enddo
+  overlap_gauss_r12 = accu
+end
+
+!---
+double precision function overlap_abs_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta)
+  BEGIN_DOC
+  ! Computes the following integral :
+  !
+  ! .. math                      ::
+  !
+  !   \int dr exp(-delta (r - D)^2 ) |(x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )|
+  !
+  END_DOC
+
+  implicit none
+  include 'constants.include.F'
+  double precision, intent(in)   :: D_center(3), delta  ! pure gaussian "D"
+  double precision, intent(in)   :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+  integer, intent(in)            :: power_A(3),power_B(3)
+
+  double precision               :: overlap_x,overlap_y,overlap_z,overlap
+  ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+  double precision               :: A_new(0:max_dim,3)! new polynom
+  double precision               :: A_center_new(3)   ! new center
+  integer                        :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+  double precision               :: alpha_new         ! new exponent
+  double precision               :: fact_a_new        ! constant factor
+  double precision               :: accu,coefx,coefy,coefz,coefxy,coefxyz,thr,dx,lower_exp_val
+  integer                        :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1
+  dim1=50
+  lower_exp_val = 40.d0
+  thr = 1.d-12
+  d(:) = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+  overlap_abs_gauss_r12 = 0.d0
+
+  ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order
+  call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new ,&
+      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+  if(fact_a_new.lt.thr)return
+  ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+  accu = 0.d0
+  do lx = 0, iorder_a_new(1)
+    coefx = A_new(lx,1)*fact_a_new
+!    if(dabs(coefx).lt.thr)cycle
+    iorder_tmp(1) = lx
+    do ly = 0, iorder_a_new(2)
+      coefy = A_new(ly,2)
+      coefxy = coefx * coefy
+      if(dabs(coefxy).lt.thr)cycle
+      iorder_tmp(2) = ly
+      do lz = 0, iorder_a_new(3)
+        coefz = A_new(lz,3)
+        coefxyz = coefxy * coefz
+        if(dabs(coefxyz).lt.thr)cycle
+        iorder_tmp(3) = lz
+        call overlap_x_abs(A_center_new(1),B_center(1),alpha_new,beta,iorder_tmp(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)
+        call overlap_x_abs(A_center_new(2),B_center(2),alpha_new,beta,iorder_tmp(2),power_B(2),overlap_y,lower_exp_val,dx,dim1)
+        call overlap_x_abs(A_center_new(3),B_center(3),alpha_new,beta,iorder_tmp(3),power_B(3),overlap_z,lower_exp_val,dx,dim1)
+        accu += dabs(coefxyz * overlap_x * overlap_y * overlap_z)
+      enddo
+    enddo
+  enddo
+  overlap_abs_gauss_r12= accu
+end
+
+!---
+
+! TODO apply Gaussian product three times first
+subroutine overlap_gauss_r12_v(D_center, LD_D, delta, A_center, B_center, power_A, power_B, alpha, beta, rvec, LD_rvec, n_points)
+
+  BEGIN_DOC
+  !
+  ! Computes the following integral :
+  !
+  !   \int dr exp(-delta (r - D)^2) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2)
+  !   using an array of D_centers
+  !
+  ! n_points: nb of integrals
+  !
+  END_DOC
+
+  implicit none
+
+  include 'constants.include.F'
+
+  integer,          intent(in)  :: LD_D, LD_rvec, n_points
+  integer,          intent(in)  :: power_A(3), power_B(3)
+  double precision, intent(in)  :: D_center(LD_D,3), delta
+  double precision, intent(in)  :: A_center(3), B_center(3), alpha, beta 
+  double precision, intent(out) :: rvec(LD_rvec)
+
+  integer                       :: maxab
+  integer                       :: d(3), i, lx, ly, lz, iorder_tmp(3), ipoint
+  double precision              :: overlap_x, overlap_y, overlap_z
+  double precision              :: alpha_new
+  double precision              :: accu, thr, coefxy
+  integer,          allocatable :: iorder_a_new(:)
+  double precision, allocatable :: overlap(:)
+  double precision, allocatable :: A_new(:,:,:), A_center_new(:,:)
+  double precision, allocatable :: fact_a_new(:)
+
+  thr  = 1.d-10
+  d(:) = 0
+
+  maxab = maxval(power_A(1:3))
+
+  allocate(A_new(n_points,0:maxab,3), A_center_new(n_points,3), fact_a_new(n_points), iorder_a_new(3), overlap(n_points))
+
+  call give_explicit_poly_and_gaussian_v(A_new, maxab, A_center_new, alpha_new, fact_a_new, iorder_a_new, delta, alpha, d, power_A, D_center, LD_D, A_center, n_points)
+
+  rvec(:) = 0.d0
+
+  do lx = 0, iorder_a_new(1)
+    iorder_tmp(1) = lx
+
+    do ly = 0, iorder_a_new(2)
+      iorder_tmp(2) = ly
+
+      do lz = 0, iorder_a_new(3)
+        iorder_tmp(3) = lz
+
+        call overlap_gaussian_xyz_v(A_center_new, B_center, alpha_new, beta, iorder_tmp, power_B, overlap, n_points)
+
+        do ipoint = 1, n_points
+          rvec(ipoint) = rvec(ipoint) + A_new(ipoint,lx,1) * A_new(ipoint,ly,2) * A_new(ipoint,lz,3) * overlap(ipoint)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  do ipoint = 1, n_points
+    rvec(ipoint) = rvec(ipoint) * fact_a_new(ipoint)
+  enddo
+
+  deallocate(A_new, A_center_new, fact_a_new, iorder_a_new, overlap)
+
+end subroutine overlap_gauss_r12_v
+
+!---
+
+subroutine overlap_gauss_xyz_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta, gauss_ints)
+
+  BEGIN_DOC
+  ! Computes the following integral :
+  !
+  ! .. math::
+  !
+  !   gauss_ints(m) = \int dr exp(-delta (r - D)^2 ) * x/y/z (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  ! with m == 1 ==> x, m == 2 ==> y, m == 3 ==> z
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), delta  ! pure gaussian "D"
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+ double precision, intent(out)   :: gauss_ints(3)
+
+ double precision  :: overlap_x,overlap_y,overlap_z,overlap
+ ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+ double precision  :: A_new(0:max_dim,3)! new polynom
+ double precision  :: A_center_new(3)   ! new center
+ integer           :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+ integer           :: power_B_new(3)
+ double precision  :: alpha_new         ! new exponent
+ double precision  :: fact_a_new        ! constant factor
+ double precision  :: coefx,coefy,coefz,coefxy,coefxyz,thr
+ integer           :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1,m
+ dim1=100
+ thr = 1.d-10
+ d = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+
+ ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order
+ call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new , &
+                                      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+ ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+ gauss_ints = 0.d0
+ do lx = 0, iorder_a_new(1)
+  coefx = A_new(lx,1)
+  if(dabs(coefx).lt.thr)cycle
+  iorder_tmp(1) = lx
+  do ly = 0, iorder_a_new(2)
+   coefy = A_new(ly,2)
+   coefxy = coefx * coefy
+   if(dabs(coefxy).lt.thr)cycle
+   iorder_tmp(2) = ly
+   do lz = 0, iorder_a_new(3)
+    coefz = A_new(lz,3)
+    coefxyz = coefxy * coefz
+    if(dabs(coefxyz).lt.thr)cycle
+    iorder_tmp(3) = lz
+    do m = 1, 3
+     ! change (x-Bx)^bx --> (x-Bx)^(bx+1) + Bx(x-Bx)^bx
+     power_B_new = power_B
+     power_B_new(m) += 1 ! (x-Bx)^(bx+1)
+     call overlap_gaussian_xyz(A_center_new,B_center,alpha_new,beta,iorder_tmp,power_B_new,overlap_x,overlap_y,overlap_z,overlap,dim1)
+     gauss_ints(m) += coefxyz * overlap
+
+     power_B_new = power_B
+     call overlap_gaussian_xyz(A_center_new,B_center,alpha_new,beta,iorder_tmp,power_B_new,overlap_x,overlap_y,overlap_z,overlap,dim1)
+     gauss_ints(m) += coefxyz * overlap * B_center(m) ! Bx (x-Bx)^(bx)
+    enddo
+   enddo
+  enddo
+ enddo
+ gauss_ints *= fact_a_new
+end
+
+double precision function overlap_gauss_xyz_r12_specific(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,mx)
+  BEGIN_DOC
+  ! Computes the following integral :
+  !
+  ! .. math::
+  !
+  !    \int dr exp(-delta (r - D)^2 ) * x/y/z (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  ! with mx == 1 ==> x, mx == 2 ==> y, mx == 3 ==> z
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), delta  ! pure gaussian "D"
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3),mx
+
+ double precision  :: overlap_x,overlap_y,overlap_z,overlap
+ ! First you multiply the usual gaussian "A" with the gaussian exp(-delta (r - D)^2 )
+ double precision  :: A_new(0:max_dim,3)! new polynom
+ double precision  :: A_center_new(3)   ! new center
+ integer           :: iorder_a_new(3)   ! i_order(i) = order of the new polynom ==> should be equal to power_A
+ integer           :: power_B_new(3)
+ double precision  :: alpha_new         ! new exponent
+ double precision  :: fact_a_new        ! constant factor
+ double precision  :: coefx,coefy,coefz,coefxy,coefxyz,thr
+ integer           :: d(3),i,lx,ly,lz,iorder_tmp(3),dim1,m
+ dim1=100
+ thr = 1.d-10
+ d = 0 ! order of the polynom for the gaussian exp(-delta (r - D)^2 )  == 0
+
+ ! New gaussian/polynom defined by :: new pol new center new expo   cst fact new order
+ call give_explicit_poly_and_gaussian(A_new , A_center_new , alpha_new, fact_a_new , iorder_a_new , &
+                                      delta,alpha,d,power_A,D_center,A_center,n_pt_max_integrals)
+ ! The new gaussian exp(-delta (r - D)^2 ) (x-A_x)^a \exp(-\alpha (x-A_x)^2
+ overlap_gauss_xyz_r12_specific = 0.d0
+ do lx = 0, iorder_a_new(1)
+  coefx = A_new(lx,1)
+  if(dabs(coefx).lt.thr)cycle
+  iorder_tmp(1) = lx
+  do ly = 0, iorder_a_new(2)
+   coefy = A_new(ly,2)
+   coefxy = coefx * coefy
+   if(dabs(coefxy).lt.thr)cycle
+   iorder_tmp(2) = ly
+   do lz = 0, iorder_a_new(3)
+    coefz = A_new(lz,3)
+    coefxyz = coefxy * coefz
+    if(dabs(coefxyz).lt.thr)cycle
+    iorder_tmp(3) = lz
+    m = mx
+    ! change (x-Bx)^bx --> (x-Bx)^(bx+1) + Bx(x-Bx)^bx
+    power_B_new = power_B
+    power_B_new(m) += 1 ! (x-Bx)^(bx+1)
+    call overlap_gaussian_xyz(A_center_new,B_center,alpha_new,beta,iorder_tmp,power_B_new,overlap_x,overlap_y,overlap_z,overlap,dim1)
+    overlap_gauss_xyz_r12_specific += coefxyz * overlap
+
+    power_B_new = power_B
+    call overlap_gaussian_xyz(A_center_new,B_center,alpha_new,beta,iorder_tmp,power_B_new,overlap_x,overlap_y,overlap_z,overlap,dim1)
+    overlap_gauss_xyz_r12_specific += coefxyz * overlap * B_center(m) ! Bx (x-Bx)^(bx)
+   enddo
+  enddo
+ enddo
+ overlap_gauss_xyz_r12_specific *= fact_a_new
+end

src/ao_many_one_e_ints/stg_gauss_int.irp.f

@@ -0,0 +1,121 @@
+double precision function ovlp_stg_gauss_int_phi_ij(D_center,gam,delta,A_center,B_center,power_A,power_B,alpha,beta)
+  BEGIN_DOC
+  ! Computes the following integral : 
+  !
+  ! .. math::
+  ! 
+  !   \int dr exp(-gam (r - D)) exp(-delta * (r -D)^2) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  END_DOC
+
+ implicit none
+ double precision, intent(in)    :: D_center(3), gam ! pure Slater "D" in r-r_D
+ double precision, intent(in)    :: delta            ! gaussian        in r-r_D
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+
+ integer :: i
+ double precision :: integral,gama_gauss
+ double precision, allocatable :: expos_slat(:)
+ allocate(expos_slat(n_max_fit_slat))
+ double precision :: overlap_gauss_r12
+ ovlp_stg_gauss_int_phi_ij = 0.d0
+ call expo_fit_slater_gam(gam,expos_slat)
+ do i = 1, n_max_fit_slat
+  gama_gauss = expos_slat(i)+delta 
+  integral = overlap_gauss_r12(D_center,gama_gauss,A_center,B_center,power_A,power_B,alpha,beta)
+  ovlp_stg_gauss_int_phi_ij += coef_fit_slat_gauss(i) * integral
+ enddo
+end
+
+
+double precision function erf_mu_stg_gauss_int_phi_ij(D_center,gam,delta,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu)
+  BEGIN_DOC
+  ! Computes the following integral : 
+  !
+  ! .. math::
+  ! 
+  !   \int dr exp(-gam(r - D)-delta(r - D)^2) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !   \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
+  !
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), gam ! pure Slater "D" in r-r_D
+ double precision, intent(in)    :: delta            ! gaussian        in r-r_D
+ double precision, intent(in)    :: C_center(3),mu      ! coulomb center "C" and "mu" in the erf(mu*x)/x function
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+
+ integer :: i
+ double precision :: NAI_pol_mult_erf_gauss_r12
+ double precision :: integral,gama_gauss
+ double precision, allocatable :: expos_slat(:)
+ allocate(expos_slat(n_max_fit_slat))
+ erf_mu_stg_gauss_int_phi_ij = 0.d0
+ call expo_fit_slater_gam(gam,expos_slat)
+ do i = 1, n_max_fit_slat
+  gama_gauss = expos_slat(i) + delta
+  integral = NAI_pol_mult_erf_gauss_r12(D_center,gama_gauss,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu)
+  erf_mu_stg_gauss_int_phi_ij += coef_fit_slat_gauss(i) * integral
+ enddo
+end
+
+double precision function overlap_stg_gauss(D_center,gam,A_center,B_center,power_A,power_B,alpha,beta)
+  BEGIN_DOC
+  ! Computes the following integral : 
+  !
+  ! .. math::
+  ! 
+  !   \int dr exp(-gam (r - D)) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  END_DOC
+
+ implicit none
+ double precision, intent(in)    :: D_center(3), gam ! pure Slater "D" 
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+ 
+ integer :: i
+ double precision :: expos_slat(n_max_fit_slat),integral,delta
+ double precision :: overlap_gauss_r12
+ overlap_stg_gauss = 0.d0
+ call expo_fit_slater_gam(gam,expos_slat)
+ do i = 1, n_max_fit_slat
+  delta = expos_slat(i) 
+  integral = overlap_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta)
+  overlap_stg_gauss += coef_fit_slat_gauss(i) * integral
+ enddo
+end
+
+double precision function erf_mu_stg_gauss(D_center,gam,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu)
+  BEGIN_DOC
+  ! Computes the following integral : 
+  !
+  ! .. math::
+  ! 
+  !   \int dr exp(-gam(r - D)) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !   \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
+  !
+  END_DOC
+
+ implicit none
+  include 'constants.include.F'
+ double precision, intent(in)    :: D_center(3), gam    ! pure Slater "D" 
+ double precision, intent(in)    :: C_center(3),mu      ! coulomb center "C" and "mu" in the erf(mu*x)/x function
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+
+ 
+ integer :: i
+ double precision :: expos_slat(n_max_fit_slat),integral,delta
+ double precision :: NAI_pol_mult_erf_gauss_r12
+ erf_mu_stg_gauss = 0.d0
+ call expo_fit_slater_gam(gam,expos_slat)
+ do i = 1, n_max_fit_slat
+  delta = expos_slat(i) 
+  integral = NAI_pol_mult_erf_gauss_r12(D_center,delta,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu)
+  erf_mu_stg_gauss += coef_fit_slat_gauss(i) * integral
+ enddo
+end

src/ao_many_one_e_ints/taylor_exp.irp.f

@@ -0,0 +1,101 @@
+double precision function exp_dl(x,n)
+ implicit none
+ double precision, intent(in) :: x
+ integer         , intent(in) :: n
+ integer :: i
+ exp_dl = 1.d0
+ do i = 1, n
+  exp_dl += fact_inv(i) * x**dble(i)
+ enddo
+end
+
+subroutine exp_dl_rout(x,n, array)
+ implicit none
+ double precision, intent(in) :: x
+ integer         , intent(in) :: n
+ double precision, intent(out)::  array(0:n)
+ integer :: i
+ double precision :: accu
+ accu = 1.d0
+ array(0) = 1.d0
+ do i = 1, n
+  accu += fact_inv(i) * x**dble(i)
+  array(i) = accu
+ enddo
+end
+
+subroutine exp_dl_ovlp_stg_phi_ij(zeta,D_center,gam,delta,A_center,B_center,power_A,power_B,alpha,beta,n_taylor,array_ints,integral_taylor,exponent_exp)
+  BEGIN_DOC
+  ! Computes the following integrals : 
+  !
+  ! .. math::
+  ! 
+  !   array(i) = \int dr EXP{exponent_exp * [exp(-gam*i (r - D)) exp(-delta*i * (r -D)^2)] (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !
+  !
+  ! and gives back the Taylor expansion of the exponential in integral_taylor
+  END_DOC
+
+ implicit none
+ double precision, intent(in)    :: zeta             ! prefactor of the argument of the exp(-zeta*x)
+ integer, intent(in)             :: n_taylor         ! order of the Taylor expansion of the exponential
+ double precision, intent(in)    :: D_center(3), gam ! pure Slater "D" in r-r_D
+ double precision, intent(in)    :: delta            ! gaussian        in r-r_D
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ double precision, intent(in)    :: exponent_exp
+ integer, intent(in)             :: power_A(3),power_B(3)
+ double precision, intent(out)   :: array_ints(0:n_taylor),integral_taylor
+
+ integer :: i,dim1
+ double precision :: delta_exp,gam_exp,ovlp_stg_gauss_int_phi_ij
+ double precision :: overlap_x,overlap_y,overlap_z,overlap
+ dim1=100
+ call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
+ array_ints(0) = overlap
+ integral_taylor = array_ints(0)
+ do i = 1, n_taylor
+  delta_exp = dble(i) * delta
+  gam_exp   = dble(i) * gam
+  array_ints(i) = ovlp_stg_gauss_int_phi_ij(D_center,gam_exp,delta_exp,A_center,B_center,power_A,power_B,alpha,beta)
+  integral_taylor += (-zeta*exponent_exp)**dble(i) * fact_inv(i) * array_ints(i)
+ enddo
+
+end
+
+subroutine exp_dl_erf_stg_phi_ij(zeta,D_center,gam,delta,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu,n_taylor,array_ints,integral_taylor)
+  BEGIN_DOC
+  ! Computes the following integrals : 
+  !
+  ! .. math::
+  ! 
+  !   array(i) = \int dr exp(-gam*i (r - D)) exp(-delta*i * (r -D)^2) (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
+  !   \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
+  !
+  !
+  ! and gives back the Taylor expansion of the exponential in integral_taylor
+  END_DOC
+
+ implicit none
+ integer, intent(in)             :: n_taylor         ! order of the Taylor expansion of the exponential
+ double precision, intent(in)    :: zeta             ! prefactor of the argument of the exp(-zeta*x)
+ double precision, intent(in)    :: D_center(3), gam ! pure Slater "D" in r-r_D
+ double precision, intent(in)    :: delta            ! gaussian        in r-r_D
+ double precision, intent(in)    :: C_center(3),mu      ! coulomb center "C" and "mu" in the erf(mu*x)/x function
+ double precision, intent(in)    :: A_center(3),B_center(3),alpha,beta ! gaussian/polynoms "A" and "B"
+ integer, intent(in)             :: power_A(3),power_B(3)
+ double precision, intent(out)   :: array_ints(0:n_taylor),integral_taylor
+
+ integer :: i,dim1
+ double precision :: delta_exp,gam_exp,NAI_pol_mult_erf,erf_mu_stg_gauss_int_phi_ij
+ dim1=100
+ 
+ array_ints(0) = NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_max_integrals,mu)
+ integral_taylor = array_ints(0)
+ do i = 1, n_taylor
+  delta_exp = dble(i) * delta
+  gam_exp   = dble(i) * gam
+  array_ints(i) = erf_mu_stg_gauss_int_phi_ij(D_center,gam_exp,delta_exp,A_center,B_center,power_A,power_B,alpha,beta,C_center,mu)
+  integral_taylor += (-zeta)**dble(i) * fact_inv(i) * array_ints(i)
+ enddo
+
+end

src/ao_many_one_e_ints/xyz_grad_xyz_ao_pol.irp.f

@@ -0,0 +1,343 @@
+ BEGIN_PROVIDER [double precision, coef_xyz_ao, (2,3,ao_num)]
+&BEGIN_PROVIDER [integer, power_xyz_ao, (2,3,ao_num)]
+ implicit none
+ BEGIN_DOC
+! coefficient for the basis function :: (x * phi_i(r), y * phi_i(r), * z_phi(r))
+!
+! x * (x - A_x)^a_x = A_x (x - A_x)^a_x + 1 * (x - A_x)^{a_x+1}
+ END_DOC
+ integer :: i,j,k,num_ao,power_ao(1:3)
+ double precision :: center_ao(1:3)
+ do i = 1, ao_num
+  power_ao(1:3)= ao_power(i,1:3) 
+  num_ao = ao_nucl(i)
+  center_ao(1:3) = nucl_coord(num_ao,1:3)
+  do j = 1, 3
+   coef_xyz_ao(1,j,i) = center_ao(j) ! A_x (x - A_x)^a_x
+   power_xyz_ao(1,j,i)= power_ao(j)
+   coef_xyz_ao(2,j,i) = 1.d0         ! 1 * (x - A_x)^a_{x+1}
+   power_xyz_ao(2,j,i)= power_ao(j) + 1
+  enddo
+ enddo
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision, ao_coef_ord_grad_transp, (2,3,ao_prim_num_max,ao_num) ]
+&BEGIN_PROVIDER [ integer, power_ord_grad_transp, (2,3,ao_num) ]
+  implicit none
+  BEGIN_DOC
+  ! grad AO in terms of polynoms and coefficients 
+  ! 
+  ! WARNING !!!! SOME polynoms might be negative !!!!! 
+  !
+  ! WHEN IT IS THE CASE, coefficients are ZERO 
+  END_DOC
+  integer                        :: i,j,power_ao(3), m,kk
+  do j=1, ao_num
+    power_ao(1:3)= ao_power(j,1:3) 
+    do m = 1, 3
+     power_ord_grad_transp(1,m,j) = power_ao(m) - 1
+     power_ord_grad_transp(2,m,j) = power_ao(m) + 1
+    enddo
+    do i=1, ao_prim_num_max
+     do m = 1, 3
+       ao_coef_ord_grad_transp(1,m,i,j) = ao_coef_normalized_ordered(j,i) * dble(power_ao(m)) ! a_x * c_i 
+       ao_coef_ord_grad_transp(2,m,i,j) = -2.d0 * ao_coef_normalized_ordered(j,i) * ao_expo_ordered_transp(i,j) ! -2 * c_i * alpha_i 
+       do kk = 1, 2
+        if(power_ord_grad_transp(kk,m,j).lt.0)then
+         ao_coef_ord_grad_transp(kk,m,i,j) = 0.d0
+        endif
+       enddo
+     enddo
+    enddo
+  enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_coef_ord_xyz_grad_transp, (4,3,ao_prim_num_max,ao_num) ]
+&BEGIN_PROVIDER [ integer, power_ord_xyz_grad_transp, (4,3,ao_num) ]
+  implicit none
+  BEGIN_DOC
+  ! x * d/dx of an AO in terms of polynoms and coefficients 
+  !
+  ! WARNING !!!! SOME polynoms might be negative !!!!! 
+  !
+  ! WHEN IT IS THE CASE, coefficients are ZERO 
+  END_DOC
+  integer                        :: i,j,power_ao(3), m,num_ao,kk
+  double precision :: center_ao(1:3)
+  do j=1, ao_num
+   power_ao(1:3)= ao_power(j,1:3) 
+   num_ao = ao_nucl(j)
+   center_ao(1:3) = nucl_coord(num_ao,1:3)
+   do m = 1, 3
+     power_ord_xyz_grad_transp(1,m,j)   = power_ao(m) - 1
+     power_ord_xyz_grad_transp(2,m,j)   = power_ao(m)
+     power_ord_xyz_grad_transp(3,m,j)   = power_ao(m) + 1
+     power_ord_xyz_grad_transp(4,m,j)   = power_ao(m) + 2
+     do kk = 1, 4
+      if(power_ord_xyz_grad_transp(kk,m,j).lt.0)then
+       power_ord_xyz_grad_transp(kk,m,j) = -1
+      endif
+     enddo
+   enddo
+   do i=1, ao_prim_num_max
+    do m = 1, 3
+     ao_coef_ord_xyz_grad_transp(1,m,i,j) = dble(power_ao(m)) * ao_coef_normalized_ordered(j,i) * center_ao(m)
+     ao_coef_ord_xyz_grad_transp(2,m,i,j) = dble(power_ao(m)) * ao_coef_normalized_ordered(j,i) 
+     ao_coef_ord_xyz_grad_transp(3,m,i,j) = -2.d0 * ao_coef_normalized_ordered(j,i) * ao_expo_ordered_transp(i,j) * center_ao(m)
+     ao_coef_ord_xyz_grad_transp(4,m,i,j) = -2.d0 * ao_coef_normalized_ordered(j,i) * ao_expo_ordered_transp(i,j) 
+     do kk = 1, 4
+      if(power_ord_xyz_grad_transp(kk,m,j).lt.0)then
+       ao_coef_ord_xyz_grad_transp(kk,m,i,j) = 0.d0
+      endif
+     enddo
+    enddo
+   enddo
+  enddo
+
+END_PROVIDER
+
+subroutine xyz_grad_phi_ao(r,i_ao,xyz_grad_phi)
+ implicit none
+ integer, intent(in) :: i_ao
+ double precision, intent(in) :: r(3)
+ double precision, intent(out):: xyz_grad_phi(3) ! x * d/dx phi i, y * d/dy phi_i, z * d/dz phi_
+ double precision :: center_ao(3),beta
+ double precision :: accu(3,4),dr(3),r2,pol_usual(3)
+ integer :: m,power_ao(3),num_ao,j_prim
+ power_ao(1:3)= ao_power(i_ao,1:3) 
+ num_ao = ao_nucl(i_ao)
+ center_ao(1:3) = nucl_coord(num_ao,1:3)
+ dr(1) = (r(1) - center_ao(1))
+ dr(2) = (r(2) - center_ao(2))
+ dr(3) = (r(3) - center_ao(3))
+ r2 = 0.d0
+ do m = 1, 3
+  r2 += dr(m)*dr(m)
+ enddo
+ ! computes the gaussian part 
+ accu = 0.d0
+ do j_prim =1,ao_prim_num(i_ao)
+   beta = ao_expo_ordered_transp(j_prim,i_ao)
+   if(dabs(beta*r2).gt.50.d0)cycle
+   do m = 1, 3
+    accu(m,1) += ao_coef_ord_xyz_grad_transp(1,m,j_prim,i_ao) * dexp(-beta*r2) 
+    accu(m,2) += ao_coef_ord_xyz_grad_transp(2,m,j_prim,i_ao) * dexp(-beta*r2) 
+    accu(m,3) += ao_coef_ord_xyz_grad_transp(3,m,j_prim,i_ao) * dexp(-beta*r2) 
+    accu(m,4) += ao_coef_ord_xyz_grad_transp(4,m,j_prim,i_ao) * dexp(-beta*r2) 
+   enddo
+ enddo
+ ! computes the polynom part
+ pol_usual = 0.d0
+ pol_usual(1) = dr(2)**dble(power_ao(2)) * dr(3)**dble(power_ao(3)) 
+ pol_usual(2) = dr(1)**dble(power_ao(1)) * dr(3)**dble(power_ao(3)) 
+ pol_usual(3) = dr(1)**dble(power_ao(1)) * dr(2)**dble(power_ao(2)) 
+
+ xyz_grad_phi = 0.d0
+ do m = 1, 3
+  xyz_grad_phi(m) += accu(m,2) * pol_usual(m) * dr(m)**dble(power_ord_xyz_grad_transp(2,m,i_ao))
+  xyz_grad_phi(m) += accu(m,3) * pol_usual(m) * dr(m)**dble(power_ord_xyz_grad_transp(3,m,i_ao))
+  xyz_grad_phi(m) += accu(m,4) * pol_usual(m) * dr(m)**dble(power_ord_xyz_grad_transp(4,m,i_ao))
+  if(power_ord_xyz_grad_transp(1,m,i_ao).lt.0)cycle
+  xyz_grad_phi(m) += accu(m,1) * pol_usual(m) * dr(m)**dble(power_ord_xyz_grad_transp(1,m,i_ao))
+ enddo
+end
+
+subroutine grad_phi_ao(r,i_ao,grad_xyz_phi)
+ implicit none
+ integer, intent(in) :: i_ao
+ double precision, intent(in) :: r(3)
+ double precision, intent(out):: grad_xyz_phi(3) ! x * phi i, y * phi_i, z * phi_
+ double precision :: center_ao(3),beta
+ double precision :: accu(3,2),dr(3),r2,pol_usual(3)
+ integer :: m,power_ao(3),num_ao,j_prim
+ power_ao(1:3)= ao_power(i_ao,1:3) 
+ num_ao = ao_nucl(i_ao)
+ center_ao(1:3) = nucl_coord(num_ao,1:3)
+ dr(1) = (r(1) - center_ao(1))
+ dr(2) = (r(2) - center_ao(2))
+ dr(3) = (r(3) - center_ao(3))
+ r2 = 0.d0
+ do m = 1, 3
+  r2 += dr(m)*dr(m)
+ enddo
+ ! computes the gaussian part 
+ accu = 0.d0
+ do j_prim =1,ao_prim_num(i_ao)
+   beta = ao_expo_ordered_transp(j_prim,i_ao)
+   if(dabs(beta*r2).gt.50.d0)cycle
+   do m = 1, 3
+    accu(m,1) += ao_coef_ord_grad_transp(1,m,j_prim,i_ao) * dexp(-beta*r2) 
+    accu(m,2) += ao_coef_ord_grad_transp(2,m,j_prim,i_ao) * dexp(-beta*r2) 
+   enddo
+ enddo
+ ! computes the polynom part
+ pol_usual = 0.d0
+ pol_usual(1) = dr(2)**dble(power_ao(2)) * dr(3)**dble(power_ao(3)) 
+ pol_usual(2) = dr(1)**dble(power_ao(1)) * dr(3)**dble(power_ao(3)) 
+ pol_usual(3) = dr(1)**dble(power_ao(1)) * dr(2)**dble(power_ao(2)) 
+ do m = 1, 3
+  grad_xyz_phi(m)  = accu(m,2) * pol_usual(m) * dr(m)**dble(power_ord_grad_transp(2,m,i_ao))
+  if(power_ao(m)==0)cycle
+  grad_xyz_phi(m) += accu(m,1) * pol_usual(m) * dr(m)**dble(power_ord_grad_transp(1,m,i_ao))
+ enddo
+end
+
+subroutine xyz_phi_ao(r,i_ao,xyz_phi)
+ implicit none
+ integer, intent(in) :: i_ao
+ double precision, intent(in) :: r(3)
+ double precision, intent(out):: xyz_phi(3) ! x * phi i, y * phi_i, z * phi_i
+ double precision :: center_ao(3),beta
+ double precision :: accu,dr(3),r2,pol_usual(3)
+ integer :: m,power_ao(3),num_ao
+ power_ao(1:3)= ao_power(i_ao,1:3) 
+ num_ao = ao_nucl(i_ao)
+ center_ao(1:3) = nucl_coord(num_ao,1:3)
+ dr(1) = (r(1) - center_ao(1))
+ dr(2) = (r(2) - center_ao(2))
+ dr(3) = (r(3) - center_ao(3))
+ r2 = 0.d0
+ do m = 1, 3
+  r2 += dr(m)*dr(m)
+ enddo
+ ! computes the gaussian part 
+ accu = 0.d0
+ do m=1,ao_prim_num(i_ao)
+   beta = ao_expo_ordered_transp(m,i_ao)
+   if(dabs(beta*r2).gt.50.d0)cycle
+   accu += ao_coef_normalized_ordered_transp(m,i_ao) * dexp(-beta*r2)
+ enddo
+ ! computes the polynom part
+ pol_usual = 0.d0
+ pol_usual(1) = dr(2)**dble(power_ao(2)) * dr(3)**dble(power_ao(3)) 
+ pol_usual(2) = dr(1)**dble(power_ao(1)) * dr(3)**dble(power_ao(3)) 
+ pol_usual(3) = dr(1)**dble(power_ao(1)) * dr(2)**dble(power_ao(2)) 
+ do m = 1, 3
+  xyz_phi(m) = accu * pol_usual(m) * dr(m)**(dble(power_ao(m))) * ( coef_xyz_ao(1,m,i_ao) + coef_xyz_ao(2,m,i_ao) * dr(m) )
+ enddo
+end
+
+
+subroutine test_pol_xyz
+ implicit none
+ integer :: ipoint,i,j,m,jpoint
+ double precision :: r1(3),derf_mu_x
+ double precision :: weight1,r12,xyz_phi(3),grad_phi(3),xyz_grad_phi(3)
+ double precision, allocatable :: aos_array(:),aos_grad_array(:,:)
+ double precision :: num_xyz_phi(3),num_grad_phi(3),num_xyz_grad_phi(3)
+ double precision :: accu_xyz_phi(3),accu_grad_phi(3),accu_xyz_grad_phi(3)
+ double precision :: meta_accu_xyz_phi(3),meta_accu_grad_phi(3),meta_accu_xyz_grad_phi(3)
+ allocate(aos_array(ao_num),aos_grad_array(3,ao_num))
+ meta_accu_xyz_phi     = 0.d0
+ meta_accu_grad_phi    = 0.d0
+ meta_accu_xyz_grad_phi= 0.d0
+ do i = 1, ao_num
+  accu_xyz_phi     = 0.d0
+  accu_grad_phi    = 0.d0
+  accu_xyz_grad_phi= 0.d0
+
+  do ipoint = 1, n_points_final_grid
+   r1(:) = final_grid_points(:,ipoint)
+   weight1 = final_weight_at_r_vector(ipoint)
+   call give_all_aos_and_grad_at_r(r1,aos_array,aos_grad_array)
+   do m = 1, 3
+    num_xyz_phi(m)      = r1(m) *  aos_array(i)  
+    num_grad_phi(m)     = aos_grad_array(m,i)  
+    num_xyz_grad_phi(m) = r1(m) *  aos_grad_array(m,i) 
+   enddo
+   call xyz_phi_ao(r1,i,xyz_phi)
+   call grad_phi_ao(r1,i,grad_phi)
+   call xyz_grad_phi_ao(r1,i,xyz_grad_phi)
+   do m = 1, 3
+    accu_xyz_phi(m)      += weight1 * dabs(num_xyz_phi(m)      -  xyz_phi(m)     )
+    accu_grad_phi(m)     += weight1 * dabs(num_grad_phi(m)     -  grad_phi(m)    )
+    accu_xyz_grad_phi(m) += weight1 * dabs(num_xyz_grad_phi(m) -  xyz_grad_phi(m))
+   enddo
+  enddo
+  print*,''
+  print*,''
+  print*,'i,',i
+  print*,''
+  do m = 1, 3
+!    print*, 'm, accu_xyz_phi(m)  ' ,m, accu_xyz_phi(m)  
+!    print*, 'm, accu_grad_phi(m) ' ,m, accu_grad_phi(m)         
+    print*, 'm, accu_xyz_grad_phi' ,m, accu_xyz_grad_phi(m)
+  enddo
+  do m = 1, 3
+   meta_accu_xyz_phi(m) += dabs(accu_xyz_phi(m))
+   meta_accu_grad_phi(m) += dabs(accu_grad_phi(m))
+   meta_accu_xyz_grad_phi(m) += dabs(accu_xyz_grad_phi(m))
+  enddo
+ enddo
+  do m = 1, 3
+!    print*, 'm, meta_accu_xyz_phi(m)  ' ,m, meta_accu_xyz_phi(m)  
+!    print*, 'm, meta_accu_grad_phi(m) ' ,m, meta_accu_grad_phi(m)         
+    print*, 'm, meta_accu_xyz_grad_phi' ,m, meta_accu_xyz_grad_phi(m)
+  enddo
+
+
+
+end
+
+subroutine test_ints_semi_bis
+ implicit none
+ integer :: ipoint,i,j,m
+ double precision :: r1(3), aos_grad_array_r1(3, ao_num), aos_array_r1(ao_num)
+ double precision :: C_center(3), weight1,mu_in,r12,derf_mu_x,dxyz_ints(3),NAI_pol_mult_erf_ao
+ double precision :: ao_mat(ao_num,ao_num),ao_xmat(3,ao_num,ao_num),accu1, accu2(3)
+ mu_in = 0.5d0
+ C_center = 0.d0
+ C_center(1) = 0.25d0
+ C_center(3) = 1.12d0
+ C_center(2) = -1.d0
+ ao_mat = 0.d0
+ ao_xmat = 0.d0
+ do ipoint = 1, n_points_final_grid
+  r1(1) = final_grid_points(1,ipoint)
+  r1(2) = final_grid_points(2,ipoint)
+  r1(3) = final_grid_points(3,ipoint)
+  call give_all_aos_and_grad_at_r(r1,aos_array_r1,aos_grad_array_r1)
+  weight1 = final_weight_at_r_vector(ipoint)
+  r12 = (r1(1) - C_center(1))**2.d0 + (r1(2) - C_center(2))**2.d0 + (r1(3) - C_center(3))**2.d0 
+  r12 = dsqrt(r12)
+  do i = 1, ao_num
+   do j = 1, ao_num
+    ao_mat(j,i)  += aos_array_r1(i) * aos_array_r1(j) * weight1 * derf_mu_x(mu_in,r12)
+    do m = 1, 3
+     ao_xmat(m,j,i) += r1(m) * aos_array_r1(j) * aos_grad_array_r1(m,i) * weight1 * derf_mu_x(mu_in,r12)
+    enddo
+   enddo
+  enddo
+ enddo
+
+ accu1 = 0.d0
+ accu2 = 0.d0
+ accu1relat = 0.d0
+ accu2relat = 0.d0
+ double precision :: accu1relat, accu2relat(3)
+ double precision :: contrib(3)
+ do i = 1, ao_num
+  do j = 1, ao_num
+   call phi_j_erf_mu_r_xyz_dxyz_phi(i,j,mu_in, C_center, dxyz_ints)
+   print*,''
+   print*,'i,j',i,j
+   print*,dxyz_ints(:)
+   print*,ao_xmat(:,j,i)
+   do m = 1, 3
+    contrib(m) = dabs(ao_xmat(m,j,i) - dxyz_ints(m))
+    accu2(m) += contrib(m)
+    if(dabs(ao_xmat(m,j,i)).gt.1.d-10)then
+     accu2relat(m) += dabs(ao_xmat(m,j,i) - dxyz_ints(m))/dabs(ao_xmat(m,j,i))
+    endif
+   enddo
+    print*,contrib
+  enddo
+   print*,''
+ enddo
+ print*,'accu2relat = '
+ print*, accu2relat /dble(ao_num * ao_num)
+
+end
+
+

src/bi_ortho_mos/EZFIO.cfg

@@ -0,0 +1,11 @@
+[mo_r_coef]
+type: double precision
+doc: right-coefficient of the i-th |AO| on the j-th |MO|
+interface: ezfio
+size: (ao_basis.ao_num,mo_basis.mo_num)
+
+[mo_l_coef]
+type: double precision
+doc: right-coefficient of the i-th |AO| on the j-th |MO|
+interface: ezfio
+size: (ao_basis.ao_num,mo_basis.mo_num)

src/bi_ortho_mos/NEED

@@ -0,0 +1,3 @@
+mo_basis
+becke_numerical_grid
+dft_utils_in_r

src/bi_ortho_mos/bi_density.irp.f

@@ -0,0 +1,70 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, TCSCF_bi_ort_dm_ao_alpha, (ao_num, ao_num) ]
+
+  BEGIN_DOC
+  ! TCSCF_bi_ort_dm_ao_alpha(i,j) = <Chi_0| a^dagger_i,alpha a_j,alpha |Phi_0> where i,j are AO basis. 
+  !
+  ! This is the equivalent of the alpha density of the HF Slater determinant, but with a couple of bi-orthonormal Slater determinant |Chi_0> and |Phi_0>
+  END_DOC
+
+  implicit none
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0               &
+            , mo_l_coef, size(mo_l_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
+            !, mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+            , 0.d0, TCSCF_bi_ort_dm_ao_alpha, size(TCSCF_bi_ort_dm_ao_alpha, 1) )
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, TCSCF_bi_ort_dm_ao_beta, (ao_num, ao_num) ]
+
+  BEGIN_DOC
+  ! TCSCF_bi_ort_dm_ao_beta(i,j) = <Chi_0| a^dagger_i,beta a_j,beta |Phi_0> where i,j are AO basis. 
+  !
+  ! This is the equivalent of the beta density of the HF Slater determinant, but with a couple of bi-orthonormal Slater determinant |Chi_0> and |Phi_0>
+  END_DOC
+
+  implicit none
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0                &
+            , mo_l_coef, size(mo_l_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
+            !, mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+            , 0.d0, TCSCF_bi_ort_dm_ao_beta, size(TCSCF_bi_ort_dm_ao_beta, 1) )
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, TCSCF_bi_ort_dm_ao, (ao_num, ao_num) ]
+
+  BEGIN_DOC
+  ! TCSCF_bi_ort_dm_ao(i,j) = <Chi_0| a^dagger_i,beta+alpha a_j,beta+alpha |Phi_0> where i,j are AO basis. 
+  !
+  ! This is the equivalent of the total electronic density of the HF Slater determinant, but with a couple of bi-orthonormal Slater determinant |Chi_0> and |Phi_0>
+  END_DOC
+
+  implicit none
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  ASSERT(size(TCSCF_bi_ort_dm_ao, 1) == size(TCSCF_bi_ort_dm_ao_alpha, 1))
+
+  if(elec_alpha_num==elec_beta_num) then
+    TCSCF_bi_ort_dm_ao = TCSCF_bi_ort_dm_ao_alpha + TCSCF_bi_ort_dm_ao_alpha
+  else
+    ASSERT(size(TCSCF_bi_ort_dm_ao, 1) == size(TCSCF_bi_ort_dm_ao_beta, 1))
+    TCSCF_bi_ort_dm_ao = TCSCF_bi_ort_dm_ao_alpha + TCSCF_bi_ort_dm_ao_beta
+  endif
+
+END_PROVIDER
+
+! ---
+

src/bi_ortho_mos/bi_ort_mos_in_r.irp.f

@@ -0,0 +1,137 @@
+
+! TODO: left & right MO without duplicate AO calculation
+
+! ---
+
+BEGIN_PROVIDER[double precision, mos_r_in_r_array, (mo_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  ! mos_in_r_array(i,j) = value of the ith RIGHT mo on the jth grid point
+  END_DOC
+
+  implicit none
+  integer          :: i, j
+  double precision :: mos_array(mo_num), r(3)
+
+ !$OMP PARALLEL DO &
+ !$OMP DEFAULT (NONE)  &
+ !$OMP PRIVATE (i, j, r, mos_array) & 
+ !$OMP SHARED (mos_r_in_r_array, n_points_final_grid, mo_num, final_grid_points)
+  do i = 1, n_points_final_grid
+    r(1) = final_grid_points(1,i)
+    r(2) = final_grid_points(2,i)
+    r(3) = final_grid_points(3,i)
+    call give_all_mos_r_at_r(r, mos_array)
+    do j = 1, mo_num
+      mos_r_in_r_array(j,i) = mos_array(j)
+    enddo
+  enddo
+ !$OMP END PARALLEL DO
+ 
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER[double precision, mos_r_in_r_array_transp, (n_points_final_grid, mo_num)]
+
+  BEGIN_DOC
+  ! mos_r_in_r_array_transp(i,j) = value of the jth mo on the ith grid point
+  END_DOC
+
+  implicit none
+  integer :: i,j
+
+  do i = 1, n_points_final_grid
+    do j = 1, mo_num
+      mos_r_in_r_array_transp(i,j) = mos_r_in_r_array(j,i) 
+    enddo
+  enddo
+ 
+END_PROVIDER
+
+! ---
+
+subroutine give_all_mos_r_at_r(r, mos_r_array)
+
+  BEGIN_DOC
+  ! mos_r_array(i) = ith RIGHT MO function evaluated at "r"
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: mos_r_array(mo_num)
+  double precision              :: aos_array(ao_num)
+
+  call give_all_aos_at_r(r, aos_array)
+  call dgemv('N', mo_num, ao_num, 1.d0, mo_r_coef_transp, mo_num, aos_array, 1, 0.d0, mos_r_array, 1)
+
+end subroutine give_all_mos_r_at_r
+
+! ---
+
+BEGIN_PROVIDER[double precision, mos_l_in_r_array, (mo_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  ! mos_in_r_array(i,j) = value of the ith LEFT mo on the jth grid point
+  END_DOC
+
+  implicit none
+  integer          :: i, j
+  double precision :: mos_array(mo_num), r(3)
+
+ !$OMP PARALLEL DO &
+ !$OMP DEFAULT (NONE)  &
+ !$OMP PRIVATE (i,r,mos_array,j) & 
+ !$OMP SHARED(mos_l_in_r_array,n_points_final_grid,mo_num,final_grid_points)
+  do i = 1, n_points_final_grid
+    r(1) = final_grid_points(1,i)
+    r(2) = final_grid_points(2,i)
+    r(3) = final_grid_points(3,i)
+    call give_all_mos_l_at_r(r, mos_array)
+    do j = 1, mo_num
+      mos_l_in_r_array(j,i) = mos_array(j)
+    enddo
+  enddo
+ !$OMP END PARALLEL DO
+ 
+END_PROVIDER
+
+! ---
+
+subroutine give_all_mos_l_at_r(r, mos_l_array)
+
+  BEGIN_DOC
+  ! mos_l_array(i) = ith LEFT MO function evaluated at "r"
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: mos_l_array(mo_num)
+  double precision              :: aos_array(ao_num)
+
+  call give_all_aos_at_r(r, aos_array)
+  call dgemv('N', mo_num, ao_num, 1.d0, mo_l_coef_transp, mo_num, aos_array, 1, 0.d0, mos_l_array, 1)
+
+end subroutine give_all_mos_l_at_r
+
+! ---
+
+BEGIN_PROVIDER[double precision, mos_l_in_r_array_transp,(n_points_final_grid,mo_num)]
+
+  BEGIN_DOC
+  ! mos_l_in_r_array_transp(i,j) = value of the jth mo on the ith grid point
+  END_DOC
+
+  implicit none
+  integer :: i, j
+
+  do i = 1, n_points_final_grid
+    do j = 1, mo_num
+      mos_l_in_r_array_transp(i,j) = mos_l_in_r_array(j,i) 
+    enddo
+  enddo
+ 
+END_PROVIDER
+
+! ---
+

src/bi_ortho_mos/grad_bi_ort_mos_in_r.irp.f

@@ -0,0 +1,100 @@
+ BEGIN_PROVIDER[double precision, mos_r_grad_in_r_array,(mo_num,n_points_final_grid,3)]
+ implicit none
+ BEGIN_DOC
+ ! mos_r_grad_in_r_array(i,j,k)          = value of the kth component of the gradient of ith RIGHT mo on the jth grid point
+ !
+ ! k = 1 : x, k= 2, y, k  3, z
+ END_DOC
+ integer :: m
+ mos_r_grad_in_r_array = 0.d0
+ do m=1,3
+  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_r_coef_transp,mo_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_r_grad_in_r_array(1,1,m),mo_num)
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER[double precision, mos_r_grad_in_r_array_transp,(3,mo_num,n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! mos_r_grad_in_r_array_transp(i,j,k)   = value of the kth component of the gradient of jth RIGHT mo on the ith grid point
+ !
+ ! k = 1 : x, k= 2, y, k  3, z
+ END_DOC
+ integer :: m
+ integer  :: i,j
+ mos_r_grad_in_r_array_transp = 0.d0
+ do i = 1, n_points_final_grid
+  do j = 1, mo_num
+   do m = 1, 3
+     mos_r_grad_in_r_array_transp(m,j,i) = mos_r_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo                                                                                                                                                                                 
+ END_PROVIDER
+
+ BEGIN_PROVIDER[double precision, mos_r_grad_in_r_array_transp_bis,(3,n_points_final_grid,mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! mos_r_grad_in_r_array_transp(i,j,k)   = value of the ith component of the gradient on the jth grid point of jth RIGHT MO 
+ END_DOC
+ integer :: m
+ integer  :: i,j
+ mos_r_grad_in_r_array_transp_bis = 0.d0
+ do j = 1, mo_num
+  do i = 1, n_points_final_grid
+   do m = 1, 3
+     mos_r_grad_in_r_array_transp_bis(m,i,j) = mos_r_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo                                                                                                                                                                                 
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER[double precision, mos_l_grad_in_r_array,(mo_num,n_points_final_grid,3)]
+ implicit none
+ BEGIN_DOC
+ ! mos_l_grad_in_r_array(i,j,k)          = value of the kth component of the gradient of ith RIGHT mo on the jth grid point
+ !
+ ! k = 1 : x, k= 2, y, k  3, z
+ END_DOC
+ integer :: m
+ mos_l_grad_in_r_array = 0.d0
+ do m=1,3
+  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_r_coef_transp,mo_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_l_grad_in_r_array(1,1,m),mo_num)
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER[double precision, mos_l_grad_in_r_array_transp,(3,mo_num,n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! mos_l_grad_in_r_array_transp(i,j,k)   = value of the kth component of the gradient of jth RIGHT mo on the ith grid point
+ !
+ ! k = 1 : x, k= 2, y, k  3, z
+ END_DOC
+ integer :: m
+ integer  :: i,j
+ mos_l_grad_in_r_array_transp = 0.d0
+ do i = 1, n_points_final_grid
+  do j = 1, mo_num
+   do m = 1, 3
+     mos_l_grad_in_r_array_transp(m,j,i) = mos_l_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo                                                                                                                                                                                 
+ END_PROVIDER
+
+ BEGIN_PROVIDER[double precision, mos_l_grad_in_r_array_transp_bis,(3,n_points_final_grid,mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! mos_l_grad_in_r_array_transp(i,j,k)   = value of the ith component of the gradient on the jth grid point of jth RIGHT MO 
+ END_DOC
+ integer :: m
+ integer  :: i,j
+ mos_l_grad_in_r_array_transp_bis = 0.d0
+ do j = 1, mo_num
+  do i = 1, n_points_final_grid
+   do m = 1, 3
+     mos_l_grad_in_r_array_transp_bis(m,i,j) = mos_l_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo                                                                                                                                                                                 
+ END_PROVIDER

src/bi_ortho_mos/mos_rl.irp.f

@@ -0,0 +1,224 @@
+
+! ---
+
+subroutine ao_to_mo_bi_ortho(A_ao, LDA_ao, A_mo, LDA_mo)
+
+  BEGIN_DOC
+  !
+  ! Transform A from the |AO| basis to the BI ORTHONORMAL MOS 
+  !
+  ! $C_L^\dagger.A_{ao}.C_R$ where C_L and C_R are the LEFT and RIGHT MO coefs
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: LDA_ao, LDA_mo
+  double precision, intent(in)  :: A_ao(LDA_ao,ao_num)
+  double precision, intent(out) :: A_mo(LDA_mo,mo_num)
+  double precision, allocatable :: T(:,:)
+
+  allocate ( T(ao_num,mo_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+
+  ! T = A_ao x mo_r_coef
+  call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0      &
+            , A_ao, LDA_ao, mo_r_coef, size(mo_r_coef, 1) &
+            , 0.d0, T, size(T, 1) )
+
+  ! A_mo = mo_l_coef.T x T
+  call dgemm( 'T', 'N', mo_num, mo_num, ao_num, 1.d0       &
+            , mo_l_coef, size(mo_l_coef, 1), T, size(T, 1) &
+            , 0.d0, A_mo, LDA_mo )
+
+!  call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-12)
+  deallocate(T)
+
+end subroutine ao_to_mo_bi_ortho
+
+! ---
+
+subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao)
+
+  BEGIN_DOC
+  !
+  ! mo_l_coef.T x     A_ao   x mo_r_coef = A_mo
+  ! mo_l_coef.T x ao_overlap x mo_r_coef =  I
+  !
+  ! ==> A_ao = (ao_overlap x mo_r_coef) x A_mo x (ao_overlap x mo_l_coef).T
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: LDA_ao, LDA_mo
+  double precision, intent(in)  :: A_mo(LDA_mo,mo_num)
+  double precision, intent(out) :: A_ao(LDA_ao,ao_num)
+  double precision, allocatable :: tmp_1(:,:), tmp_2(:,:)
+
+  ! ao_overlap x mo_r_coef
+  allocate( tmp_1(ao_num,mo_num) )
+  call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0                         &
+            , ao_overlap, size(ao_overlap, 1), mo_r_coef, size(mo_r_coef, 1) &
+            , 0.d0, tmp_1, size(tmp_1, 1) )
+
+  ! (ao_overlap x mo_r_coef) x A_mo
+  allocate( tmp_2(ao_num,mo_num) )
+  call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
+            , tmp_1, size(tmp_1, 1), A_mo, LDA_mo    &
+            , 0.d0, tmp_2, size(tmp_2, 1) )
+  
+  ! ao_overlap x mo_l_coef
+  tmp_1 = 0.d0
+  call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0                         &
+            , ao_overlap, size(ao_overlap, 1), mo_l_coef, size(mo_l_coef, 1) &
+            , 0.d0, tmp_1, size(tmp_1, 1) )
+
+  ! (ao_overlap x mo_r_coef) x A_mo x (ao_overlap x mo_l_coef).T
+  call dgemm( 'N', 'T', ao_num, ao_num, mo_num, 1.d0       &
+            , tmp_2, size(tmp_2, 1), tmp_1, size(tmp_1, 1) &
+            , 0.d0, A_ao, LDA_ao )
+  
+  deallocate(tmp_1, tmp_2)
+
+end subroutine mo_to_ao_bi_ortho
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_r_coef, (ao_num, mo_num) ]
+
+  BEGIN_DOC
+  !
+  ! Molecular right-orbital coefficients on |AO| basis set
+  !
+  END_DOC
+
+  implicit none
+  integer :: i, j
+  logical :: exists
+
+  PROVIDE ezfio_filename
+
+  if (mpi_master) then
+    call ezfio_has_bi_ortho_mos_mo_r_coef(exists)
+  endif
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read mo_r_coef with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      call ezfio_get_bi_ortho_mos_mo_r_coef(mo_r_coef)
+      write(*,*) 'Read mo_r_coef'
+    endif
+    IRP_IF MPI
+      call MPI_BCAST(mo_r_coef, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read mo_r_coef with MPI'
+      endif
+    IRP_ENDIF
+  else
+
+    print*, 'mo_r_coef are mo_coef'
+    do i = 1, mo_num
+      do j = 1, ao_num
+        mo_r_coef(j,i) = mo_coef(j,i)
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_l_coef, (ao_num, mo_num) ]
+
+  BEGIN_DOC
+  !
+  ! Molecular left-orbital coefficients on |AO| basis set
+  !
+  END_DOC
+
+  implicit none
+  integer :: i, j
+  logical :: exists
+
+  PROVIDE ezfio_filename
+
+  if (mpi_master) then
+    call ezfio_has_bi_ortho_mos_mo_l_coef(exists)
+  endif
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read mo_l_coef with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      call ezfio_get_bi_ortho_mos_mo_l_coef(mo_l_coef)
+      write(*,*) 'Read mo_l_coef'
+    endif
+    IRP_IF MPI
+      call MPI_BCAST(mo_l_coef, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read mo_l_coef with MPI'
+      endif
+    IRP_ENDIF
+  else
+
+    print*, 'mo_l_coef are mo_coef'
+    do i = 1, mo_num
+      do j = 1, ao_num
+        mo_l_coef(j,i) = mo_coef(j,i)
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_r_coef_transp, (mo_num, ao_num)]
+
+  implicit none
+  integer :: j, m
+  do j = 1, mo_num
+    do m = 1, ao_num
+      mo_r_coef_transp(j,m) = mo_r_coef(m,j)
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_l_coef_transp, (mo_num, ao_num)]
+
+  implicit none
+  integer :: j, m
+  do j = 1, mo_num
+    do m = 1, ao_num
+      mo_l_coef_transp(j,m) = mo_l_coef(m,j)
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+

src/bi_ortho_mos/overlap.irp.f

@@ -0,0 +1,160 @@
+
+
+ BEGIN_PROVIDER [ double precision, overlap_bi_ortho, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, overlap_diag_bi_ortho, (mo_num)]
+
+  BEGIN_DOC
+  ! Overlap matrix between the RIGHT and LEFT MOs. Should be the identity matrix 
+  END_DOC
+
+  implicit none
+  integer                       :: i, k, m, n
+  double precision              :: accu_d, accu_nd 
+  double precision, allocatable :: tmp(:,:)
+  
+ ! TODO : re do the DEGEMM
+
+  overlap_bi_ortho = 0.d0
+  do i = 1, mo_num
+    do k = 1, mo_num
+      do m = 1, ao_num
+        do n = 1, ao_num
+          overlap_bi_ortho(k,i) += ao_overlap(n,m) * mo_l_coef(n,k) * mo_r_coef(m,i)
+        enddo
+      enddo
+    enddo
+  enddo
+
+!  allocate( tmp(mo_num,ao_num) )
+!
+!  ! tmp <-- L.T x S_ao
+!  call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0                         & 
+!            , mo_l_coef, size(mo_l_coef, 1), ao_overlap, size(ao_overlap, 1) &
+!            , 0.d0, tmp, size(tmp, 1) )
+!
+!  ! S <-- tmp x R
+!  call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0           & 
+!            , tmp, size(tmp, 1), mo_r_coef, size(mo_r_coef, 1) &
+!            , 0.d0, overlap_bi_ortho, size(overlap_bi_ortho, 1) )
+!
+!  deallocate( tmp )
+
+  do i = 1, mo_num
+    overlap_diag_bi_ortho(i) = overlap_bi_ortho(i,i)
+  enddo
+
+  accu_d  = 0.d0
+  accu_nd = 0.d0
+  do i = 1, mo_num
+    do k = 1, mo_num
+      if(i==k) then
+        accu_d += dabs(overlap_bi_ortho(k,i))
+      else
+        accu_nd += dabs(overlap_bi_ortho(k,i))
+      endif
+    enddo 
+  enddo
+  accu_d = accu_d/dble(mo_num)
+  accu_nd = accu_nd/dble(mo_num**2-mo_num)
+  if(dabs(accu_d-1.d0).gt.1.d-10.or.dabs(accu_nd).gt.1.d-10)then
+    print*,'Warning !!!'
+    print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0)
+    print*,'And bi orthogonality is off by an average of ',accu_nd
+    print*,'****************'
+    print*,'Overlap matrix betwee mo_l_coef and mo_r_coef  '
+    do i = 1, mo_num
+      write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
+    enddo
+  endif
+  print*,'Average trace of overlap_bi_ortho (should be 1.)'
+  print*,'accu_d  = ',accu_d
+  print*,'Sum of off diagonal terms of overlap_bi_ortho (should be zero)'
+  print*,'accu_nd = ',accu_nd
+  print*,'****************'
+ 
+END_PROVIDER 
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, overlap_mo_r, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, overlap_mo_l, (mo_num, mo_num)]
+
+  BEGIN_DOC
+  ! overlap_mo_r_mo(j,i) = <MO_i|MO_R_j>
+  END_DOC
+
+  implicit none
+  integer :: i, j, p, q
+
+  overlap_mo_r = 0.d0
+  overlap_mo_l = 0.d0
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do p = 1, ao_num
+        do q = 1, ao_num
+          overlap_mo_r(j,i) += mo_r_coef(q,i) * mo_r_coef(p,j) * ao_overlap(q,p) 
+          overlap_mo_l(j,i) += mo_l_coef(q,i) * mo_l_coef(p,j) * ao_overlap(q,p)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, overlap_mo_r_mo, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, overlap_mo_l_mo, (mo_num, mo_num)]
+
+  BEGIN_DOC
+  ! overlap_mo_r_mo(j,i) = <MO_j|MO_R_i>
+  END_DOC
+
+  implicit none
+  integer :: i, j, p, q
+
+  overlap_mo_r_mo = 0.d0
+  overlap_mo_l_mo = 0.d0
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do p = 1, ao_num
+        do q = 1, ao_num
+          overlap_mo_r_mo(j,i) += mo_coef(p,j) * mo_r_coef(q,i) * ao_overlap(q,p)
+          overlap_mo_l_mo(j,i) += mo_coef(p,j) * mo_l_coef(q,i) * ao_overlap(q,p)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, angle_left_right, (mo_num)]
+&BEGIN_PROVIDER [ double precision, max_angle_left_right]
+
+  BEGIN_DOC
+  ! angle_left_right(i) = angle between the left-eigenvector chi_i and the right-eigenvector phi_i
+  END_DOC
+
+  implicit none
+  integer          :: i, j
+  double precision :: left, right, arg
+  double precision :: angle(mo_num)
+
+  do i = 1, mo_num
+    left  = overlap_mo_l(i,i)
+    right = overlap_mo_r(i,i)
+    arg = min(overlap_bi_ortho(i,i)/(left*right),1.d0)
+    arg = max(arg, -1.d0)
+    angle_left_right(i) = dacos(arg) * 180.d0/dacos(-1.d0)
+  enddo
+
+  angle(1:mo_num) = dabs(angle_left_right(1:mo_num))
+  max_angle_left_right = maxval(angle)
+
+END_PROVIDER 
+
+! ---
+
+