5idx arrays need O(N3) temp arrays

src/bi_ort_ints/bi_ort_ints.irp.f

@@ -8,8 +8,8 @@ program bi_ort_ints
   my_n_pt_a_grid = 14
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 ! call test_3e
+  call test_5idx2
  call test_5idx
-!  call test_5idx2
 end
 
 subroutine test_5idx2
@@ -60,6 +60,8 @@ subroutine test_5idx
  k = 1
  n = 0
  accu = 0.d0
+ PROVIDE three_e_5_idx_direct_bi_ort_old
+
  do i = 1, mo_num
   do k = 1, mo_num
    do j = 1, mo_num
@@ -69,28 +71,28 @@ subroutine test_5idx
 !         stop
 !      endif
 
-!      new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
+      new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
-!      ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
+      ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
+      contrib = dabs(new - ref)
+      accu += contrib
+      if(contrib .gt. 1.d-10)then
+       print*,'direct'
+       print*,i,k,j,l,m
+       print*,ref,new,contrib
+       stop
+      endif
+!
+!      new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
+!      ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
 !      contrib = dabs(new - ref)
 !      accu += contrib
 !      if(contrib .gt. 1.d-10)then
-!       print*,'direct'
+!       print*,'exch12'
 !       print*,i,k,j,l,m
 !       print*,ref,new,contrib
 !       stop
 !      endif
 !
-      new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
-      accu += contrib
-      if(contrib .gt. 1.d-10)then
-       print*,'exch12'
-       print*,i,k,j,l,m
-       print*,ref,new,contrib
-       stop
-      endif
-
 !
 !      new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
 !      ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
@@ -135,7 +137,7 @@ subroutine test_5idx
 !       print*,ref,new,contrib
 !       stop
 !      endif
-!
+
      enddo
     enddo
    enddo

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -15,7 +15,7 @@ end
   !
   ! matrix element of the -L  three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
   !
-  ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
+  ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> :: : notice that i is the RIGHT MO and k is the LEFT MO
   !
   ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
   !
@@ -26,78 +26,97 @@ end
   integer                        :: i, j, k, m, l
   double precision               :: wall1, wall0
   integer                        :: ipoint
-  double precision, allocatable :: grad_mli(:,:,:), orb_mat(:,:,:)
+  double precision, allocatable  :: grad_mli(:,:), orb_mat(:,:,:)
-  double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:)
+  double precision, allocatable  :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:)
-  double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:)
+  double precision, allocatable  :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:)
-  double precision, allocatable :: tmp_mat(:,:,:,:)
+  double precision, allocatable  :: tmp_mat(:,:,:)
-  allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
 
   provide mos_r_in_r_array_transp mos_l_in_r_array_transp
   PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
 
+  call print_memory_usage
   print *, ' Providing the three_e_5_idx_bi_ort ...'
   call wall_time(wall0)
 
- do m = 1, mo_num
+  three_e_5_idx_direct_bi_ort (:,:,:,:,:) = 0.d0
+  three_e_5_idx_cycle_1_bi_ort(:,:,:,:,:) = 0.d0
+  three_e_5_idx_cycle_2_bi_ort(:,:,:,:,:) = 0.d0
+  three_e_5_idx_exch23_bi_ort (:,:,:,:,:) = 0.d0
+  three_e_5_idx_exch13_bi_ort (:,:,:,:,:) = 0.d0
 
-  allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
+  call print_memory_usage
+
+  allocate(tmp_mat(mo_num,mo_num,mo_num))
   allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
- !$OMP PARALLEL                     &
+
- !$OMP DEFAULT (NONE)               &
+  !$OMP PARALLEL DO PRIVATE (i,l,ipoint)
- !$OMP PRIVATE (i,l,ipoint) &
- !$OMP SHARED (m,mo_num,n_points_final_grid, &
- !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
- !$OMP         int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
- !$OMP         grad_mli, orb_mat)
- !$OMP DO COLLAPSE(2)
   do i=1,mo_num
     do l=1,mo_num
       do ipoint=1, n_points_final_grid
 
-         grad_mli(ipoint,l,i) = final_weight_at_r_vector(ipoint) * ( &
+        orb_mat(ipoint,l,i) = final_weight_at_r_vector(ipoint)       &
-               int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) + &
+            * mos_l_in_r_array_transp(ipoint,l)                      &
-               int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) + &
+            * mos_r_in_r_array_transp(ipoint,i)
-               int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) )
 
-         orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
-
-       enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, 1.d0, &
-      orb_mat, n_points_final_grid,  &
-      grad_mli, n_points_final_grid,  0.d0, &
-      tmp_mat, mo_num*mo_num)
-
-  !$OMP PARALLEL DO PRIVATE(i,j,k,l)
-  do i = 1, mo_num
-    do k = 1, mo_num
-      do j = 1, mo_num
-        do l = 1, mo_num
-            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j)
-        enddo
       enddo
     enddo
   enddo
   !$OMP END PARALLEL DO
 
-  deallocate(orb_mat,grad_mli)
+  tmp_mat = 0.d0
+  call print_memory_usage
+!
+  do m = 1, mo_num
+
+    allocate(grad_mli(n_points_final_grid,mo_num))
+
+    do i=1,mo_num
+      !$OMP PARALLEL DO PRIVATE (l,ipoint)
+      do l=1,mo_num
+        do ipoint=1, n_points_final_grid
+
+          grad_mli(ipoint,l) =                                       &
+              int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) +&
+              int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) +&
+              int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
+
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, n_points_final_grid, 1.d0,&
+          orb_mat, n_points_final_grid,                              &
+          grad_mli, n_points_final_grid,  0.d0,                      &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL PRIVATE(j,k,l)
+      !$OMP DO
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do l = 1, mo_num
+            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k)
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP DO
+      do j = 1, mo_num
+        do l = 1, mo_num
+          do k = 1, mo_num
+            three_e_5_idx_direct_bi_ort(m,k,i,l,j) = three_e_5_idx_direct_bi_ort(m,k,i,l,j) - tmp_mat(l,j,k)
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+    enddo
+
+    deallocate(grad_mli)
 
     allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
-  allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
+    allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
-  allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
 
- !$OMP PARALLEL                     &
+    !$OMP PARALLEL DO PRIVATE (i,l,ipoint)
- !$OMP DEFAULT (NONE)               &
- !$OMP PRIVATE (i,l,ipoint) &
- !$OMP SHARED (m,mo_num,n_points_final_grid, &
- !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
- !$OMP         int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
- !$OMP         rm_grad_ik, lm_grad_ik, rk_grad_im, lk_grad_mi)
- !$OMP DO COLLAPSE(2)
     do i=1,mo_num
       do l=1,mo_num
         do ipoint=1, n_points_final_grid
@@ -106,72 +125,6 @@ end
           lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint)
           lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint)
 
-         rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
-         rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
-         rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
-
-         rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
-         rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
-         rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
-
-       enddo
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-  call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
-      lm_grad_ik, 3*n_points_final_grid,  &
-      rm_grad_ik, 3*n_points_final_grid,  0.d0, &
-      tmp_mat, mo_num*mo_num)
-
-  !$OMP PARALLEL DO PRIVATE(i,j,k,l)
-  do i = 1, mo_num
-    do k = 1, mo_num
-      do j = 1, mo_num
-        do l = 1, mo_num
-            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END PARALLEL DO
-
-  call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
-      lm_grad_ik, 3*n_points_final_grid,  &
-      rk_grad_im, 3*n_points_final_grid,  0.d0, &
-      tmp_mat, mo_num*mo_num)
-
-  !$OMP PARALLEL DO PRIVATE(i,j,k,l)
-  do i = 1, mo_num
-    do k = 1, mo_num
-      do j = 1, mo_num
-        do l = 1, mo_num
-            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,j,k)
-            three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = - tmp_mat(k,j,i,l)
-            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = - tmp_mat(k,i,j,l)
-            three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = - tmp_mat(l,j,i,k)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$OMP END PARALLEL DO
-
-  deallocate(lm_grad_ik)
-
-  allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
-
- !$OMP PARALLEL                     &
- !$OMP DEFAULT (NONE)               &
- !$OMP PRIVATE (i,l,ipoint) &
- !$OMP SHARED (m,mo_num,n_points_final_grid, &
- !$OMP         mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
- !$OMP         int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
- !$OMP         lk_grad_mi)
- !$OMP DO COLLAPSE(2)
-  do i=1,mo_num
-    do l=1,mo_num
-       do ipoint=1, n_points_final_grid
-
           lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint)
           lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint)
           lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint)
@@ -179,54 +132,109 @@ end
         enddo
       enddo
     enddo
-  !$OMP END DO
+    !$OMP END PARALLEL DO
-  !$OMP END PARALLEL
 
-  call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
+
+    allocate(rm_grad_ik(n_points_final_grid,3,mo_num))
+    allocate(rk_grad_im(n_points_final_grid,3,mo_num))
+
+    do i=1,mo_num
+      !$OMP PARALLEL DO PRIVATE (l,ipoint)
+      do l=1,mo_num
+        do ipoint=1, n_points_final_grid
+
+          rm_grad_ik(ipoint,1,l) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
+          rm_grad_ik(ipoint,2,l) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
+          rm_grad_ik(ipoint,3,l) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
+
+          rk_grad_im(ipoint,1,l) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
+          rk_grad_im(ipoint,2,l) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
+          rk_grad_im(ipoint,3,l) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
+
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
+          lm_grad_ik, 3*n_points_final_grid,                         &
+          rm_grad_ik, 3*n_points_final_grid,  0.d0,                  &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do l = 1, mo_num
+            three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k)
+          enddo
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
+          lm_grad_ik, 3*n_points_final_grid,                         &
+          rk_grad_im, 3*n_points_final_grid,  0.d0,                  &
+          tmp_mat, mo_num*mo_num)
+
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do l = 1, mo_num
+            three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) = three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) - tmp_mat(l,k,j)
+            three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) = three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) - tmp_mat(k,j,l)
+            three_e_5_idx_exch23_bi_ort (m,i,j,k,l) = three_e_5_idx_exch23_bi_ort (m,i,j,k,l) - tmp_mat(k,l,j)
+            three_e_5_idx_exch13_bi_ort (m,l,j,i,k) = three_e_5_idx_exch13_bi_ort (m,l,j,i,k) - tmp_mat(l,j,k)
+          enddo
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
           lk_grad_mi, 3*n_points_final_grid,                         &
           rm_grad_ik, 3*n_points_final_grid,  0.d0,                  &
           tmp_mat, mo_num*mo_num)
 
-  !$OMP PARALLEL DO PRIVATE(i,j,k,l)
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
-  do i = 1, mo_num
       do k = 1, mo_num
         do j = 1, mo_num
           do l = 1, mo_num
-            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(k,j,l,i)
+            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(k,j,l)
-            three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
+            three_e_5_idx_cycle_2_bi_ort(m,l,i,k,j) = three_e_5_idx_cycle_2_bi_ort(m,l,i,k,j) - tmp_mat(l,j,k)
-            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(l,j,k,i)
+            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(l,j,k)
-            three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(k,i,l,j)
+            three_e_5_idx_exch13_bi_ort (m,l,i,k,j) = three_e_5_idx_exch13_bi_ort (m,l,i,k,j) - tmp_mat(k,j,l)
-        enddo
           enddo
         enddo
       enddo
       !$OMP END PARALLEL DO
 
-  call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
+      call dgemm('T','N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0,&
           lk_grad_mi, 3*n_points_final_grid,                         &
           rk_grad_im, 3*n_points_final_grid,  0.d0,                  &
           tmp_mat, mo_num*mo_num)
 
-  !$OMP PARALLEL DO PRIVATE(i,j,k,l)
+      !$OMP PARALLEL DO PRIVATE(j,k,l)
-  do i = 1, mo_num
       do k = 1, mo_num
         do j = 1, mo_num
           do l = 1, mo_num
-            three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(l,j,i,k)
+            three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) = three_e_5_idx_cycle_1_bi_ort(m,l,j,i,k) - tmp_mat(l,j,k)
-            three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(k,i,j,l)
+            three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) = three_e_5_idx_cycle_2_bi_ort(m,i,j,k,l) - tmp_mat(k,l,j)
-            three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(k,j,i,l)
+            three_e_5_idx_exch23_bi_ort (m,i,j,k,l) = three_e_5_idx_exch23_bi_ort (m,i,j,k,l) - tmp_mat(k,j,l)
-            three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(l,i,j,k)
+            three_e_5_idx_exch13_bi_ort (m,l,j,i,k) = three_e_5_idx_exch13_bi_ort (m,l,j,i,k) - tmp_mat(l,k,j)
-        enddo
           enddo
         enddo
       enddo
       !$OMP END PARALLEL DO
 
-  deallocate(lk_grad_mi)
+    enddo
     deallocate(rm_grad_ik)
     deallocate(rk_grad_im)
+    deallocate(lk_grad_mi)
+    deallocate(lm_grad_ik)
+
   enddo
 
+  deallocate(orb_mat)
+
   call wall_time(wall1)
   print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0