OPTIM E_3e for TCSCF

src/bi_ort_ints/no_dressing.irp.f

@@ -1,7 +1,7 @@
 
 ! ---
 
-BEGIN_PROVIDER [double precision, noL_0e]
+BEGIN_PROVIDER [double precision, noL_0e_v0]
 
   implicit none
   integer                       :: i, j, k
@@ -10,7 +10,7 @@ BEGIN_PROVIDER [double precision, noL_0e]
   double precision, allocatable :: tmp(:)
 
   call wall_time(t0)
-  print*, " Providing noL_0e ..."
+  print*, " Providing noL_0e_v0 ..."
 
   if(elec_alpha_num .eq. elec_beta_num) then
 
@@ -40,7 +40,7 @@ BEGIN_PROVIDER [double precision, noL_0e]
     !$OMP END DO
     !$OMP END PARALLEL
 
-    noL_0e = -1.d0 * (sum(tmp)) / 6.d0
+    noL_0e_v0 = -1.d0 * (sum(tmp)) / 6.d0
 
     deallocate(tmp)
 
@@ -94,9 +94,8 @@ BEGIN_PROVIDER [double precision, noL_0e]
           call give_integrals_3_body_bi_ort(i, j, k, j, k, i, I_ijk_jki)
           call give_integrals_3_body_bi_ort(i, j, k, i, k, j, I_ijk_ikj)
           call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik)
-          call give_integrals_3_body_bi_ort(i, j, k, k, j, i, I_ijk_kji)
 
-          tmp(i) = tmp(i) + 6.d0 * (2.d0 * I_ijk_ijk + I_ijk_jki - I_ijk_ikj - I_ijk_jik - I_ijk_kji)
+          tmp(i) = tmp(i) + 6.d0 * (2.d0 * I_ijk_ijk + I_ijk_jki - I_ijk_ikj - 2.d0 * I_ijk_jik)
         enddo ! k
 
         do k = elec_beta_num+1, elec_alpha_num
@@ -104,26 +103,25 @@ BEGIN_PROVIDER [double precision, noL_0e]
           call give_integrals_3_body_bi_ort(i, j, k, i, j, k, I_ijk_ijk)
           call give_integrals_3_body_bi_ort(i, j, k, j, k, i, I_ijk_jki)
           call give_integrals_3_body_bi_ort(i, j, k, i, k, j, I_ijk_ikj)
-          call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik)
           call give_integrals_3_body_bi_ort(i, j, k, k, j, i, I_ijk_kji)
 
-          tmp(i) = tmp(i) + 3.d0 * (2.d0 * I_ijk_ijk + 2.d0 * I_ijk_jki - I_ijk_ikj - I_ijk_jik - 2.d0 * I_ijk_kji)
+          tmp(i) = tmp(i) + 6.d0 * (I_ijk_ijk + I_ijk_jki - I_ijk_ikj - I_ijk_kji)
         enddo ! k
       enddo ! j
     enddo ! i
     !$OMP END DO
     !$OMP END PARALLEL
 
-    noL_0e = -1.d0 * (sum(tmp)) / 6.d0
+    noL_0e_v0 = -1.d0 * (sum(tmp)) / 6.d0
 
     deallocate(tmp)
 
   endif
 
   call wall_time(t1)
-  print*, " Wall time for noL_0e (min) = ", (t1 - t0)/60.d0
+  print*, " Wall time for noL_0e_v0 (min) = ", (t1 - t0)/60.d0
 
-  print*, " noL_0e = ", noL_0e
+  print*, " noL_0e_v0 = ", noL_0e_v0
 
 END_PROVIDER
 
@@ -322,6 +320,403 @@ END_PROVIDER
 
 ! ---
 
+BEGIN_PROVIDER [double precision, noL_0e]
+
+  implicit none
+  integer                       :: i, j, k, ipoint
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp(:)
+  double precision, allocatable :: tmp_L(:,:), tmp_R(:,:)
+  double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:)
+  double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:)
+
+
+  call wall_time(t0)
+  print*, " Providing noL_0e ..."
+
+  if(elec_alpha_num .eq. elec_beta_num) then
+
+    allocate(tmp(elec_beta_num))
+    allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
+
+    !$OMP PARALLEL                                                 &
+    !$OMP DEFAULT(NONE)                                            &
+    !$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R)                      &
+    !$OMP SHARED(elec_beta_num, n_points_final_grid,               & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
+    !$OMP        int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
+
+    !$OMP DO
+    do j = 1, elec_beta_num
+
+      tmp_L = 0.d0
+      tmp_R = 0.d0
+      do i = 1, elec_beta_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+
+          tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+        enddo
+      enddo
+
+      tmp(j) = 0.d0
+      do ipoint = 1, n_points_final_grid
+        tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+      enddo
+    enddo ! j
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+    noL_0e = -2.d0 * sum(tmp)
+
+    deallocate(tmp)
+    deallocate(tmp_L, tmp_R)
+
+    ! ---
+
+    allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
+    tmp_O = 0.d0
+    tmp_J = 0.d0
+
+    !$OMP PARALLEL                                                  &
+    !$OMP DEFAULT(NONE)                                             &
+    !$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv)                &
+    !$OMP SHARED(elec_beta_num, n_points_final_grid,                & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,  &
+    !$OMP        int2_grad1_u12_bimo_t, tmp_O, tmp_J)
+
+    allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
+    tmp_O_priv = 0.d0
+    tmp_J_priv = 0.d0
+  
+    !$OMP DO 
+    do i = 1, elec_beta_num
+      do ipoint = 1, n_points_final_grid
+        tmp_O_priv(ipoint)   = tmp_O_priv(ipoint)   + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+        tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
+        tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
+        tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP CRITICAL
+    tmp_O = tmp_O + tmp_O_priv
+    tmp_J = tmp_J + tmp_J_priv
+    !$OMP END CRITICAL
+
+    deallocate(tmp_O_priv, tmp_J_priv)
+    !$OMP END PARALLEL
+
+    allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
+    tmp_M = 0.d0
+    tmp_S = 0.d0
+
+    !$OMP PARALLEL                                                 &
+    !$OMP DEFAULT(NONE)                                            &
+    !$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv)            &
+    !$OMP SHARED(elec_beta_num, n_points_final_grid,               & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
+    !$OMP        int2_grad1_u12_bimo_t, tmp_M, tmp_S)
+
+    allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
+    tmp_M_priv = 0.d0
+    tmp_S_priv = 0.d0
+  
+    !$OMP DO COLLAPSE(2)
+    do i = 1, elec_beta_num
+      do j = 1, elec_beta_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+
+          tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                      + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                      + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP CRITICAL
+    tmp_M = tmp_M + tmp_M_priv
+    tmp_S = tmp_S + tmp_S_priv
+    !$OMP END CRITICAL
+
+    deallocate(tmp_M_priv, tmp_S_priv)
+    !$OMP END PARALLEL
+
+    allocate(tmp(n_points_final_grid))
+
+    do ipoint = 1, n_points_final_grid
+
+      tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
+
+      tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint)              &
+                                                       - 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
+                                                                + tmp_J(ipoint,2) * tmp_M(ipoint,2) &
+                                                                + tmp_J(ipoint,3) * tmp_M(ipoint,3)))
+    enddo
+
+    noL_0e = noL_0e -2.d0 * (sum(tmp))
+
+    deallocate(tmp)
+
+  else
+
+    allocate(tmp(elec_alpha_num))
+    allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
+
+    !$OMP PARALLEL                                                   &
+    !$OMP DEFAULT(NONE)                                              &
+    !$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R)                        &
+    !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+    !$OMP        int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
+
+    !$OMP DO
+    do j = 1, elec_beta_num
+
+      tmp_L = 0.d0
+      tmp_R = 0.d0
+      do i = elec_beta_num+1, elec_alpha_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_L(ipoint,1) = tmp_L(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,2) = tmp_L(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,3) = tmp_L(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+
+          tmp_R(ipoint,1) = tmp_R(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,2) = tmp_R(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,3) = tmp_R(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+        enddo
+      enddo
+
+      tmp(j) = 0.d0
+      do ipoint = 1, n_points_final_grid
+        tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+      enddo
+
+      do i = 1, elec_beta_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+
+          tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+        enddo
+      enddo
+
+      do ipoint = 1, n_points_final_grid
+        tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+      enddo
+    enddo ! j
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+    ! ---
+
+    !$OMP PARALLEL                                                   &
+    !$OMP DEFAULT(NONE)                                              &
+    !$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R)                        &
+    !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+    !$OMP        int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
+
+    !$OMP DO
+    do j = elec_beta_num+1, elec_alpha_num
+
+      tmp_L = 0.d0
+      tmp_R = 0.d0
+      do i = 1, elec_alpha_num
+        do ipoint = 1, n_points_final_grid
+          tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+          tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+
+          tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+        enddo
+      enddo
+
+      tmp(j) = 0.d0
+      do ipoint = 1, n_points_final_grid
+        tmp(j) = tmp(j) + 0.5d0 * final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+      enddo
+    enddo ! j
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+    noL_0e = -2.d0 * sum(tmp)
+
+    deallocate(tmp)
+    deallocate(tmp_L, tmp_R)
+
+    ! ---
+
+    allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
+    tmp_O = 0.d0
+    tmp_J = 0.d0
+
+    !$OMP PARALLEL                                                   &
+    !$OMP DEFAULT(NONE)                                              &
+    !$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv)                 &
+    !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+    !$OMP        int2_grad1_u12_bimo_t, tmp_O, tmp_J)
+
+    allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
+    tmp_O_priv = 0.d0
+    tmp_J_priv = 0.d0
+  
+    !$OMP DO 
+    do i = 1, elec_beta_num
+      do ipoint = 1, n_points_final_grid
+        tmp_O_priv(ipoint)   = tmp_O_priv(ipoint)   + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+        tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
+        tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
+        tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP DO 
+    do i = elec_beta_num+1, elec_alpha_num
+      do ipoint = 1, n_points_final_grid
+        tmp_O_priv(ipoint)   = tmp_O_priv(ipoint)   + 0.5d0 * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+        tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,i)
+        tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,i)
+        tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,i)
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP CRITICAL
+    tmp_O = tmp_O + tmp_O_priv
+    tmp_J = tmp_J + tmp_J_priv
+    !$OMP END CRITICAL
+
+    deallocate(tmp_O_priv, tmp_J_priv)
+    !$OMP END PARALLEL
+
+    ! ---
+
+    allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
+    tmp_M = 0.d0
+    tmp_S = 0.d0
+
+    !$OMP PARALLEL                                                   &
+    !$OMP DEFAULT(NONE)                                              &
+    !$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv)              &
+    !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+    !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+    !$OMP        int2_grad1_u12_bimo_t, tmp_M, tmp_S)
+
+    allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
+    tmp_M_priv = 0.d0
+    tmp_S_priv = 0.d0
+  
+    !$OMP DO COLLAPSE(2)
+    do i = 1, elec_beta_num
+      do j = 1, elec_beta_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+
+          tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                      + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                      + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP DO COLLAPSE(2)
+    do i = elec_beta_num+1, elec_alpha_num
+      do j = 1, elec_beta_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+
+          tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
+          tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
+
+          tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                      + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                      + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP DO COLLAPSE(2)
+    do i = elec_beta_num+1, elec_alpha_num
+      do j = elec_beta_num+1, elec_alpha_num
+        do ipoint = 1, n_points_final_grid
+
+          tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+          tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+
+          tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                      + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                      + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO NOWAIT
+
+    !$OMP CRITICAL
+    tmp_M = tmp_M + tmp_M_priv
+    tmp_S = tmp_S + tmp_S_priv
+    !$OMP END CRITICAL
+
+    deallocate(tmp_M_priv, tmp_S_priv)
+    !$OMP END PARALLEL
+
+    allocate(tmp(n_points_final_grid))
+
+    do ipoint = 1, n_points_final_grid
+
+      tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
+
+      tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint)              &
+                                                       - 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
+                                                                + tmp_J(ipoint,2) * tmp_M(ipoint,2) &
+                                                                + tmp_J(ipoint,3) * tmp_M(ipoint,3)))
+    enddo
+
+    noL_0e = noL_0e -2.d0 * (sum(tmp))
+
+    deallocate(tmp)
+
+  endif
+
+  call wall_time(t1)
+  print*, " Wall time for noL_0e (min) = ", (t1 - t0)/60.d0
+
+  print*, " noL_0e = ", noL_0e
+
+END_PROVIDER
+
+! ---
+
 BEGIN_PROVIDER [double precision, noL_1e, (mo_num, mo_num)]
 
   implicit none
@@ -1028,7 +1423,7 @@ BEGIN_PROVIDER [double precision, noL_2e, (mo_num, mo_num, mo_num, mo_num)]
 
     call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 4*n_points_final_grid, 0.5d0       &
               , tmp1(1,1,1,1), 4*n_points_final_grid, tmp2(1,1,1,1), 4*n_points_final_grid &
-              , 1.d0, tmp(1,1,1,1), mo_num*mo_num)
+              , 0.d0, tmp(1,1,1,1), mo_num*mo_num)
 
     deallocate(tmp1, tmp2)
 
@@ -1178,7 +1573,7 @@ BEGIN_PROVIDER [double precision, noL_2e, (mo_num, mo_num, mo_num, mo_num)]
 
     call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 4*n_points_final_grid, 0.5d0       &
               , tmp1(1,1,1,1), 4*n_points_final_grid, tmp2(1,1,1,1), 4*n_points_final_grid &
-              , 1.d0, tmp(1,1,1,1), mo_num*mo_num)
+              , 0.d0, tmp(1,1,1,1), mo_num*mo_num)
 
     deallocate(tmp1, tmp2)
 

src/tc_bi_ortho/test_tc_bi_ortho.irp.f

@@ -37,9 +37,10 @@ program tc_bi_ortho
   !call test_no()
   !call test_no_v0()
 
-  !call test_no_0()
+  call test_noL_0e()
-  call test_no_1()
+  call test_noL_1e()
-  call test_no_2()
+  !call test_noL_2e_v0()
+  call test_noL_2e()
 
 end
 
@@ -319,7 +320,7 @@ subroutine test_no_v0()
   print*, ' accu (%) = ', 100.d0*accu/norm
 
   return
-end subroutine test_no_0
+end
 
 ! ---
 
@@ -365,7 +366,7 @@ subroutine test_no()
   print*, ' accu (%) = ', 100.d0*accu/norm
 
   return
-end subroutine test_no
+end
 
 ! ---
 
@@ -502,19 +503,28 @@ end
 
 ! ---
 
-subroutine test_no_0()
+subroutine test_noL_0e()
 
   implicit none
-  double precision :: accu, norm
+  double precision :: accu, norm, thr
 
-  print*, ' testing no_0 ...'
+  thr = 1d-8
+
+  print*, ' testing noL_0e ...'
 
   PROVIDE noL_0e_naive
+  PROVIDE noL_0e_v0
   PROVIDE noL_0e
 
   accu = dabs(noL_0e_naive - noL_0e)
   norm = dabs(noL_0e_naive)
 
+  if(accu .gt. thr) then
+    print*, ' problem on noL_0e'
+    print*, noL_0e_naive, noL_0e
+    stop
+  endif
+
   print*, ' accu (%) = ', 100.d0*accu/norm
 
   return
@@ -522,16 +532,17 @@ end
 
 ! ---
 
-subroutine test_no_1()
+subroutine test_noL_1e()
 
   implicit none
   integer          :: i, j
   double precision :: accu, contrib, new, ref, thr, norm
 
-  print*, ' testing no_1 ...'
+  print*, ' testing noL_1e ...'
 
   PROVIDE noL_1e_naive
   PROVIDE noL_1e
+  PROVIDE energy_1e_noL_HF
 
   thr = 1d-8
 
@@ -557,24 +568,68 @@ subroutine test_no_1()
 
   print*, ' accu (%) = ', 100.d0*accu/norm
 
-  PROVIDE energy_1e_noL_HF
+  return
+end
+
+! ---
+
+subroutine test_noL_2e_v0()
+
+  implicit none
+  integer          :: i, j, k, l
+  double precision :: accu, contrib, new, ref, thr, norm
+
+  print*, ' testing noL_2e_v0 ...'
+
+  PROVIDE noL_2e_naive
+  PROVIDE noL_2e_v0
+  PROVIDE energy_2e_noL_HF
+
+  thr = 1d-8
+
+  accu = 0.d0
+  norm = 0.d0
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do k = 1, mo_num
+        do l = 1, mo_num
+
+          new = noL_2e_v0   (l,k,j,i)
+          ref = noL_2e_naive(l,k,j,i)
+          contrib = dabs(new - ref)
+          if(contrib .gt. thr) then
+            print*, ' problem on noL_2e_v0'
+            print*, l, k, j, i
+            print*, ref, new, contrib
+            stop
+          endif
+          
+          accu += contrib
+          norm += dabs(ref)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  print*, ' accu (%) = ', 100.d0*accu/norm
 
   return
 end
 
 ! ---
 
-subroutine test_no_2()
+
+subroutine test_noL_2e()
 
   implicit none
   integer          :: i, j, k, l
   double precision :: accu, contrib, new, ref, thr, norm
 
-  print*, ' testing no_2 ...'
+  print*, ' testing noL_2e ...'
 
   PROVIDE noL_2e_naive
   PROVIDE noL_2e
-  !PROVIDE energy_2e_noL_HF
+  PROVIDE energy_2e_noL_HF
 
   thr = 1d-8
 

src/tc_scf/fock_three_hermit.irp.f

@@ -78,13 +78,16 @@ end
 
 ! ---
 
-! TODO DGEMM
 BEGIN_PROVIDER [double precision, diag_three_elem_hf]
 
   implicit none
-  integer          :: i, j, k, ipoint, mm
+  integer                       :: i, j, k, ipoint, mm
-  double precision :: contrib, weight, four_third, one_third, two_third, exchange_int_231
+  double precision              :: contrib, weight, four_third, one_third, two_third, exchange_int_231
-  double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb
+  double precision              :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb
+  double precision, allocatable :: tmp(:)
+  double precision, allocatable :: tmp_L(:,:), tmp_R(:,:)
+  double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:)
+  double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:)
 
   PROVIDE mo_l_coef mo_r_coef
 
@@ -131,14 +134,397 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf]
 
     else
 
-      provide mo_l_coef mo_r_coef
+      ! ------------
-      call give_aaa_contrib(integral_aaa)
+      ! SLOW VERSION
-      call give_aab_contrib(integral_aab)
+      ! ------------
-      call give_abb_contrib(integral_abb)
+
-      call give_bbb_contrib(integral_bbb)
+      !call give_aaa_contrib(integral_aaa)
-      diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb
+      !call give_aab_contrib(integral_aab)
-!      print*,'integral_aaa + integral_aab + integral_abb + integral_bbb'
+      !call give_abb_contrib(integral_abb)
-!      print*,integral_aaa , integral_aab , integral_abb , integral_bbb
+      !call give_bbb_contrib(integral_bbb)
+      !diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb
+
+      ! ------------
+      ! ------------
+
+      PROVIDE int2_grad1_u12_bimo_t
+      PROVIDE mos_l_in_r_array_transp
+      PROVIDE mos_r_in_r_array_transp
+
+      if(elec_alpha_num .eq. elec_beta_num) then
+    
+        allocate(tmp(elec_beta_num))
+        allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
+    
+        !$OMP PARALLEL                                                 &
+        !$OMP DEFAULT(NONE)                                            &
+        !$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R)                      &
+        !$OMP SHARED(elec_beta_num, n_points_final_grid,               & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
+        !$OMP        int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
+    
+        !$OMP DO
+        do j = 1, elec_beta_num
+    
+          tmp_L = 0.d0
+          tmp_R = 0.d0
+          do i = 1, elec_beta_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+    
+              tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+            enddo
+          enddo
+    
+          tmp(j) = 0.d0
+          do ipoint = 1, n_points_final_grid
+            tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+          enddo
+        enddo ! j
+        !$OMP END DO
+        !$OMP END PARALLEL
+    
+        diag_three_elem_hf = -2.d0 * sum(tmp)
+    
+        deallocate(tmp)
+        deallocate(tmp_L, tmp_R)
+    
+        ! ---
+    
+        allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
+        tmp_O = 0.d0
+        tmp_J = 0.d0
+    
+        !$OMP PARALLEL                                                  &
+        !$OMP DEFAULT(NONE)                                             &
+        !$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv)                &
+        !$OMP SHARED(elec_beta_num, n_points_final_grid,                & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,  &
+        !$OMP        int2_grad1_u12_bimo_t, tmp_O, tmp_J)
+    
+        allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
+        tmp_O_priv = 0.d0
+        tmp_J_priv = 0.d0
+      
+        !$OMP DO 
+        do i = 1, elec_beta_num
+          do ipoint = 1, n_points_final_grid
+            tmp_O_priv(ipoint)   = tmp_O_priv(ipoint)   + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+            tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
+            tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
+            tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP CRITICAL
+        tmp_O = tmp_O + tmp_O_priv
+        tmp_J = tmp_J + tmp_J_priv
+        !$OMP END CRITICAL
+    
+        deallocate(tmp_O_priv, tmp_J_priv)
+        !$OMP END PARALLEL
+    
+        allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
+        tmp_M = 0.d0
+        tmp_S = 0.d0
+    
+        !$OMP PARALLEL                                                 &
+        !$OMP DEFAULT(NONE)                                            &
+        !$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv)            &
+        !$OMP SHARED(elec_beta_num, n_points_final_grid,               & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
+        !$OMP        int2_grad1_u12_bimo_t, tmp_M, tmp_S)
+    
+        allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
+        tmp_M_priv = 0.d0
+        tmp_S_priv = 0.d0
+      
+        !$OMP DO COLLAPSE(2)
+        do i = 1, elec_beta_num
+          do j = 1, elec_beta_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+    
+              tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                          + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                          + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+            enddo
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP CRITICAL
+        tmp_M = tmp_M + tmp_M_priv
+        tmp_S = tmp_S + tmp_S_priv
+        !$OMP END CRITICAL
+    
+        deallocate(tmp_M_priv, tmp_S_priv)
+        !$OMP END PARALLEL
+    
+        allocate(tmp(n_points_final_grid))
+    
+        do ipoint = 1, n_points_final_grid
+    
+          tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
+    
+          tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint)              &
+                                                           - 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
+                                                                    + tmp_J(ipoint,2) * tmp_M(ipoint,2) &
+                                                                    + tmp_J(ipoint,3) * tmp_M(ipoint,3)))
+        enddo
+    
+        diag_three_elem_hf = diag_three_elem_hf -2.d0 * (sum(tmp))
+    
+        deallocate(tmp)
+    
+      else
+    
+        allocate(tmp(elec_alpha_num))
+        allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
+    
+        !$OMP PARALLEL                                                   &
+        !$OMP DEFAULT(NONE)                                              &
+        !$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R)                        &
+        !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+        !$OMP        int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
+    
+        !$OMP DO
+        do j = 1, elec_beta_num
+    
+          tmp_L = 0.d0
+          tmp_R = 0.d0
+          do i = elec_beta_num+1, elec_alpha_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_L(ipoint,1) = tmp_L(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,2) = tmp_L(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,3) = tmp_L(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+    
+              tmp_R(ipoint,1) = tmp_R(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,2) = tmp_R(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,3) = tmp_R(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+            enddo
+          enddo
+    
+          tmp(j) = 0.d0
+          do ipoint = 1, n_points_final_grid
+            tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+          enddo
+    
+          do i = 1, elec_beta_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+    
+              tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+            enddo
+          enddo
+    
+          do ipoint = 1, n_points_final_grid
+            tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+          enddo
+        enddo ! j
+        !$OMP END DO
+        !$OMP END PARALLEL
+    
+        ! ---
+    
+        !$OMP PARALLEL                                                   &
+        !$OMP DEFAULT(NONE)                                              &
+        !$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R)                        &
+        !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+        !$OMP        int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
+    
+        !$OMP DO
+        do j = elec_beta_num+1, elec_alpha_num
+    
+          tmp_L = 0.d0
+          tmp_R = 0.d0
+          do i = 1, elec_alpha_num
+            do ipoint = 1, n_points_final_grid
+              tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
+    
+              tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
+            enddo
+          enddo
+    
+          tmp(j) = 0.d0
+          do ipoint = 1, n_points_final_grid
+            tmp(j) = tmp(j) + 0.5d0 * final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
+          enddo
+        enddo ! j
+        !$OMP END DO
+        !$OMP END PARALLEL
+    
+        diag_three_elem_hf = -2.d0 * sum(tmp)
+    
+        deallocate(tmp)
+        deallocate(tmp_L, tmp_R)
+    
+        ! ---
+    
+        allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
+        tmp_O = 0.d0
+        tmp_J = 0.d0
+    
+        !$OMP PARALLEL                                                   &
+        !$OMP DEFAULT(NONE)                                              &
+        !$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv)                 &
+        !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+        !$OMP        int2_grad1_u12_bimo_t, tmp_O, tmp_J)
+    
+        allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
+        tmp_O_priv = 0.d0
+        tmp_J_priv = 0.d0
+      
+        !$OMP DO 
+        do i = 1, elec_beta_num
+          do ipoint = 1, n_points_final_grid
+            tmp_O_priv(ipoint)   = tmp_O_priv(ipoint)   + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+            tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
+            tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
+            tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP DO 
+        do i = elec_beta_num+1, elec_alpha_num
+          do ipoint = 1, n_points_final_grid
+            tmp_O_priv(ipoint)   = tmp_O_priv(ipoint)   + 0.5d0 * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+            tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,i)
+            tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,i)
+            tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,i)
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP CRITICAL
+        tmp_O = tmp_O + tmp_O_priv
+        tmp_J = tmp_J + tmp_J_priv
+        !$OMP END CRITICAL
+    
+        deallocate(tmp_O_priv, tmp_J_priv)
+        !$OMP END PARALLEL
+    
+        ! ---
+    
+        allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
+        tmp_M = 0.d0
+        tmp_S = 0.d0
+    
+        !$OMP PARALLEL                                                   &
+        !$OMP DEFAULT(NONE)                                              &
+        !$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv)              &
+        !$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, & 
+        !$OMP        mos_l_in_r_array_transp, mos_r_in_r_array_transp,   &
+        !$OMP        int2_grad1_u12_bimo_t, tmp_M, tmp_S)
+    
+        allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
+        tmp_M_priv = 0.d0
+        tmp_S_priv = 0.d0
+      
+        !$OMP DO COLLAPSE(2)
+        do i = 1, elec_beta_num
+          do j = 1, elec_beta_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+    
+              tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                          + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                          + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+            enddo
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP DO COLLAPSE(2)
+        do i = elec_beta_num+1, elec_alpha_num
+          do j = 1, elec_beta_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+    
+              tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
+              tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
+    
+              tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                          + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                          + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+            enddo
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP DO COLLAPSE(2)
+        do i = elec_beta_num+1, elec_alpha_num
+          do j = elec_beta_num+1, elec_alpha_num
+            do ipoint = 1, n_points_final_grid
+    
+              tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+              tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
+    
+              tmp_S_priv(ipoint)   = tmp_S_priv(ipoint)   + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+                                                          + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+                                                          + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
+            enddo
+          enddo
+        enddo
+        !$OMP END DO NOWAIT
+    
+        !$OMP CRITICAL
+        tmp_M = tmp_M + tmp_M_priv
+        tmp_S = tmp_S + tmp_S_priv
+        !$OMP END CRITICAL
+    
+        deallocate(tmp_M_priv, tmp_S_priv)
+        !$OMP END PARALLEL
+    
+        allocate(tmp(n_points_final_grid))
+    
+        do ipoint = 1, n_points_final_grid
+    
+          tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
+    
+          tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint)              &
+                                                           - 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
+                                                                    + tmp_J(ipoint,2) * tmp_M(ipoint,2) &
+                                                                    + tmp_J(ipoint,3) * tmp_M(ipoint,3)))
+        enddo
+    
+        diag_three_elem_hf = diag_three_elem_hf - 2.d0 * (sum(tmp))
+    
+        deallocate(tmp)
+    
+      endif
+
 
     endif
 
@@ -374,3 +760,7 @@ BEGIN_PROVIDER [ double precision, fock_3_w_kl_wla_phi_k, (n_points_final_grid,3
  enddo
 END_PROVIDER 
 
+
+
+
+