normal ordering: aaa-DGEMM OK

src/tc_bi_ortho/normal_ordered.irp.f

@@ -47,12 +47,13 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
 
     PROVIDE no_aba_contraction
     PROVIDE no_aab_contraction
+    PROVIDE no_aaa_contraction
 
     !$OMP PARALLEL                                                              &
     !$OMP DEFAULT (NONE)                                                        &
     !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aab, hthree_aaa)  & 
     !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth, &
-    !$OMP         no_aba_contraction,no_aab_contraction)
+    !$OMP         no_aba_contraction, no_aab_contraction, no_aaa_contraction)
     !$OMP DO SCHEDULE (static) 
     do hh1 = 1, n_act_orb
       h1 = list_act(hh1) 
@@ -66,33 +67,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
           do pp2 = 1, n_act_orb
             p2 = list_act(pp2)
 
-            ! all contributions from the 3-e terms to the double excitations 
+            normal_two_body_bi_orth(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + no_aab_contraction(p2,h2,p1,h1) + no_aaa_contraction(p2,h2,p1,h1)
-            ! s1:(h1-->p1), s2:(h2-->p2) from the HF reference determinant 
-
-            ! same spin double excitations : s1 == s2 
-            if((h1 < h2) .and. (p1 > p2)) then
-
-              ! same spin double excitations with same spin contributions 
-              if(Ne(2) .ge. 3) then
-                call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
-              else
-                hthree_aaa = 0.d0
-              endif
-
-            else
-
-              if(Ne(2) .ge. 3) then
-                ! same spin double excitations with same spin contributions 
-                call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
-              else
-                hthree_aaa = 0.d0
-              endif
-
-            endif
-
-            normal_two_body_bi_orth(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) &
-                                                 + no_aab_contraction(p2,h2,p1,h1) &
-                                                 + 0.5d0 * hthree_aaa
           enddo
         enddo
       enddo
@@ -616,3 +591,565 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
 END_PROVIDER
 
 ! ---
+
+BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_num)]
+
+  BEGIN_DOC
+  !
+  ! if:
+  !    h1 < h2
+  !    p1 > p2
+  !
+  !   no_aaa_contraction(p2,h2.p1,h1) =  0.5 [Ialpha(p2,h1,p1,h2) + Ibeta(p2,h1,p1,h2)]
+  !                                   = -0.5 [Ialpha(p2,h2,p1,h1) + Ibeta(p2,h2,p1,h1)]
+  !
+  ! else:
+  !
+  !   no_aaa_contraction(p2,h2.p1,h1) = 0.5 [Ialpha(p2,h2,p1,h1) + Ibeta(p2,h2,p1,h1)]
+  !
+  ! 
+  ! I(p2,h2,p1,h1) = J(p2,h2,p1,h1) - J(p1,h2,p2,h1)
+  ! J(p2,h2,p1,h1) = \sum_i [ <  i p2 p1 | i h2 h1 >
+  !                         + < p2 p1  i | i h2 h1 >
+  !                         + < p1  i p2 | i h2 h1 > ]
+  !
+  !
+  END_DOC
+
+  use bitmasks ! you need to include the bitmasks_module.f90 features
+
+  implicit none
+  integer                        :: i, ii, h1, p1, h2, p2, ipoint
+  integer                        :: Ne(2)
+  double precision               :: wall0, wall1
+  integer,           allocatable :: occ(:,:)
+  integer(bit_kind), allocatable :: key_i_core(:,:)
+  double precision,  allocatable :: tmp_2d(:,:), tmp_3d(:,:,:)
+  double precision,  allocatable :: tmp1(:,:,:), tmp2(:,:), tmp3(:,:,:)
+  double precision,  allocatable :: tmpval_1(:), tmpval_2(:), tmpvec_1(:,:), tmpvec_2(:,:), tmpvec_3(:,:)
+
+  print*,' Providing no_aaa_contraction ...'
+  call wall_time(wall0)
+
+  PROVIDE N_int
+
+  allocate(occ(N_int*bit_kind_size,2))
+  allocate(key_i_core(N_int,2))
+
+  if(core_tc_op) then
+    do i = 1, N_int
+      key_i_core(i,1) = xor(ref_bitmask(i,1), core_bitmask(i,1))
+      key_i_core(i,2) = xor(ref_bitmask(i,2), core_bitmask(i,2))
+    enddo
+    call bitstring_to_list_ab(key_i_core, occ, Ne, N_int)
+  else
+    call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
+  endif
+
+  if(Ne(2) .lt. 3) then
+
+    no_aaa_contraction = 0.d0
+
+  else
+
+    allocate(tmp_2d(mo_num,mo_num))
+    allocate(tmp_3d(mo_num,mo_num,mo_num))
+    allocate(tmp1(n_points_final_grid,3,mo_num))
+    allocate(tmp2(n_points_final_grid,mo_num))
+    allocate(tmp3(n_points_final_grid,3,mo_num))
+    allocate(tmpval_1(n_points_final_grid))
+    allocate(tmpval_2(n_points_final_grid))
+    allocate(tmpvec_1(n_points_final_grid,3))
+    allocate(tmpvec_2(n_points_final_grid,3))
+    allocate(tmpvec_3(n_points_final_grid,3))
+
+    ! purely closed shell part 
+    do ii = 1, Ne(2)
+      i = occ(ii,2)
+
+      ! to avoid tmp(N^4)
+      do h1 = 1, mo_num
+
+        ! to minimize the number of operations
+        !$OMP PARALLEL                                                  &
+        !$OMP DEFAULT (NONE)                                            &
+        !$OMP PRIVATE (ipoint)                                          &
+        !$OMP SHARED (n_points_final_grid, i, h1,                       &
+        !$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         tmpval_1, tmpval_2, tmpvec_1, tmpvec_2 )
+        !$OMP DO
+        do ipoint = 1, n_points_final_grid
+
+          tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+
+          tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,h1)
+
+          tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,i) * mos_r_in_r_array_transp(ipoint,h1)
+          tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,i) * mos_r_in_r_array_transp(ipoint,h1)
+          tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,i) * mos_r_in_r_array_transp(ipoint,h1)
+
+          tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_r_in_r_array_transp(ipoint,i)
+          tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_r_in_r_array_transp(ipoint,i)
+          tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_r_in_r_array_transp(ipoint,i)
+        enddo
+        !$OMP END DO
+        !$OMP END PARALLEL
+
+        !$OMP PARALLEL                                                &
+        !$OMP DEFAULT (NONE)                                          &
+        !$OMP PRIVATE (p1, ipoint)                                    &
+        !$OMP SHARED (mo_num, n_points_final_grid, h1, i,             &
+        !$OMP         mos_l_in_r_array_transp, int2_grad1_u12_bimo_t, &
+        !$OMP         tmpval_1, tmpvec_1, tmp1)
+        !$OMP DO 
+        do p1 = 1, mo_num
+          do ipoint = 1, n_points_final_grid
+            tmp1(ipoint,1,p1) = mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1)
+            tmp1(ipoint,2,p1) = mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1)
+            tmp1(ipoint,3,p1) = mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1)
+          enddo
+        enddo
+        !$OMP END DO
+        !$OMP END PARALLEL
+
+        call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 &
+                  , int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid        &
+                  , tmp1(1,1,1), 3*n_points_final_grid                           &
+                  , 0.d0, tmp_3d(1,1,1), mo_num*mo_num)
+
+        !$OMP PARALLEL DO PRIVATE(p1,h2,p2)
+        do p1 = 1, mo_num
+          do h2 = 1, mo_num
+            do p2 = 1, mo_num
+              no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
+            enddo
+          enddo
+        enddo
+        !$OMP END PARALLEL DO
+
+        !$OMP PARALLEL                                                &
+        !$OMP DEFAULT (NONE)                                          &
+        !$OMP PRIVATE (p2, ipoint)                                    &
+        !$OMP SHARED (mo_num, n_points_final_grid, h1, i,             &
+        !$OMP         mos_l_in_r_array_transp, int2_grad1_u12_bimo_t, &
+        !$OMP         tmpval_2, tmpvec_2, tmp1)
+        !$OMP DO 
+        do p2 = 1, mo_num
+          do ipoint = 1, n_points_final_grid
+            tmp1(ipoint,1,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,1)
+            tmp1(ipoint,2,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,2)
+            tmp1(ipoint,3,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,3)
+          enddo
+        enddo
+        !$OMP END DO
+        !$OMP END PARALLEL
+
+        call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
+                  , tmp1(1,1,1), 3*n_points_final_grid                           &
+                  , int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid        &
+                  , 0.d0, tmp_3d(1,1,1), mo_num)
+
+        !$OMP PARALLEL DO PRIVATE(p1,h2,p2)
+        do p1 = 1, mo_num
+          do h2 = 1, mo_num
+            do p2 = 1, mo_num
+              no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,p1,h2)
+            enddo
+          enddo
+        enddo
+        !$OMP END PARALLEL DO
+
+        ! to avoid tmp(N^4)
+        do p1 = 1, mo_num
+
+          !$OMP PARALLEL                                                  &
+          !$OMP DEFAULT (NONE)                                            &
+          !$OMP PRIVATE (ipoint)                                          &
+          !$OMP SHARED (n_points_final_grid, i, h1, p1,                   &
+          !$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         tmpval_1, tmpval_2, tmpvec_1, tmpvec_2, tmpvec_3)
+          !$OMP DO
+          do ipoint = 1, n_points_final_grid
+
+            tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) *                                          &
+                             ( int2_grad1_u12_bimo_t(ipoint,1,i,i) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) &
+                             + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) &
+                             + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) )
+
+            tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i)
+
+            tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_r_in_r_array_transp(ipoint,h1)
+            tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_r_in_r_array_transp(ipoint,h1)
+            tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_r_in_r_array_transp(ipoint,h1)
+
+            tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp(ipoint,p1)
+            tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp(ipoint,p1)
+            tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp(ipoint,p1)
+
+            tmpvec_3(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_l_in_r_array_transp(ipoint,i)
+            tmpvec_3(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_l_in_r_array_transp(ipoint,i)
+            tmpvec_3(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_l_in_r_array_transp(ipoint,i)
+          enddo
+          !$OMP END DO
+          !$OMP END PARALLEL
+
+          !$OMP PARALLEL                                &
+          !$OMP DEFAULT (NONE)                          &
+          !$OMP PRIVATE (h2, ipoint)                    &
+          !$OMP SHARED (mo_num, n_points_final_grid, i, &
+          !$OMP         mos_r_in_r_array_transp,        &
+          !$OMP         int2_grad1_u12_bimo_t,          &
+          !$OMP         tmp1, tmp2, tmpval_1, tmpval_2, tmpvec_1)
+          !$OMP DO 
+          do h2 = 1, mo_num
+            do ipoint = 1, n_points_final_grid
+
+              tmp2(ipoint,h2) = mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint)     & 
+                              + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) &
+                              + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) &
+                              + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3)
+
+              tmp1(ipoint,1,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h2)
+              tmp1(ipoint,2,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h2)
+              tmp1(ipoint,3,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h2)
+
+            enddo
+          enddo
+          !$OMP END DO
+          !$OMP END PARALLEL
+
+          call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 &
+                    , mos_l_in_r_array_transp(1,1), n_points_final_grid   &
+                    , tmp2(1,1), n_points_final_grid                      &
+                    , 0.d0, tmp_2d(1,1), mo_num)
+
+          !$OMP PARALLEL DO PRIVATE(h2,p2)
+          do h2 = 1, mo_num
+            do p2 = 1, mo_num
+              no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
+            enddo
+          enddo
+          !$OMP END PARALLEL DO
+
+          !$OMP PARALLEL                                    &
+          !$OMP DEFAULT (NONE)                              &
+          !$OMP PRIVATE (p2, ipoint)                        &
+          !$OMP SHARED (mo_num, n_points_final_grid, i, h1, &
+          !$OMP         int2_grad1_u12_bimo_t,              &
+          !$OMP         tmpvec_2, tmpvec_3, tmp2, tmp3)
+          !$OMP DO 
+          do p2 = 1, mo_num
+            do ipoint = 1, n_points_final_grid
+
+              tmp2(ipoint,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) &
+                              + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) &
+                              + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) 
+
+              tmp3(ipoint,1,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,h1) 
+              tmp3(ipoint,2,p2) = int2_grad1_u12_bimo_t(ipoint,2,p2,h1) 
+              tmp3(ipoint,3,p2) = int2_grad1_u12_bimo_t(ipoint,3,p2,h1) 
+            enddo
+          enddo
+          !$OMP END DO
+          !$OMP END PARALLEL
+
+          call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 &
+                    , tmp2(1,1), n_points_final_grid                      &
+                    , mos_r_in_r_array_transp(1,1), n_points_final_grid   &
+                    , 0.d0, tmp_2d(1,1), mo_num)
+
+          call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 1.d0 &
+                    , tmp3(1,1,1), 3*n_points_final_grid                    &
+                    , tmp1(1,1,1), 3*n_points_final_grid                    &
+                    , 1.d0, tmp_2d(1,1), mo_num)
+
+          !$OMP PARALLEL DO PRIVATE(h2,p2)
+          do h2 = 1, mo_num
+            do p2 = 1, mo_num
+              no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
+            enddo
+          enddo
+          !$OMP END PARALLEL DO
+
+        enddo ! p1
+      enddo ! h1
+    enddo ! i
+
+
+
+    ! purely open-shell part 
+    if(Ne(2) < Ne(1)) then
+
+      do ii = Ne(2) + 1, Ne(1)
+        i = occ(ii,1)
+
+
+        ! to avoid tmp(N^4)
+        do h1 = 1, mo_num
+
+          ! to minimize the number of operations
+          !$OMP PARALLEL                                                  &
+          !$OMP DEFAULT (NONE)                                            &
+          !$OMP PRIVATE (ipoint)                                          &
+          !$OMP SHARED (n_points_final_grid, i, h1,                       &
+          !$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         tmpval_1, tmpval_2, tmpvec_1, tmpvec_2 )
+          !$OMP DO
+          do ipoint = 1, n_points_final_grid
+
+            tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
+
+            tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,h1)
+
+            tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,i) * mos_r_in_r_array_transp(ipoint,h1)
+            tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,i) * mos_r_in_r_array_transp(ipoint,h1)
+            tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,i) * mos_r_in_r_array_transp(ipoint,h1)
+
+            tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_r_in_r_array_transp(ipoint,i)
+            tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_r_in_r_array_transp(ipoint,i)
+            tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_r_in_r_array_transp(ipoint,i)
+          enddo
+          !$OMP END DO
+          !$OMP END PARALLEL
+
+          !$OMP PARALLEL                                                &
+          !$OMP DEFAULT (NONE)                                          &
+          !$OMP PRIVATE (p1, ipoint)                                    &
+          !$OMP SHARED (mo_num, n_points_final_grid, h1, i,             &
+          !$OMP         mos_l_in_r_array_transp, int2_grad1_u12_bimo_t, &
+          !$OMP         tmpval_1, tmpvec_1, tmp1)
+          !$OMP DO 
+          do p1 = 1, mo_num
+            do ipoint = 1, n_points_final_grid
+              tmp1(ipoint,1,p1) = mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,1) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1)
+              tmp1(ipoint,2,p1) = mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,2) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1)
+              tmp1(ipoint,3,p1) = mos_l_in_r_array_transp(ipoint,p1) * tmpvec_1(ipoint,3) + tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1)
+            enddo
+          enddo
+          !$OMP END DO
+          !$OMP END PARALLEL
+
+          call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 &
+                    , int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid         &
+                    , tmp1(1,1,1), 3*n_points_final_grid                            &
+                    , 0.d0, tmp_3d(1,1,1), mo_num*mo_num)
+
+          !$OMP PARALLEL DO PRIVATE(p1,h2,p2)
+          do p1 = 1, mo_num
+            do h2 = 1, mo_num
+              do p2 = 1, mo_num
+                no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
+              enddo
+            enddo
+          enddo
+          !$OMP END PARALLEL DO
+
+          !$OMP PARALLEL                                                &
+          !$OMP DEFAULT (NONE)                                          &
+          !$OMP PRIVATE (p2, ipoint)                                    &
+          !$OMP SHARED (mo_num, n_points_final_grid, h1, i,             &
+          !$OMP         mos_l_in_r_array_transp, int2_grad1_u12_bimo_t, &
+          !$OMP         tmpval_2, tmpvec_2, tmp1)
+          !$OMP DO 
+          do p2 = 1, mo_num
+            do ipoint = 1, n_points_final_grid
+              tmp1(ipoint,1,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,1)
+              tmp1(ipoint,2,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,2)
+              tmp1(ipoint,3,p2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p2,i) + mos_l_in_r_array_transp(ipoint,p2) * tmpvec_2(ipoint,3)
+            enddo
+          enddo
+          !$OMP END DO
+          !$OMP END PARALLEL
+
+          call dgemm( 'T', 'N', mo_num, mo_num*mo_num, 3*n_points_final_grid, 0.5d0 &
+                    , tmp1(1,1,1), 3*n_points_final_grid                            &
+                    , int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid         &
+                    , 0.d0, tmp_3d(1,1,1), mo_num)
+
+          !$OMP PARALLEL DO PRIVATE(p1,h2,p2)
+          do p1 = 1, mo_num
+            do h2 = 1, mo_num
+              do p2 = 1, mo_num
+                no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,p1,h2)
+              enddo
+            enddo
+          enddo
+          !$OMP END PARALLEL DO
+
+          ! to avoid tmp(N^4)
+          do p1 = 1, mo_num
+
+            !$OMP PARALLEL                                                  &
+            !$OMP DEFAULT (NONE)                                            &
+            !$OMP PRIVATE (ipoint)                                          &
+            !$OMP SHARED (n_points_final_grid, i, h1, p1,                   &
+            !$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         tmpval_1, tmpval_2, tmpvec_1, tmpvec_2, tmpvec_3)
+            !$OMP DO
+            do ipoint = 1, n_points_final_grid
+
+              tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) *                                          &
+                               ( int2_grad1_u12_bimo_t(ipoint,1,i,i) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) &
+                               + int2_grad1_u12_bimo_t(ipoint,2,i,i) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) &
+                               + int2_grad1_u12_bimo_t(ipoint,3,i,i) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) )
+
+              tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,p1) * mos_r_in_r_array_transp(ipoint,i)
+
+              tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_r_in_r_array_transp(ipoint,h1)
+              tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_r_in_r_array_transp(ipoint,h1)
+              tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_r_in_r_array_transp(ipoint,h1)
+
+              tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_l_in_r_array_transp(ipoint,p1)
+              tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_l_in_r_array_transp(ipoint,p1)
+              tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_l_in_r_array_transp(ipoint,p1)
+
+              tmpvec_3(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmpvec_3(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i) * mos_l_in_r_array_transp(ipoint,i)
+              tmpvec_3(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i) * mos_l_in_r_array_transp(ipoint,i)
+            enddo
+            !$OMP END DO
+            !$OMP END PARALLEL
+
+            !$OMP PARALLEL                                &
+            !$OMP DEFAULT (NONE)                          &
+            !$OMP PRIVATE (h2, ipoint)                    &
+            !$OMP SHARED (mo_num, n_points_final_grid, i, &
+            !$OMP         mos_r_in_r_array_transp,        &
+            !$OMP         int2_grad1_u12_bimo_t,          &
+            !$OMP         tmp1, tmp2, tmpval_1, tmpval_2, tmpvec_1)
+            !$OMP DO 
+            do h2 = 1, mo_num
+              do ipoint = 1, n_points_final_grid
+
+                tmp2(ipoint,h2) = mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint)     & 
+                                + int2_grad1_u12_bimo_t(ipoint,1,i,h2) * tmpvec_1(ipoint,1) &
+                                + int2_grad1_u12_bimo_t(ipoint,2,i,h2) * tmpvec_1(ipoint,2) &
+                                + int2_grad1_u12_bimo_t(ipoint,3,i,h2) * tmpvec_1(ipoint,3)
+
+                tmp1(ipoint,1,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h2)
+                tmp1(ipoint,2,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h2)
+                tmp1(ipoint,3,h2) = tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h2)
+
+              enddo
+            enddo
+            !$OMP END DO
+            !$OMP END PARALLEL
+
+            call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 &
+                      , mos_l_in_r_array_transp(1,1), n_points_final_grid    &
+                      , tmp2(1,1), n_points_final_grid                       &
+                      , 0.d0, tmp_2d(1,1), mo_num)
+
+            !$OMP PARALLEL DO PRIVATE(h2,p2)
+            do h2 = 1, mo_num
+              do p2 = 1, mo_num
+                no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
+              enddo
+            enddo
+            !$OMP END PARALLEL DO
+
+            !$OMP PARALLEL                                    &
+            !$OMP DEFAULT (NONE)                              &
+            !$OMP PRIVATE (p2, ipoint)                        &
+            !$OMP SHARED (mo_num, n_points_final_grid, i, h1, &
+            !$OMP         int2_grad1_u12_bimo_t,              &
+            !$OMP         tmpvec_2, tmpvec_3, tmp2, tmp3)
+            !$OMP DO 
+            do p2 = 1, mo_num
+              do ipoint = 1, n_points_final_grid
+
+                tmp2(ipoint,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,i) * tmpvec_2(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,p2,h1) * tmpvec_3(ipoint,1) &
+                                + int2_grad1_u12_bimo_t(ipoint,2,p2,i) * tmpvec_2(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,p2,h1) * tmpvec_3(ipoint,2) &
+                                + int2_grad1_u12_bimo_t(ipoint,3,p2,i) * tmpvec_2(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,p2,h1) * tmpvec_3(ipoint,3) 
+
+                tmp3(ipoint,1,p2) = int2_grad1_u12_bimo_t(ipoint,1,p2,h1) 
+                tmp3(ipoint,2,p2) = int2_grad1_u12_bimo_t(ipoint,2,p2,h1) 
+                tmp3(ipoint,3,p2) = int2_grad1_u12_bimo_t(ipoint,3,p2,h1) 
+              enddo
+            enddo
+            !$OMP END DO
+            !$OMP END PARALLEL
+
+            call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 &
+                      , tmp2(1,1), n_points_final_grid                       &
+                      , mos_r_in_r_array_transp(1,1), n_points_final_grid    &
+                      , 0.d0, tmp_2d(1,1), mo_num)
+
+            call dgemm( 'T', 'N', mo_num, mo_num, 3*n_points_final_grid, 0.5d0 &
+                      , tmp3(1,1,1), 3*n_points_final_grid                     &
+                      , tmp1(1,1,1), 3*n_points_final_grid                     &
+                      , 1.d0, tmp_2d(1,1), mo_num)
+
+            !$OMP PARALLEL DO PRIVATE(h2,p2)
+            do h2 = 1, mo_num
+              do p2 = 1, mo_num
+                no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
+              enddo
+            enddo
+            !$OMP END PARALLEL DO
+
+          enddo ! p1
+        enddo ! h1
+      enddo !i
+    endif
+
+    deallocate(tmp_2d, tmp_3d)
+    deallocate(tmp1, tmp2, tmp3)
+    deallocate(tmpval_1, tmpval_2)
+    deallocate(tmpvec_1, tmpvec_2, tmpvec_3)
+
+    no_aaa_contraction = -0.5d0 * no_aaa_contraction
+
+    !$OMP PARALLEL                 &
+    !$OMP DEFAULT (NONE)           &
+    !$OMP PRIVATE (h1, h2, p1, p2) & 
+    !$OMP SHARED (no_aaa_contraction, mo_num)
+
+    !$OMP DO 
+    do h1 = 1, mo_num
+      do h2 = 1, mo_num
+        do p1 = 1, mo_num
+          do p2 = p1, mo_num
+            no_aaa_contraction(p2,h2,p1,h1) -= no_aaa_contraction(p1,h2,p2,h1)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+
+    !$OMP DO 
+    do h1 = 1, mo_num
+      do h2 = 1, mo_num
+        do p1 = 2, mo_num
+          do p2 = 1, p1-1
+            no_aaa_contraction(p2,h2,p1,h1) = -no_aaa_contraction(p1,h2,p2,h1)
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+
+    !$OMP DO 
+    do h1 = 1, mo_num-1
+      do h2 = h1+1, mo_num
+        do p1 = 2, mo_num
+          do p2 = 1, p1-1
+            no_aaa_contraction(p2,h2,p1,h1) *= -1.d0
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END PARALLEL
+
+  endif
+
+  call wall_time(wall1)
+  print*,' Wall time for no_aaa_contraction', wall1-wall0
+
+END_PROVIDER
+
+! ---

src/tc_bi_ortho/normal_ordered_old.irp.f

@@ -89,6 +89,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_old, (mo_num, mo_num,
                hthree_aaa = 0.d0
              endif
             endif
+
             normal_two_body_bi_orth_old(p2,h2,p1,h1) = 0.5d0*(hthree_aba + hthree_aab + hthree_aaa)
           enddo
         enddo
@@ -350,7 +351,8 @@ subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
     call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
     int_exc_23 = -1.d0 * integral
 
-    hthree +=  1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
+    !hthree +=  1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
+    hthree +=  0.5d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
   enddo
 
   return