diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f
index c06d64bb..93ac459e 100644
--- a/src/ao_many_one_e_ints/listj1b.irp.f
+++ b/src/ao_many_one_e_ints/listj1b.irp.f
@@ -336,7 +336,8 @@ END_PROVIDER
         r2 = dx*dx + dy*dy + dz*dz
         
         ii = ii + 1
-        List_all_comb_b3_coef(  ii) = dexp(-tmp1*tmp2*tmp4*r2)
+        ! x 2 to avoid doing integrals twice
+        List_all_comb_b3_coef(  ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2)
         List_all_comb_b3_expo(  ii) = tmp3
         List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj)
         List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj)
diff --git a/src/non_h_ints_mu/debug_fit.irp.f b/src/non_h_ints_mu/debug_fit.irp.f
index 547fd198..146028d5 100644
--- a/src/non_h_ints_mu/debug_fit.irp.f
+++ b/src/non_h_ints_mu/debug_fit.irp.f
@@ -266,6 +266,8 @@ subroutine test_list_b3()
     acc_ij = dabs(i_exc - i_num)
     if(acc_ij .gt. eps_ij) then
       print *, ' problem in list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
       print *, ' analyt = ', i_exc
       print *, ' numeri = ', i_num
       print *, ' diff   = ', acc_ij
diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f
index b4ddb606..3f1a9bf5 100644
--- a/src/non_h_ints_mu/grad_squared.irp.f
+++ b/src/non_h_ints_mu/grad_squared.irp.f
@@ -403,7 +403,9 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
     ! ---
 
-    if((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
+    if(((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) .and. use_ipp) then
+
+      print*, " going through Manu's IPP"
 
       ! an additional term is added here directly instead of 
       ! being added in int2_grad1_u12_square_ao for performance
diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f
index efa9be43..d5995ae5 100644
--- a/src/non_h_ints_mu/tc_integ.irp.f
+++ b/src/non_h_ints_mu/tc_integ.irp.f
@@ -219,28 +219,50 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
     !$OMP END DO
     !$OMP END PARALLEL
 
-  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4))  then
 
-    ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance
-    !PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
-    PROVIDE u12sq_j1bsq grad12_j12
+    if(use_ipp) then
 
-    int2_grad1_u12_square_ao = 0.d0
-    !$OMP PARALLEL               &
-    !$OMP DEFAULT (NONE)         &
-    !$OMP PRIVATE (i, j, ipoint) &
-    !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12)
-    !$OMP DO SCHEDULE (static)
-    do ipoint = 1, n_points_final_grid
-      do j = 1, ao_num
-        do i = 1, ao_num
-          !int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
-          int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+      ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance
+      PROVIDE u12sq_j1bsq grad12_j12
+
+      int2_grad1_u12_square_ao = 0.d0
+      !$OMP PARALLEL               &
+      !$OMP DEFAULT (NONE)         &
+      !$OMP PRIVATE (i, j, ipoint) &
+      !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12)
+      !$OMP DO SCHEDULE (static)
+      do ipoint = 1, n_points_final_grid
+        do j = 1, ao_num
+          do i = 1, ao_num
+            int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+          enddo
         enddo
       enddo
-    enddo
-    !$OMP END DO
-    !$OMP END PARALLEL
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    else
+
+      PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+
+      int2_grad1_u12_square_ao = 0.d0
+      !$OMP PARALLEL               &
+      !$OMP DEFAULT (NONE)         &
+      !$OMP PRIVATE (i, j, ipoint) &
+      !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12, u12_grad1_u12_j1b_grad1_j1b)
+      !$OMP DO SCHEDULE (static)
+      do ipoint = 1, n_points_final_grid
+        do j = 1, ao_num
+          do i = 1, ao_num
+            int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    endif
 
   elseif(j1b_type .ge. 100) then
 
diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg
index 41c07d0b..3a26a6eb 100644
--- a/src/tc_keywords/EZFIO.cfg
+++ b/src/tc_keywords/EZFIO.cfg
@@ -220,4 +220,9 @@ doc: If |true|, save the bi-ortho wave functions in a sorted way
 interface: ezfio,provider,ocaml
 default: True
 
+[use_ipp]
+type: logical
+doc: If |true|, use Manu IPP
+interface: ezfio,provider,ocaml
+default: True