Starting to transpose

src/ao_many_one_e_ints/ao_gaus_gauss.irp.f

@@ -150,6 +150,58 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
 
 end function overlap_gauss_r12_ao
 
+! --
+
+subroutine overlap_gauss_r12_ao_v(D_center, delta, i, j, resv, n_points)
+
+  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, n_points
+  double precision, intent(in) :: D_center(3,n_points), delta
+  double precision, intent(out) :: resv(n_points)
+
+  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(:)
+
+  double precision, external   :: overlap_gauss_r12
+  integer                        :: ipoint
+
+  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, delta, A_center, B_center, power_A, power_B, alpha, beta, analytical_j, n_points)
+      do ipoint=1, n_points
+        resv(ipoint) = resv(ipoint) + coef*analytical_j(ipoint)
+      enddo
+    enddo
+  enddo
+  deallocate(analytical_j)
+
+end
+
 ! ---
 
 double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, delta, i, j)
@@ -241,9 +293,6 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, delta, i, j,
   double precision             :: gama, gama_inv
   double precision             :: bg, dg, bdg
 
-  double precision, external   :: overlap_gauss_r12, overlap_gauss_r12_ao
-
-  double precision, external   :: overlap_gauss_r12_ao_with1s
   integer :: ipoint
 
   double precision, allocatable :: fact_g(:), G_center(:,:), analytical_j(:)
@@ -255,9 +304,7 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, delta, i, j,
   ASSERT(beta .gt. 0.d0)
 
   if(beta .lt. 1d-10) then
-    do ipoint=1,n_points
+    call overlap_gauss_r12_ao_v(D_center, delta, i, j, resv, n_points)
-      resv(ipoint) = overlap_gauss_r12_ao(D_center(1,ipoint), delta, i, j)
-    enddo
     return
   endif
 
@@ -317,6 +364,7 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, delta, i, j,
 
     enddo
   enddo
+  deallocate (fact_g, G_center, analytical_j )
 
 
 end

src/utils/one_e_integration.irp.f

@@ -151,4 +151,67 @@ subroutine overlap_x_abs(A_center, B_center, alpha, beta, power_A, power_B, over
 end
 
 
+! ---
 
+subroutine overlap_gaussian_xyz_v(A_center,B_center,alpha,beta,power_A,&
+      power_B,overlap_x,overlap_y,overlap_z,overlap,dim, n_points)
+  implicit none
+  BEGIN_DOC
+  !.. math::
+  !
+  !   S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2)  (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\
+  !   S = S_x S_y S_z
+  !
+  END_DOC
+  include 'constants.include.F'
+  integer,intent(in)             :: dim, n_points
+  double precision,intent(in)    :: A_center(3),B_center(3)  ! center of the x1 functions
+  double precision, intent(in)   :: alpha,beta
+  integer,intent(in)             :: power_A(3), power_B(3) ! power of the x1 functions
+  double precision, intent(out)  :: overlap_x(n_points),overlap_y(n_points),overlap_z(n_points),overlap(n_points)
+  double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,p
+  double precision               :: F_integral_tab(0:max_dim)
+  integer                        :: iorder_p(3)
+
+  call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
+   if(fact_p.lt.1d-20)then
+     overlap_x = 1.d-10
+     overlap_y = 1.d-10
+     overlap_z = 1.d-10
+     overlap = 1.d-10
+     return
+   endif
+  integer                        :: nmax
+  double precision               :: F_integral
+  nmax = maxval(iorder_p)
+  do i = 0,nmax
+    F_integral_tab(i) = F_integral(i,p)
+  enddo
+  overlap_x = P_new(0,1) * F_integral_tab(0)
+  overlap_y = P_new(0,2) * F_integral_tab(0)
+  overlap_z = P_new(0,3) * F_integral_tab(0)
+
+  integer                        :: i
+  do i = 1,iorder_p(1)
+    overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i)
+  enddo
+  call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1))
+  overlap_x *= fact_p
+
+  do i = 1,iorder_p(2)
+    overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i)
+  enddo
+  call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2))
+  overlap_y *= fact_p
+
+  do i = 1,iorder_p(3)
+    overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i)
+  enddo
+  call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3))
+  overlap_z *= fact_p
+
+  overlap = overlap_x * overlap_y * overlap_z
+
+end
+
+! ---