v0 of new Jast added

src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f

@@ -90,7 +90,7 @@ 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|
+  ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
   END_DOC
 
   implicit none
@@ -119,7 +119,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_num, n_points_final_grid)]
 
   BEGIN_DOC
-  ! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
+  ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
   END_DOC
 
   implicit none

src/bi_ort_ints/semi_num_ints_mo.irp.f

@@ -1,81 +1,154 @@
-BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu, ( mo_num, mo_num,n_points_final_grid)]
+
- implicit none
+! ---
- BEGIN_DOC
+
-! mo_v_ki_bi_ortho_erf_rk_cst_mu(k,i,ip) = int dr chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1 )/(2|r - R_ip|) on the BI-ORTHO MO basis 
+! TODO :: optimization : transform into a DGEMM
-! 
+
-! where phi_k(r) is a LEFT MOs and phi_i(r) is a RIGHT MO
+BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, n_points_final_grid)]
-!
+
-! R_ip = the "ip"-th point of the DFT Grid
+  BEGIN_DOC
- END_DOC
+  !
- integer :: ipoint
+  ! mo_v_ki_bi_ortho_erf_rk_cst_mu(k,i,ip) = int dr chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1 )/(2|r - R_ip|) on the BI-ORTHO MO basis 
- !$OMP PARALLEL                  &
+  ! 
- !$OMP DEFAULT (NONE)            &
+  ! where phi_k(r) is a LEFT MOs and phi_i(r) is a RIGHT MO
+  !
+  ! R_ip = the "ip"-th point of the DFT Grid
+  !
+  END_DOC
+
+  implicit none
+  integer :: ipoint
+ !$OMP PARALLEL         &
+ !$OMP DEFAULT (NONE)   &
  !$OMP PRIVATE (ipoint) & 
  !$OMP SHARED (n_points_final_grid,v_ij_erf_rk_cst_mu,mo_v_ki_bi_ortho_erf_rk_cst_mu)
  !$OMP DO SCHEDULE (dynamic)
-! TODO :: optimization : transform into a DGEMM
+  do ipoint = 1, n_points_final_grid
- do ipoint = 1, n_points_final_grid
+    call ao_to_mo_bi_ortho( v_ij_erf_rk_cst_mu            (1,1,ipoint), size(v_ij_erf_rk_cst_mu,             1) &
-   call ao_to_mo_bi_ortho(v_ij_erf_rk_cst_mu(1,1,ipoint),size(v_ij_erf_rk_cst_mu,1),mo_v_ki_bi_ortho_erf_rk_cst_mu(1,1,ipoint),size(mo_v_ki_bi_ortho_erf_rk_cst_mu,1))
+                          , mo_v_ki_bi_ortho_erf_rk_cst_mu(1,1,ipoint), size(mo_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
- enddo
+  enddo
  !$OMP END DO
  !$OMP END PARALLEL
- mo_v_ki_bi_ortho_erf_rk_cst_mu = mo_v_ki_bi_ortho_erf_rk_cst_mu * 0.5d0
-END_PROVIDER 
 
-
+  mo_v_ki_bi_ortho_erf_rk_cst_mu = mo_v_ki_bi_ortho_erf_rk_cst_mu * 0.5d0
-BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu_transp, ( n_points_final_grid,mo_num, mo_num)]
- implicit none
- BEGIN_DOC
-! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the BI-ORTHO MO basis
- END_DOC
- integer :: ipoint,i,j
- do i = 1, mo_num
-  do j = 1, mo_num
-   do ipoint = 1, n_points_final_grid
-    mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ki_bi_ortho_erf_rk_cst_mu(j,i,ipoint)
-   enddo
-  enddo
- enddo
-! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu
-END_PROVIDER 
-
-BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, ( mo_num, mo_num,3,n_points_final_grid)]
- implicit none
- BEGIN_DOC
-! mo_x_v_ki_bi_ortho_erf_rk_cst_mu(k,i,m,ip) = int dr x(m) * chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1)/2|r - R_ip| on the BI-ORTHO MO basis 
-!
-! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => x(m) = x, m=2 => x(m) = y, m=3 => x(m) = z,
-!
-! R_ip = the "ip"-th point of the DFT Grid
- END_DOC
- integer :: ipoint,m
- !$OMP PARALLEL                  &
- !$OMP DEFAULT (NONE)            &
- !$OMP PRIVATE (ipoint,m) & 
- !$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ki_bi_ortho_erf_rk_cst_mu)
- !$OMP DO SCHEDULE (dynamic)
-! TODO :: optimization : transform into a DGEMM
- do ipoint = 1, n_points_final_grid
-  do m = 1, 3
-   call ao_to_mo_bi_ortho(x_v_ij_erf_rk_cst_mu_transp(1,1,m,ipoint),size(x_v_ij_erf_rk_cst_mu_transp,1),mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,m,ipoint),size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu,1))
-  enddo
- enddo
- !$OMP END DO
- !$OMP END PARALLEL
- mo_x_v_ki_bi_ortho_erf_rk_cst_mu = 0.5d0 * mo_x_v_ki_bi_ortho_erf_rk_cst_mu
 
 END_PROVIDER 
 
 ! ---
-BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, 3, mo_num, mo_num)]
+
+BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, mo_num, mo_num)]
+
+  BEGIN_DOC
+  !
+  ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the BI-ORTHO MO basis
+  !
+  END_DOC
+
   implicit none
-  integer :: i, j, m, ipoint
+  integer :: ipoint, i, j
+
   do i = 1, mo_num
     do j = 1, mo_num
-      do m = 1, 3
+      do ipoint = 1, n_points_final_grid
-        do ipoint = 1, n_points_final_grid
+        mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ki_bi_ortho_erf_rk_cst_mu(j,i,ipoint)
-          mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,m,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,m,ipoint)
+      enddo
-        enddo
+    enddo
+  enddo
+
+! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu
+
+END_PROVIDER 
+
+! ---
+
+! TODO :: optimization : transform into a DGEMM
+
+BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, 3, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! mo_x_v_ki_bi_ortho_erf_rk_cst_mu(k,i,m,ip) = int dr x(m) * chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1)/2|r - R_ip| on the BI-ORTHO MO basis 
+  !
+  ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => x(m) = x, m=2 => x(m) = y, m=3 => x(m) = z,
+  !
+  ! R_ip = the "ip"-th point of the DFT Grid
+  !
+  END_DOC
+
+  implicit none
+  integer :: ipoint
+
+ !$OMP PARALLEL         &
+ !$OMP DEFAULT (NONE)   &
+ !$OMP PRIVATE (ipoint) & 
+ !$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ki_bi_ortho_erf_rk_cst_mu)
+ !$OMP DO SCHEDULE (dynamic)
+  do ipoint = 1, n_points_final_grid
+
+    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,1,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
+                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,1,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+
+    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,2,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
+                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,2,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+
+    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,3,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
+                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,3,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  mo_x_v_ki_bi_ortho_erf_rk_cst_mu = 0.5d0 * mo_x_v_ki_bi_ortho_erf_rk_cst_mu
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int2_grad1_u12_bimo(:,k,i,ipoint) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \chi_k(r2) \phi_i(r2) 
+  !
+  END_DOC
+
+  implicit none
+  integer :: ipoint
+
+ !$OMP PARALLEL         &
+ !$OMP DEFAULT (NONE)   &
+ !$OMP PRIVATE (ipoint) & 
+ !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao,int2_grad1_u12_bimo)
+ !$OMP DO SCHEDULE (dynamic)
+  do ipoint = 1, n_points_final_grid
+
+    call ao_to_mo_bi_ortho( int2_grad1_u12_ao  (1,1,1,ipoint), size(int2_grad1_u12_ao  , 1) &
+                          , int2_grad1_u12_bimo(1,1,1,ipoint), size(int2_grad1_u12_bimo, 1) )
+
+    call ao_to_mo_bi_ortho( int2_grad1_u12_ao  (2,1,1,ipoint), size(int2_grad1_u12_ao  , 1) &
+                          , int2_grad1_u12_bimo(2,1,1,ipoint), size(int2_grad1_u12_bimo, 1) )
+
+    call ao_to_mo_bi_ortho( int2_grad1_u12_ao  (3,1,1,ipoint), size(int2_grad1_u12_ao  , 1) &
+                          , int2_grad1_u12_bimo(3,1,1,ipoint), size(int2_grad1_u12_bimo, 1) )
+
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, 3, mo_num, mo_num)]
+
+  implicit none
+  integer :: i, j, ipoint
+
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do ipoint = 1, n_points_final_grid
+        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,1,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,1,ipoint)
+        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,2,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,2,ipoint)
+        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,3,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,3,ipoint)
       enddo
     enddo
   enddo
@@ -83,14 +156,15 @@ END_PROVIDER
 
 ! ---
 
-
 BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, 3, mo_num, mo_num)]
+
   BEGIN_DOC
-  ! x_W_ki_bi_ortho_erf_rk(ip,m,k,i) = \int dr chi_k(r) (1 - erf(mu |r-R_ip|)) (x(m)-X(m)_ip) phi_i(r) ON THE BI-ORTHO MO BASIS 
+  !
-!
+  ! x_W_ki_bi_ortho_erf_rk(ip,m,k,i) = \int dr chi_k(r) \frac{(1 - erf(mu |r-R_ip|))}{2|r-R_ip|} (x(m)-R_ip(m)) phi_i(r) ON THE BI-ORTHO MO BASIS 
-! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z,
+  !
-!
+  ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z,
-! R_ip = the "ip"-th point of the DFT Grid
+  !
+  ! R_ip = the "ip"-th point of the DFT Grid
   END_DOC
  
   implicit none
@@ -100,7 +174,7 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid,
   double precision :: xyz
   double precision :: wall0, wall1
  
-  print*,'providing x_W_ki_bi_ortho_erf_rk ...'
+  print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
   call wall_time(wall0)
 
  !$OMP PARALLEL                   &
@@ -126,7 +200,7 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid,
  ! FREE mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp
 
   call wall_time(wall1)
-  print*,'time to provide x_W_ki_bi_ortho_erf_rk = ',wall1 - wall0
+  print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
 
 END_PROVIDER 
 

src/bi_ort_ints/three_body_ints_bi_ort.irp.f

@@ -53,26 +53,55 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n
 
 END_PROVIDER 
 
-subroutine give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral)
+! ---
- implicit none
+
- double precision, intent(out) :: integral
+subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
- integer, intent(in) :: n,l,k,m,j,i
+
- double precision :: weight
+  BEGIN_DOC
- BEGIN_DOC
+  !
-! <n l k|-L|m j i> with a BI ORTHONORMAL ORBITALS 
+  ! < n l k | -L | m j i > with a BI-ORTHONORMAL ORBITALS 
- END_DOC
+  !
- integer :: ipoint,mm
+  END_DOC
- integral = 0.d0
+
- do mm = 1, 3
+  implicit none
-  do ipoint = 1, n_points_final_grid
+  integer,           intent(in) :: n,l,k,m,j,i
-   weight = final_weight_at_r_vector(ipoint)                                                                          
+  double precision, intent(out) :: integral
-   integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & 
+  integer                       :: ipoint
-                      * x_W_ki_bi_ortho_erf_rk(ipoint,mm,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,mm,l,j) 
+  double precision              :: weight
-   integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & 
+
-                      * x_W_ki_bi_ortho_erf_rk(ipoint,mm,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,mm,k,i) 
+  integral = 0.d0
-   integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & 
+  do ipoint = 1, n_points_final_grid
-                      * x_W_ki_bi_ortho_erf_rk(ipoint,mm,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,mm,k,i) 
+    weight = final_weight_at_r_vector(ipoint)                                                                          
-  enddo
+
- enddo
+!    integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & 
-end
+!              * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,1,l,j)  &
+!                + x_W_ki_bi_ortho_erf_rk(ipoint,2,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j)  &
+!                + x_W_ki_bi_ortho_erf_rk(ipoint,3,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j)  )
+!    integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & 
+!              * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,1,k,i)  &
+!                + x_W_ki_bi_ortho_erf_rk(ipoint,2,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i)  &
+!                + x_W_ki_bi_ortho_erf_rk(ipoint,3,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i)  )
+!    integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
+!              * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,1,k,i)  &
+!                + x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i)  &
+!                + x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i)  )
+
+    integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & 
+              * ( int2_grad1_u12_bimo(1,ipoint,n,m) * int2_grad1_u12_bimo(1,ipoint,l,j)  &
+                + int2_grad1_u12_bimo(2,ipoint,n,m) * int2_grad1_u12_bimo(2,ipoint,l,j)  &
+                + int2_grad1_u12_bimo(3,ipoint,n,m) * int2_grad1_u12_bimo(3,ipoint,l,j)  )
+    integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & 
+              * ( int2_grad1_u12_bimo(1,ipoint,n,m) * int2_grad1_u12_bimo(1,ipoint,k,i)  &
+                + int2_grad1_u12_bimo(2,ipoint,n,m) * int2_grad1_u12_bimo(2,ipoint,k,i)  &
+                + int2_grad1_u12_bimo(3,ipoint,n,m) * int2_grad1_u12_bimo(3,ipoint,k,i)  )
+    integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
+              * ( int2_grad1_u12_bimo(1,ipoint,l,j) * int2_grad1_u12_bimo(1,ipoint,k,i)  &
+                + int2_grad1_u12_bimo(2,ipoint,l,j) * int2_grad1_u12_bimo(2,ipoint,k,i)  &
+                + int2_grad1_u12_bimo(3,ipoint,l,j) * int2_grad1_u12_bimo(3,ipoint,k,i)  )
+
+  enddo
+
+end subroutine give_integrals_3_body_bi_ort
+
+! ---
 

src/non_h_ints_mu/debug_integ_jmu_modif.irp.f

@@ -34,7 +34,7 @@ program debug_integ_jmu_modif
   !call test_int2_u2_j1b2()
   !call test_int2_grad1u2_grad2u2_j1b2()
 
-  !call test_grad_1_u_ij_mu()
+  !call test_int2_grad1_u12_ao()
   !call test_gradu_squared_u_ij_mu()
 
 end
@@ -287,16 +287,16 @@ end subroutine test_int2_grad1u2_grad2u2_j1b2
 
 ! ---
 
-subroutine test_grad_1_u_ij_mu()
+subroutine test_int2_grad1_u12_ao()
 
   implicit none
   integer          :: i, j, ipoint
   double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
   double precision :: integ(3)
 
-  print*, ' test_grad_1_u_ij_mu ...'
+  print*, ' test_int2_grad1_u12_ao ...'
 
-  PROVIDE grad_1_u_ij_mu 
+  PROVIDE int2_grad1_u12_ao
 
   eps_ij  = 1d-6
   acc_tot = 0.d0
@@ -305,13 +305,13 @@ subroutine test_grad_1_u_ij_mu()
     do j = 1, ao_num
       do i = 1, ao_num
 
-        call num_grad_1_u_ij_mu(i, j, ipoint, integ)
+        call num_int2_grad1_u12_ao(i, j, ipoint, integ)
 
-        i_exc  = grad_1_u_ij_mu(i,j,ipoint,1) 
+        i_exc  = int2_grad1_u12_ao(1,i,j,ipoint) 
         i_num  = integ(1)
         acc_ij = dabs(i_exc - i_num)
         if(acc_ij .gt. eps_ij) then
-          print *, ' problem in x part of grad_1_u_ij_mu on', i, j, ipoint
+          print *, ' problem in x part of int2_grad1_u12_ao on', i, j, ipoint
           print *, ' analyt integ = ', i_exc
           print *, ' numeri integ = ', i_num
           print *, ' diff         = ', acc_ij
@@ -319,11 +319,11 @@ subroutine test_grad_1_u_ij_mu()
         acc_tot += acc_ij
         normalz += dabs(i_num)
 
-        i_exc  = grad_1_u_ij_mu(i,j,ipoint,2) 
+        i_exc  = int2_grad1_u12_ao(2,i,j,ipoint) 
         i_num  = integ(2)
         acc_ij = dabs(i_exc - i_num)
         if(acc_ij .gt. eps_ij) then
-          print *, ' problem in y part of grad_1_u_ij_mu on', i, j, ipoint
+          print *, ' problem in y part of int2_grad1_u12_ao on', i, j, ipoint
           print *, ' analyt integ = ', i_exc
           print *, ' numeri integ = ', i_num
           print *, ' diff         = ', acc_ij
@@ -331,11 +331,11 @@ subroutine test_grad_1_u_ij_mu()
         acc_tot += acc_ij
         normalz += dabs(i_num)
 
-        i_exc  = grad_1_u_ij_mu(i,j,ipoint,3) 
+        i_exc  = int2_grad1_u12_ao(3,i,j,ipoint) 
         i_num  = integ(3)
         acc_ij = dabs(i_exc - i_num)
         if(acc_ij .gt. eps_ij) then
-          print *, ' problem in z part of grad_1_u_ij_mu on', i, j, ipoint
+          print *, ' problem in z part of int2_grad1_u12_ao on', i, j, ipoint
           print *, ' analyt integ = ', i_exc
           print *, ' numeri integ = ', i_num
           print *, ' diff         = ', acc_ij
@@ -352,7 +352,7 @@ subroutine test_grad_1_u_ij_mu()
   print*, ' normalz        = ', normalz
 
   return
-end subroutine test_grad_1_u_ij_mu
+end subroutine test_int2_grad1_u12_ao
 
 ! ---
 

src/non_h_ints_mu/grad_squared.irp.f

@@ -4,20 +4,34 @@
 !       TODO : strong optmization : write the loops in a different way
 !            : for each couple of AO, the gaussian product are done once for all 
 
-BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_points_final_grid)]
+BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_points_final_grid) ]
 
   BEGIN_DOC
   !
-  ! -1/2 [ (grad_1 u)^2 + (grad_2 u^2)] = - 1/4 * (1 - erf(mu*r12))^2 
+  ! if J(r1,r2) = u12:
   !
+  ! gradu_squared_u_ij_mu = -0.50 x \int r2 [ (grad_1 u12)^2 + (grad_2 u12^2)] \phi_i(2) \phi_j(2)
+  !                       = -0.25 x \int r2       (1 - erf(mu*r12))^2          \phi_i(2) \phi_j(2)
   ! and
   !   (1 - erf(mu*r12))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * r12^2)
   !
+  ! if J(r1,r2) = u12 x v1 x v2
+  !
+  ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2  v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2  v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ]
+  !                       = -0.25 x        v1^2      \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2
+  !                       + -0.50 x    (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2)      u12^2          v2^2
+  !                       + -1.00 x v1 (grad_1 v1)   \int r2 \phi_i(2) \phi_j(2)   (grad_1 u12)      v2^2 
+  !                       =                 v1^2        x   int2_grad1u2_grad2u2_j1b2
+  !                       + -0.5 x      (grad_1 v1)^2   x   int2_u2_j1b2
+  !                       + -1.0 X V1 x (grad_1 v1)   \cdot int2_u_grad1u_x_j1b
+  !
+  !
   END_DOC
  
   implicit none
   integer                    :: ipoint, i, j, m, igauss
-  double precision           :: r(3), delta, coef, tmp
+  double precision           :: r(3), delta, coef
+  double precision           :: tmp_v, tmp_x, tmp_y, tmp_z, tmp1, tmp2, tmp3, tmp4, tmp5
   double precision           :: time0, time1
   double precision, external :: overlap_gauss_r12_ao
 
@@ -27,13 +41,28 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_poin
   PROVIDE j1b_type j1b_pen
 
   if(j1b_type .eq. 3) then
-    ! v1_1b^2 \int d2 \phi_i(2) \phi_j(2) \frac{-[1 - \erf(\mu r12)]^2}{4} v2_1b^2
 
     do ipoint = 1, n_points_final_grid
-      tmp = v_1b(ipoint) * v_1b(ipoint) 
+
+      tmp_v = v_1b       (ipoint)
+      tmp_x = v_1b_grad(1,ipoint)
+      tmp_y = v_1b_grad(2,ipoint)
+      tmp_z = v_1b_grad(3,ipoint)
+
+      tmp1 = tmp_v * tmp_v
+      tmp2 = 0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
+      tmp3 = tmp_v * tmp_x
+      tmp4 = tmp_v * tmp_y
+      tmp5 = tmp_v * tmp_z
+
       do j = 1, ao_num
         do i = 1, ao_num
-          gradu_squared_u_ij_mu(j,i,ipoint) += tmp * int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
+
+          gradu_squared_u_ij_mu(j,i,ipoint) += tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) &
+                                             - tmp2 * int2_u2_j1b2             (i,j,ipoint) &
+                                             - tmp3 * int2_u_grad1u_x_j1b    (1,i,j,ipoint) &
+                                             - tmp4 * int2_u_grad1u_x_j1b    (2,i,j,ipoint) &
+                                             - tmp5 * int2_u_grad1u_x_j1b    (3,i,j,ipoint) 
         enddo
       enddo
     enddo
@@ -74,22 +103,27 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   implicit none
   integer                       :: ipoint, i, j, k, l
-  double precision              :: contrib, weight1
+  double precision              :: contrib, weight1, ao_k_r, ao_i_r
   double precision, allocatable :: ac_mat(:,:,:,:)
 
   allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
   ac_mat = 0.d0
 
   do ipoint = 1, n_points_final_grid
-    weight1 = final_weight_at_r_vector(ipoint)
+    weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
 
-    do j = 1, ao_num
+    do i = 1, ao_num
-      do l = 1, ao_num
+      ao_i_r = aos_in_r_array_transp(ipoint,i)
-        do i = 1, ao_num
+
-          do k = 1, ao_num
+      do k = 1, ao_num
-            contrib = weight1 * 0.5d0 * (aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i))
+        ao_k_r = aos_in_r_array_transp(ipoint,k)
-            ! \int dr1 phi_k(r1) phi_i(r1) . \int dr2 |\grad_1 u(r1,r2)|^2 \phi_l(r2) \phi_j(r2)
+
-            ac_mat(k,i,l,j) += gradu_squared_u_ij_mu(l,j,ipoint) * contrib
+        do j = 1, ao_num
+          do l = 1, ao_num
+
+            contrib = gradu_squared_u_ij_mu(l,j,ipoint) * ao_k_r * ao_i_r
+
+            ac_mat(k,i,l,j) += weight1 * contrib
           enddo
         enddo
       enddo

src/non_h_ints_mu/new_grad_tc.irp.f

@@ -1,15 +1,28 @@
   
 ! ---
 
-BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)]
+BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_points_final_grid)]
 
   BEGIN_DOC
   !
-  ! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 [-1 * \grad_r1 u(r1,r2)]              \phi_i(r2) \phi_j(r2) x 1s_j1b(r2)
+  ! int2_grad1_u12_ao(:,i,j,ipoint) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) 
-  !                            = \int dr2 [(r1 - r2) (erf(mu * r12)-1)/2 r_12]  \phi_i(r2) \phi_j(r2) x 1s_j1b(r2)
   !
   ! where r1 = r(ipoint)
   !
+  ! if J(r1,r2) = u12:
+  !
+  ! int2_grad1_u12_ao(:,i,j,ipoint) = 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r2) \phi_j(r2)
+  !                                 = 0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ]
+  !
+  ! if J(r1,r2) = u12 x v1 x v2
+  !
+  ! int2_grad1_u12_ao(:,i,j,ipoint) =      v1    x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
+  !                                 + \grad_1 v1 x [       \int dr2                  u12              v2 \phi_i(r2) \phi_j(r2) ] 
+  !                                 =    0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) 
+  !                                 -    0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) 
+  !                                 + v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
+  !
+  !
   END_DOC
 
   implicit none
@@ -25,10 +38,10 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin
       y = final_grid_points(2,ipoint)
       z = final_grid_points(3,ipoint)
 
-      tmp0  = v_1b       (ipoint)
+      tmp0  = 0.5d0 * v_1b(ipoint)
-      tmp_x = v_1b_grad(1,ipoint)
+      tmp_x =  v_1b_grad(1,ipoint)
-      tmp_y = v_1b_grad(2,ipoint)
+      tmp_y =  v_1b_grad(2,ipoint)
-      tmp_z = v_1b_grad(3,ipoint)
+      tmp_z =  v_1b_grad(3,ipoint)
   
       do j = 1, ao_num
         do i = 1, ao_num
@@ -36,9 +49,9 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin
           tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) 
           tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) 
 
-          grad_1_u_ij_mu(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + tmp_x * tmp2
+          int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + tmp_x * tmp2
-          grad_1_u_ij_mu(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + tmp_y * tmp2
+          int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + tmp_y * tmp2
-          grad_1_u_ij_mu(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) + tmp_z * tmp2
+          int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) + tmp_z * tmp2
         enddo
       enddo
     enddo
@@ -51,61 +64,87 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin
       z = final_grid_points(3,ipoint)
       do j = 1, ao_num
         do i = 1, ao_num
-          grad_1_u_ij_mu(i,j,ipoint,1) = v_ij_erf_rk_cst_mu(i,j,ipoint) * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)
+          int2_grad1_u12_ao(1,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)
-          grad_1_u_ij_mu(i,j,ipoint,2) = v_ij_erf_rk_cst_mu(i,j,ipoint) * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)
+          int2_grad1_u12_ao(2,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)
-          grad_1_u_ij_mu(i,j,ipoint,3) = v_ij_erf_rk_cst_mu(i,j,ipoint) * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)
+          int2_grad1_u12_ao(3,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)
         enddo
       enddo
     enddo
 
-  endif
+    int2_grad1_u12_ao *= 0.5d0
 
-  grad_1_u_ij_mu *= 0.5d0
+  endif
 
 END_PROVIDER 
 
 ! ---
 
 BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)]
- implicit none
- BEGIN_DOC
- ! tc_grad_and_lapl_ao(k,i,l,j) = <kl | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij>
- !
- ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) 
- !
- ! This is obtained by integration by parts. 
- END_DOC
- integer :: ipoint,i,j,k,l,m
- double precision :: contrib,weight1
- double precision, allocatable :: ac_mat(:,:,:,:)
- allocate(ac_mat(ao_num, ao_num, ao_num, ao_num))
- ac_mat = 0.d0
- do m = 1, 3
-  do ipoint = 1, n_points_final_grid
-   weight1 = final_weight_at_r_vector(ipoint)
-   do j = 1, ao_num
-    do l = 1, ao_num
-     do i = 1, ao_num
-      do k = 1, ao_num
-      contrib = weight1 *0.5D0* (aos_in_r_array_transp(ipoint,k) * aos_grad_in_r_array_transp_bis(ipoint,i,m)  &
-                                -aos_in_r_array_transp(ipoint,i) * aos_grad_in_r_array_transp_bis(ipoint,k,m) )
-      ! \int dr1 phi_k(r1) \grad_r1 phi_i(r1) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2)
-       ac_mat(k,i,l,j) += grad_1_u_ij_mu(l,j,ipoint,m) * contrib
-      enddo
-     enddo
-    enddo
-   enddo
-  enddo
- enddo
 
- do j = 1, ao_num
+  BEGIN_DOC
-  do l = 1, ao_num
+  !
-   do i = 1, ao_num
+  ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij >
-    do k = 1, ao_num
+  !
-     tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
+  ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) 
+  !
+  ! This is obtained by integration by parts. 
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: ipoint, i, j, k, l
+  double precision              :: contrib, weight1, contrib_x, contrib_y, contrib_z
+  double precision              :: ao_k_r, ao_k_dx, ao_k_dy, ao_k_dz
+  double precision              :: ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz
+  double precision, allocatable :: ac_mat(:,:,:,:)
+
+  allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
+  ac_mat = 0.d0
+
+  do ipoint = 1, n_points_final_grid
+    weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
+
+    do i = 1, ao_num
+      ao_i_r  = aos_in_r_array_transp         (ipoint,i)
+      ao_i_dx = aos_grad_in_r_array_transp_bis(ipoint,i,1)
+      ao_i_dy = aos_grad_in_r_array_transp_bis(ipoint,i,2)
+      ao_i_dz = aos_grad_in_r_array_transp_bis(ipoint,i,3)
+
+      do k = 1, ao_num
+        ao_k_r  = aos_in_r_array_transp         (ipoint,k)
+        ao_k_dx = aos_grad_in_r_array_transp_bis(ipoint,k,1)
+        ao_k_dy = aos_grad_in_r_array_transp_bis(ipoint,k,2)
+        ao_k_dz = aos_grad_in_r_array_transp_bis(ipoint,k,3)
+
+        do j = 1, ao_num
+          do l = 1, ao_num
+
+            contrib_x = int2_grad1_u12_ao(1,l,j,ipoint) * ( ao_k_r * ao_i_dx - ao_i_r * ao_k_dx )
+            contrib_y = int2_grad1_u12_ao(2,l,j,ipoint) * ( ao_k_r * ao_i_dy - ao_i_r * ao_k_dy )
+            contrib_z = int2_grad1_u12_ao(3,l,j,ipoint) * ( ao_k_r * ao_i_dz - ao_i_r * ao_k_dz )
+
+            contrib   = weight1 * ( contrib_x + contrib_y + contrib_z )
+
+            ac_mat(k,i,l,j) += contrib
+          enddo
+        enddo
+      enddo
     enddo
-   enddo
   enddo
- enddo
+ 
+  do j = 1, ao_num
+    do l = 1, ao_num
+      do i = 1, ao_num
+        do k = 1, ao_num
+          tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  deallocate(ac_mat)
 
 END_PROVIDER 
+
+! ---
+

src/non_h_ints_mu/numerical_integ.irp.f

@@ -1,51 +1,6 @@
 
 ! --- 
 
-!
-!            \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2)
-!
-
-!BEGIN_PROVIDER [ double precision, num_grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)]
-!
-!  implicit none
-!  
-!  integer                    :: i, j, ipoint, jpoint
-!  double precision           :: tmp, r1(3), r2(3), grad(3)
-!
-!  double precision, external :: ao_value
-!  double precision, external :: j12_nucl
-!
-!  num_grad_1_u_ij_mu = 0.d0
-!
-!  do j = 1, ao_num
-!    do i = 1, ao_num
-!
-!      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)
-!
-!        do jpoint = 1, n_points_final_grid
-!          r2(1) = final_grid_points(1,jpoint)
-!          r2(2) = final_grid_points(2,jpoint)
-!          r2(3) = final_grid_points(3,jpoint)
-!          tmp   = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
-!
-!          call grad1_j12_mu_exc(r1, r2, grad)
-!
-!          num_grad_1_u_ij_mu(i,j,ipoint,1) += tmp * (-1.d0 * grad(1))
-!          num_grad_1_u_ij_mu(i,j,ipoint,2) += tmp * (-1.d0 * grad(2))
-!          num_grad_1_u_ij_mu(i,j,ipoint,3) += tmp * (-1.d0 * grad(3))
-!        enddo
-!
-!      enddo
-!    enddo
-!  enddo
-!
-!END_PROVIDER
-
-! ---
-
 double precision function num_v_ij_u_cst_mu_j1b(i, j, ipoint)
 
   BEGIN_DOC
@@ -289,7 +244,7 @@ end subroutine num_x_v_ij_erf_rk_cst_mu_j1b
 
 ! ---
 
-subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ)
+subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ)
 
   implicit none
 
@@ -328,7 +283,7 @@ subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ)
   integ(3) = tmp_z
 
   return
-end subroutine num_grad_1_u_ij_mu
+end subroutine num_int2_grad1_u12_ao
 
 ! ---