Merge pull request #283 from AbdAmmar/dev-stable-tc-scf

Dev stable tc scf

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)

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

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

src/non_h_ints_mu/tc_integ.irp.f

@@ -221,8 +221,9 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
 
   elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4))  then
 
+    if(use_ipp) 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
 
       int2_grad1_u12_square_ao = 0.d0
@@ -234,7 +235,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
       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)
           enddo
         enddo
@@ -242,6 +242,28 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
       !$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
 
     PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra

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